]> gitweb.michael.orlitzky.com - sage.d.git/blob - mjo/eja/eja_element.py
makefile: add --long to the test command.
[sage.d.git] / mjo / eja / eja_element.py
1 from sage.matrix.constructor import matrix
2 from sage.modules.free_module import VectorSpace
3 from sage.modules.with_basis.indexed_element import IndexedFreeModuleElement
4
5 from mjo.eja.eja_operator import FiniteDimensionalEuclideanJordanAlgebraOperator
6 from mjo.eja.eja_utils import _mat2vec
7
8 class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
9 """
10 An element of a Euclidean Jordan algebra.
11 """
12
13 def __dir__(self):
14 """
15 Oh man, I should not be doing this. This hides the "disabled"
16 methods ``left_matrix`` and ``matrix`` from introspection;
17 in particular it removes them from tab-completion.
18 """
19 return filter(lambda s: s not in ['left_matrix', 'matrix'],
20 dir(self.__class__) )
21
22
23
24
25 def __pow__(self, n):
26 """
27 Return ``self`` raised to the power ``n``.
28
29 Jordan algebras are always power-associative; see for
30 example Faraut and Korányi, Proposition II.1.2 (ii).
31
32 We have to override this because our superclass uses row
33 vectors instead of column vectors! We, on the other hand,
34 assume column vectors everywhere.
35
36 SETUP::
37
38 sage: from mjo.eja.eja_algebra import random_eja
39
40 TESTS:
41
42 The definition of `x^2` is the unambiguous `x*x`::
43
44 sage: set_random_seed()
45 sage: x = random_eja().random_element()
46 sage: x*x == (x^2)
47 True
48
49 A few examples of power-associativity::
50
51 sage: set_random_seed()
52 sage: x = random_eja().random_element()
53 sage: x*(x*x)*(x*x) == x^5
54 True
55 sage: (x*x)*(x*x*x) == x^5
56 True
57
58 We also know that powers operator-commute (Koecher, Chapter
59 III, Corollary 1)::
60
61 sage: set_random_seed()
62 sage: x = random_eja().random_element()
63 sage: m = ZZ.random_element(0,10)
64 sage: n = ZZ.random_element(0,10)
65 sage: Lxm = (x^m).operator()
66 sage: Lxn = (x^n).operator()
67 sage: Lxm*Lxn == Lxn*Lxm
68 True
69
70 """
71 if n == 0:
72 return self.parent().one()
73 elif n == 1:
74 return self
75 else:
76 return (self**(n-1))*self
77
78
79 def apply_univariate_polynomial(self, p):
80 """
81 Apply the univariate polynomial ``p`` to this element.
82
83 A priori, SageMath won't allow us to apply a univariate
84 polynomial to an element of an EJA, because we don't know
85 that EJAs are rings (they are usually not associative). Of
86 course, we know that EJAs are power-associative, so the
87 operation is ultimately kosher. This function sidesteps
88 the CAS to get the answer we want and expect.
89
90 SETUP::
91
92 sage: from mjo.eja.eja_algebra import (HadamardEJA,
93 ....: random_eja)
94
95 EXAMPLES::
96
97 sage: R = PolynomialRing(QQ, 't')
98 sage: t = R.gen(0)
99 sage: p = t^4 - t^3 + 5*t - 2
100 sage: J = HadamardEJA(5)
101 sage: J.one().apply_univariate_polynomial(p) == 3*J.one()
102 True
103
104 TESTS:
105
106 We should always get back an element of the algebra::
107
108 sage: set_random_seed()
109 sage: p = PolynomialRing(AA, 't').random_element()
110 sage: J = random_eja()
111 sage: x = J.random_element()
112 sage: x.apply_univariate_polynomial(p) in J
113 True
114
115 """
116 if len(p.variables()) > 1:
117 raise ValueError("not a univariate polynomial")
118 P = self.parent()
119 R = P.base_ring()
120 # Convert the coeficcients to the parent's base ring,
121 # because a priori they might live in an (unnecessarily)
122 # larger ring for which P.sum() would fail below.
123 cs = [ R(c) for c in p.coefficients(sparse=False) ]
124 return P.sum( cs[k]*(self**k) for k in range(len(cs)) )
125
126
127 def characteristic_polynomial(self):
128 """
129 Return the characteristic polynomial of this element.
130
131 SETUP::
132
133 sage: from mjo.eja.eja_algebra import HadamardEJA
134
135 EXAMPLES:
136
137 The rank of `R^3` is three, and the minimal polynomial of
138 the identity element is `(t-1)` from which it follows that
139 the characteristic polynomial should be `(t-1)^3`::
140
141 sage: J = HadamardEJA(3)
142 sage: J.one().characteristic_polynomial()
143 t^3 - 3*t^2 + 3*t - 1
144
145 Likewise, the characteristic of the zero element in the
146 rank-three algebra `R^{n}` should be `t^{3}`::
147
148 sage: J = HadamardEJA(3)
149 sage: J.zero().characteristic_polynomial()
150 t^3
151
152 TESTS:
153
154 The characteristic polynomial of an element should evaluate
155 to zero on that element::
156
157 sage: set_random_seed()
158 sage: x = HadamardEJA(3).random_element()
159 sage: p = x.characteristic_polynomial()
160 sage: x.apply_univariate_polynomial(p)
161 0
162
163 The characteristic polynomials of the zero and unit elements
164 should be what we think they are in a subalgebra, too::
165
166 sage: J = HadamardEJA(3)
167 sage: p1 = J.one().characteristic_polynomial()
168 sage: q1 = J.zero().characteristic_polynomial()
169 sage: e0,e1,e2 = J.gens()
170 sage: A = (e0 + 2*e1 + 3*e2).subalgebra_generated_by() # dim 3
171 sage: p2 = A.one().characteristic_polynomial()
172 sage: q2 = A.zero().characteristic_polynomial()
173 sage: p1 == p2
174 True
175 sage: q1 == q2
176 True
177
178 """
179 p = self.parent().characteristic_polynomial_of()
180 return p(*self.to_vector())
181
182
183 def inner_product(self, other):
184 """
185 Return the parent algebra's inner product of myself and ``other``.
186
187 SETUP::
188
189 sage: from mjo.eja.eja_algebra import (
190 ....: ComplexHermitianEJA,
191 ....: JordanSpinEJA,
192 ....: QuaternionHermitianEJA,
193 ....: RealSymmetricEJA,
194 ....: random_eja)
195
196 EXAMPLES:
197
198 The inner product in the Jordan spin algebra is the usual
199 inner product on `R^n` (this example only works because the
200 basis for the Jordan algebra is the standard basis in `R^n`)::
201
202 sage: J = JordanSpinEJA(3)
203 sage: x = vector(QQ,[1,2,3])
204 sage: y = vector(QQ,[4,5,6])
205 sage: x.inner_product(y)
206 32
207 sage: J.from_vector(x).inner_product(J.from_vector(y))
208 32
209
210 The inner product on `S^n` is `<X,Y> = trace(X*Y)`, where
211 multiplication is the usual matrix multiplication in `S^n`,
212 so the inner product of the identity matrix with itself
213 should be the `n`::
214
215 sage: J = RealSymmetricEJA(3)
216 sage: J.one().inner_product(J.one())
217 3
218
219 Likewise, the inner product on `C^n` is `<X,Y> =
220 Re(trace(X*Y))`, where we must necessarily take the real
221 part because the product of Hermitian matrices may not be
222 Hermitian::
223
224 sage: J = ComplexHermitianEJA(3)
225 sage: J.one().inner_product(J.one())
226 3
227
228 Ditto for the quaternions::
229
230 sage: J = QuaternionHermitianEJA(3)
231 sage: J.one().inner_product(J.one())
232 3
233
234 TESTS:
235
236 Ensure that we can always compute an inner product, and that
237 it gives us back a real number::
238
239 sage: set_random_seed()
240 sage: J = random_eja()
241 sage: x,y = J.random_elements(2)
242 sage: x.inner_product(y) in RLF
243 True
244
245 """
246 P = self.parent()
247 if not other in P:
248 raise TypeError("'other' must live in the same algebra")
249
250 return P.inner_product(self, other)
251
252
253 def operator_commutes_with(self, other):
254 """
255 Return whether or not this element operator-commutes
256 with ``other``.
257
258 SETUP::
259
260 sage: from mjo.eja.eja_algebra import random_eja
261
262 EXAMPLES:
263
264 The definition of a Jordan algebra says that any element
265 operator-commutes with its square::
266
267 sage: set_random_seed()
268 sage: x = random_eja().random_element()
269 sage: x.operator_commutes_with(x^2)
270 True
271
272 TESTS:
273
274 Test Lemma 1 from Chapter III of Koecher::
275
276 sage: set_random_seed()
277 sage: u,v = random_eja().random_elements(2)
278 sage: lhs = u.operator_commutes_with(u*v)
279 sage: rhs = v.operator_commutes_with(u^2)
280 sage: lhs == rhs
281 True
282
283 Test the first polarization identity from my notes, Koecher
284 Chapter III, or from Baes (2.3)::
285
286 sage: set_random_seed()
287 sage: x,y = random_eja().random_elements(2)
288 sage: Lx = x.operator()
289 sage: Ly = y.operator()
290 sage: Lxx = (x*x).operator()
291 sage: Lxy = (x*y).operator()
292 sage: bool(2*Lx*Lxy + Ly*Lxx == 2*Lxy*Lx + Lxx*Ly)
293 True
294
295 Test the second polarization identity from my notes or from
296 Baes (2.4)::
297
298 sage: set_random_seed()
299 sage: x,y,z = random_eja().random_elements(3)
300 sage: Lx = x.operator()
301 sage: Ly = y.operator()
302 sage: Lz = z.operator()
303 sage: Lzy = (z*y).operator()
304 sage: Lxy = (x*y).operator()
305 sage: Lxz = (x*z).operator()
306 sage: bool(Lx*Lzy + Lz*Lxy + Ly*Lxz == Lzy*Lx + Lxy*Lz + Lxz*Ly)
307 True
308
309 Test the third polarization identity from my notes or from
310 Baes (2.5)::
311
312 sage: set_random_seed()
313 sage: u,y,z = random_eja().random_elements(3)
314 sage: Lu = u.operator()
315 sage: Ly = y.operator()
316 sage: Lz = z.operator()
317 sage: Lzy = (z*y).operator()
318 sage: Luy = (u*y).operator()
319 sage: Luz = (u*z).operator()
320 sage: Luyz = (u*(y*z)).operator()
321 sage: lhs = Lu*Lzy + Lz*Luy + Ly*Luz
322 sage: rhs = Luyz + Ly*Lu*Lz + Lz*Lu*Ly
323 sage: bool(lhs == rhs)
324 True
325
326 """
327 if not other in self.parent():
328 raise TypeError("'other' must live in the same algebra")
329
330 A = self.operator()
331 B = other.operator()
332 return (A*B == B*A)
333
334
335 def det(self):
336 """
337 Return my determinant, the product of my eigenvalues.
338
339 SETUP::
340
341 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
342 ....: TrivialEJA,
343 ....: random_eja)
344
345 EXAMPLES::
346
347 sage: J = JordanSpinEJA(2)
348 sage: e0,e1 = J.gens()
349 sage: x = sum( J.gens() )
350 sage: x.det()
351 0
352
353 ::
354
355 sage: J = JordanSpinEJA(3)
356 sage: e0,e1,e2 = J.gens()
357 sage: x = sum( J.gens() )
358 sage: x.det()
359 -1
360
361 The determinant of the sole element in the rank-zero trivial
362 algebra is ``1``, by three paths of reasoning. First, its
363 characteristic polynomial is a constant ``1``, so the constant
364 term in that polynomial is ``1``. Second, the characteristic
365 polynomial evaluated at zero is again ``1``. And finally, the
366 (empty) product of its eigenvalues is likewise just unity::
367
368 sage: J = TrivialEJA()
369 sage: J.zero().det()
370 1
371
372 TESTS:
373
374 An element is invertible if and only if its determinant is
375 non-zero::
376
377 sage: set_random_seed()
378 sage: x = random_eja().random_element()
379 sage: x.is_invertible() == (x.det() != 0)
380 True
381
382 Ensure that the determinant is multiplicative on an associative
383 subalgebra as in Faraut and Korányi's Proposition II.2.2::
384
385 sage: set_random_seed()
386 sage: J = random_eja().random_element().subalgebra_generated_by()
387 sage: x,y = J.random_elements(2)
388 sage: (x*y).det() == x.det()*y.det()
389 True
390 """
391 P = self.parent()
392 r = P.rank()
393
394 if r == 0:
395 # Special case, since we don't get the a0=1
396 # coefficient when the rank of the algebra
397 # is zero.
398 return P.base_ring().one()
399
400 p = P._charpoly_coefficients()[0]
401 # The _charpoly_coeff function already adds the factor of -1
402 # to ensure that _charpoly_coefficients()[0] is really what
403 # appears in front of t^{0} in the charpoly. However, we want
404 # (-1)^r times THAT for the determinant.
405 return ((-1)**r)*p(*self.to_vector())
406
407
408 def inverse(self):
409 """
410 Return the Jordan-multiplicative inverse of this element.
411
412 ALGORITHM:
413
414 We appeal to the quadratic representation as in Koecher's
415 Theorem 12 in Chapter III, Section 5.
416
417 SETUP::
418
419 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
420 ....: JordanSpinEJA,
421 ....: random_eja)
422
423 EXAMPLES:
424
425 The inverse in the spin factor algebra is given in Alizadeh's
426 Example 11.11::
427
428 sage: set_random_seed()
429 sage: J = JordanSpinEJA.random_instance()
430 sage: x = J.random_element()
431 sage: while not x.is_invertible():
432 ....: x = J.random_element()
433 sage: x_vec = x.to_vector()
434 sage: x0 = x_vec[:1]
435 sage: x_bar = x_vec[1:]
436 sage: coeff = x0.inner_product(x0) - x_bar.inner_product(x_bar)
437 sage: x_inverse = x_vec.parent()(x0.list() + (-x_bar).list())
438 sage: if not coeff.is_zero(): x_inverse = x_inverse/coeff
439 sage: x.inverse() == J.from_vector(x_inverse)
440 True
441
442 Trying to invert a non-invertible element throws an error:
443
444 sage: JordanSpinEJA(3).zero().inverse()
445 Traceback (most recent call last):
446 ...
447 ValueError: element is not invertible
448
449 TESTS:
450
451 The identity element is its own inverse::
452
453 sage: set_random_seed()
454 sage: J = random_eja()
455 sage: J.one().inverse() == J.one()
456 True
457
458 If an element has an inverse, it acts like one::
459
460 sage: set_random_seed()
461 sage: J = random_eja()
462 sage: x = J.random_element()
463 sage: (not x.is_invertible()) or (x.inverse()*x == J.one())
464 True
465
466 The inverse of the inverse is what we started with::
467
468 sage: set_random_seed()
469 sage: J = random_eja()
470 sage: x = J.random_element()
471 sage: (not x.is_invertible()) or (x.inverse().inverse() == x)
472 True
473
474 Proposition II.2.3 in Faraut and Korányi says that the inverse
475 of an element is the inverse of its left-multiplication operator
476 applied to the algebra's identity, when that inverse exists::
477
478 sage: set_random_seed()
479 sage: J = random_eja()
480 sage: x = J.random_element()
481 sage: (not x.operator().is_invertible()) or (
482 ....: x.operator().inverse()(J.one()) == x.inverse() )
483 True
484
485 Proposition II.2.4 in Faraut and Korányi gives a formula for
486 the inverse based on the characteristic polynomial and the
487 Cayley-Hamilton theorem for Euclidean Jordan algebras::
488
489 sage: set_random_seed()
490 sage: J = ComplexHermitianEJA(3)
491 sage: x = J.random_element()
492 sage: while not x.is_invertible():
493 ....: x = J.random_element()
494 sage: r = J.rank()
495 sage: a = x.characteristic_polynomial().coefficients(sparse=False)
496 sage: expected = (-1)^(r+1)/x.det()
497 sage: expected *= sum( a[i+1]*x^i for i in range(r) )
498 sage: x.inverse() == expected
499 True
500
501 """
502 if not self.is_invertible():
503 raise ValueError("element is not invertible")
504
505 if self.parent()._charpoly_coefficients.is_in_cache():
506 # We can invert using our charpoly if it will be fast to
507 # compute. If the coefficients are cached, our rank had
508 # better be too!
509 r = self.parent().rank()
510 a = self.characteristic_polynomial().coefficients(sparse=False)
511 return (-1)**(r+1)*sum(a[i+1]*self**i for i in range(r))/self.det()
512
513 return (~self.quadratic_representation())(self)
514
515
516 def is_invertible(self):
517 """
518 Return whether or not this element is invertible.
519
520 ALGORITHM:
521
522 The usual way to do this is to check if the determinant is
523 zero, but we need the characteristic polynomial for the
524 determinant. The minimal polynomial is a lot easier to get,
525 so we use Corollary 2 in Chapter V of Koecher to check
526 whether or not the paren't algebra's zero element is a root
527 of this element's minimal polynomial.
528
529 That is... unless the coefficients of our algebra's
530 "characteristic polynomial of" function are already cached!
531 In that case, we just use the determinant (which will be fast
532 as a result).
533
534 Beware that we can't use the superclass method, because it
535 relies on the algebra being associative.
536
537 SETUP::
538
539 sage: from mjo.eja.eja_algebra import random_eja
540
541 TESTS:
542
543 The identity element is always invertible::
544
545 sage: set_random_seed()
546 sage: J = random_eja()
547 sage: J.one().is_invertible()
548 True
549
550 The zero element is never invertible in a non-trivial algebra::
551
552 sage: set_random_seed()
553 sage: J = random_eja()
554 sage: (not J.is_trivial()) and J.zero().is_invertible()
555 False
556
557 """
558 if self.is_zero():
559 if self.parent().is_trivial():
560 return True
561 else:
562 return False
563
564 if self.parent()._charpoly_coefficients.is_in_cache():
565 # The determinant will be quicker than computing the minimal
566 # polynomial from scratch, most likely.
567 return (not self.det().is_zero())
568
569 # In fact, we only need to know if the constant term is non-zero,
570 # so we can pass in the field's zero element instead.
571 zero = self.base_ring().zero()
572 p = self.minimal_polynomial()
573 return not (p(zero) == zero)
574
575
576 def is_primitive_idempotent(self):
577 """
578 Return whether or not this element is a primitive (or minimal)
579 idempotent.
580
581 A primitive idempotent is a non-zero idempotent that is not
582 the sum of two other non-zero idempotents. Remark 2.7.15 in
583 Baes shows that this is what he refers to as a "minimal
584 idempotent."
585
586 An element of a Euclidean Jordan algebra is a minimal idempotent
587 if it :meth:`is_idempotent` and if its Peirce subalgebra
588 corresponding to the eigenvalue ``1`` has dimension ``1`` (Baes,
589 Proposition 2.7.17).
590
591 SETUP::
592
593 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
594 ....: RealSymmetricEJA,
595 ....: TrivialEJA,
596 ....: random_eja)
597
598 WARNING::
599
600 This method is sloooooow.
601
602 EXAMPLES:
603
604 The spectral decomposition of a non-regular element should always
605 contain at least one non-minimal idempotent::
606
607 sage: J = RealSymmetricEJA(3)
608 sage: x = sum(J.gens())
609 sage: x.is_regular()
610 False
611 sage: [ c.is_primitive_idempotent()
612 ....: for (l,c) in x.spectral_decomposition() ]
613 [False, True]
614
615 On the other hand, the spectral decomposition of a regular
616 element should always be in terms of minimal idempotents::
617
618 sage: J = JordanSpinEJA(4)
619 sage: x = sum( i*J.gens()[i] for i in range(len(J.gens())) )
620 sage: x.is_regular()
621 True
622 sage: [ c.is_primitive_idempotent()
623 ....: for (l,c) in x.spectral_decomposition() ]
624 [True, True]
625
626 TESTS:
627
628 The identity element is minimal only in an EJA of rank one::
629
630 sage: set_random_seed()
631 sage: J = random_eja()
632 sage: J.rank() == 1 or not J.one().is_primitive_idempotent()
633 True
634
635 A non-idempotent cannot be a minimal idempotent::
636
637 sage: set_random_seed()
638 sage: J = JordanSpinEJA(4)
639 sage: x = J.random_element()
640 sage: (not x.is_idempotent()) and x.is_primitive_idempotent()
641 False
642
643 Proposition 2.7.19 in Baes says that an element is a minimal
644 idempotent if and only if it's idempotent with trace equal to
645 unity::
646
647 sage: set_random_seed()
648 sage: J = JordanSpinEJA(4)
649 sage: x = J.random_element()
650 sage: expected = (x.is_idempotent() and x.trace() == 1)
651 sage: actual = x.is_primitive_idempotent()
652 sage: actual == expected
653 True
654
655 Primitive idempotents must be non-zero::
656
657 sage: set_random_seed()
658 sage: J = random_eja()
659 sage: J.zero().is_idempotent()
660 True
661 sage: J.zero().is_primitive_idempotent()
662 False
663
664 As a consequence of the fact that primitive idempotents must
665 be non-zero, there are no primitive idempotents in a trivial
666 Euclidean Jordan algebra::
667
668 sage: J = TrivialEJA()
669 sage: J.one().is_idempotent()
670 True
671 sage: J.one().is_primitive_idempotent()
672 False
673
674 """
675 if not self.is_idempotent():
676 return False
677
678 if self.is_zero():
679 return False
680
681 (_,_,J1) = self.parent().peirce_decomposition(self)
682 return (J1.dimension() == 1)
683
684
685 def is_nilpotent(self):
686 """
687 Return whether or not some power of this element is zero.
688
689 ALGORITHM:
690
691 We use Theorem 5 in Chapter III of Koecher, which says that
692 an element ``x`` is nilpotent if and only if ``x.operator()``
693 is nilpotent. And it is a basic fact of linear algebra that
694 an operator on an `n`-dimensional space is nilpotent if and
695 only if, when raised to the `n`th power, it equals the zero
696 operator (for example, see Axler Corollary 8.8).
697
698 SETUP::
699
700 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
701 ....: random_eja)
702
703 EXAMPLES::
704
705 sage: J = JordanSpinEJA(3)
706 sage: x = sum(J.gens())
707 sage: x.is_nilpotent()
708 False
709
710 TESTS:
711
712 The identity element is never nilpotent, except in a trivial EJA::
713
714 sage: set_random_seed()
715 sage: J = random_eja()
716 sage: J.one().is_nilpotent() and not J.is_trivial()
717 False
718
719 The additive identity is always nilpotent::
720
721 sage: set_random_seed()
722 sage: random_eja().zero().is_nilpotent()
723 True
724
725 """
726 P = self.parent()
727 zero_operator = P.zero().operator()
728 return self.operator()**P.dimension() == zero_operator
729
730
731 def is_regular(self):
732 """
733 Return whether or not this is a regular element.
734
735 SETUP::
736
737 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
738 ....: random_eja)
739
740 EXAMPLES:
741
742 The identity element always has degree one, but any element
743 linearly-independent from it is regular::
744
745 sage: J = JordanSpinEJA(5)
746 sage: J.one().is_regular()
747 False
748 sage: e0, e1, e2, e3, e4 = J.gens() # e0 is the identity
749 sage: for x in J.gens():
750 ....: (J.one() + x).is_regular()
751 False
752 True
753 True
754 True
755 True
756
757 TESTS:
758
759 The zero element should never be regular, unless the parent
760 algebra has dimension less than or equal to one::
761
762 sage: set_random_seed()
763 sage: J = random_eja()
764 sage: J.dimension() <= 1 or not J.zero().is_regular()
765 True
766
767 The unit element isn't regular unless the algebra happens to
768 consist of only its scalar multiples::
769
770 sage: set_random_seed()
771 sage: J = random_eja()
772 sage: J.dimension() <= 1 or not J.one().is_regular()
773 True
774
775 """
776 return self.degree() == self.parent().rank()
777
778
779 def degree(self):
780 """
781 Return the degree of this element, which is defined to be
782 the degree of its minimal polynomial.
783
784 ALGORITHM:
785
786 For now, we skip the messy minimal polynomial computation
787 and instead return the dimension of the vector space spanned
788 by the powers of this element. The latter is a bit more
789 straightforward to compute.
790
791 SETUP::
792
793 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
794 ....: random_eja)
795
796 EXAMPLES::
797
798 sage: J = JordanSpinEJA(4)
799 sage: J.one().degree()
800 1
801 sage: e0,e1,e2,e3 = J.gens()
802 sage: (e0 - e1).degree()
803 2
804
805 In the spin factor algebra (of rank two), all elements that
806 aren't multiples of the identity are regular::
807
808 sage: set_random_seed()
809 sage: J = JordanSpinEJA.random_instance()
810 sage: n = J.dimension()
811 sage: x = J.random_element()
812 sage: x.degree() == min(n,2) or (x == x.coefficient(0)*J.one())
813 True
814
815 TESTS:
816
817 The zero and unit elements are both of degree one in nontrivial
818 algebras::
819
820 sage: set_random_seed()
821 sage: J = random_eja()
822 sage: d = J.zero().degree()
823 sage: (J.is_trivial() and d == 0) or d == 1
824 True
825 sage: d = J.one().degree()
826 sage: (J.is_trivial() and d == 0) or d == 1
827 True
828
829 Our implementation agrees with the definition::
830
831 sage: set_random_seed()
832 sage: x = random_eja().random_element()
833 sage: x.degree() == x.minimal_polynomial().degree()
834 True
835
836 """
837 if self.is_zero() and not self.parent().is_trivial():
838 # The minimal polynomial of zero in a nontrivial algebra
839 # is "t"; in a trivial algebra it's "1" by convention
840 # (it's an empty product).
841 return 1
842 return self.subalgebra_generated_by().dimension()
843
844
845 def left_matrix(self):
846 """
847 Our parent class defines ``left_matrix`` and ``matrix``
848 methods whose names are misleading. We don't want them.
849 """
850 raise NotImplementedError("use operator().matrix() instead")
851
852 matrix = left_matrix
853
854
855 def minimal_polynomial(self):
856 """
857 Return the minimal polynomial of this element,
858 as a function of the variable `t`.
859
860 ALGORITHM:
861
862 We restrict ourselves to the associative subalgebra
863 generated by this element, and then return the minimal
864 polynomial of this element's operator matrix (in that
865 subalgebra). This works by Baes Proposition 2.3.16.
866
867 SETUP::
868
869 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
870 ....: RealSymmetricEJA,
871 ....: TrivialEJA,
872 ....: random_eja)
873
874 EXAMPLES:
875
876 Keeping in mind that the polynomial ``1`` evaluates the identity
877 element (also the zero element) of the trivial algebra, it is clear
878 that the polynomial ``1`` is the minimal polynomial of the only
879 element in a trivial algebra::
880
881 sage: J = TrivialEJA()
882 sage: J.one().minimal_polynomial()
883 1
884 sage: J.zero().minimal_polynomial()
885 1
886
887 TESTS:
888
889 The minimal polynomial of the identity and zero elements are
890 always the same, except in trivial algebras where the minimal
891 polynomial of the unit/zero element is ``1``::
892
893 sage: set_random_seed()
894 sage: J = random_eja()
895 sage: mu = J.one().minimal_polynomial()
896 sage: t = mu.parent().gen()
897 sage: mu + int(J.is_trivial())*(t-2)
898 t - 1
899 sage: mu = J.zero().minimal_polynomial()
900 sage: t = mu.parent().gen()
901 sage: mu + int(J.is_trivial())*(t-1)
902 t
903
904 The degree of an element is (by one definition) the degree
905 of its minimal polynomial::
906
907 sage: set_random_seed()
908 sage: x = random_eja().random_element()
909 sage: x.degree() == x.minimal_polynomial().degree()
910 True
911
912 The minimal polynomial and the characteristic polynomial coincide
913 and are known (see Alizadeh, Example 11.11) for all elements of
914 the spin factor algebra that aren't scalar multiples of the
915 identity. We require the dimension of the algebra to be at least
916 two here so that said elements actually exist::
917
918 sage: set_random_seed()
919 sage: n_max = max(2, JordanSpinEJA._max_random_instance_size())
920 sage: n = ZZ.random_element(2, n_max)
921 sage: J = JordanSpinEJA(n)
922 sage: y = J.random_element()
923 sage: while y == y.coefficient(0)*J.one():
924 ....: y = J.random_element()
925 sage: y0 = y.to_vector()[0]
926 sage: y_bar = y.to_vector()[1:]
927 sage: actual = y.minimal_polynomial()
928 sage: t = PolynomialRing(J.base_ring(),'t').gen(0)
929 sage: expected = t^2 - 2*y0*t + (y0^2 - norm(y_bar)^2)
930 sage: bool(actual == expected)
931 True
932
933 The minimal polynomial should always kill its element::
934
935 sage: set_random_seed()
936 sage: x = random_eja().random_element()
937 sage: p = x.minimal_polynomial()
938 sage: x.apply_univariate_polynomial(p)
939 0
940
941 The minimal polynomial is invariant under a change of basis,
942 and in particular, a re-scaling of the basis::
943
944 sage: set_random_seed()
945 sage: n_max = RealSymmetricEJA._max_random_instance_size()
946 sage: n = ZZ.random_element(1, n_max)
947 sage: J1 = RealSymmetricEJA(n)
948 sage: J2 = RealSymmetricEJA(n,normalize_basis=False)
949 sage: X = random_matrix(AA,n)
950 sage: X = X*X.transpose()
951 sage: x1 = J1(X)
952 sage: x2 = J2(X)
953 sage: x1.minimal_polynomial() == x2.minimal_polynomial()
954 True
955
956 """
957 if self.is_zero():
958 # We would generate a zero-dimensional subalgebra
959 # where the minimal polynomial would be constant.
960 # That might be correct, but only if *this* algebra
961 # is trivial too.
962 if not self.parent().is_trivial():
963 # Pretty sure we know what the minimal polynomial of
964 # the zero operator is going to be. This ensures
965 # consistency of e.g. the polynomial variable returned
966 # in the "normal" case without us having to think about it.
967 return self.operator().minimal_polynomial()
968
969 A = self.subalgebra_generated_by(orthonormalize_basis=False)
970 return A(self).operator().minimal_polynomial()
971
972
973
974 def to_matrix(self):
975 """
976 Return an (often more natural) representation of this element as a
977 matrix.
978
979 Every finite-dimensional Euclidean Jordan Algebra is a direct
980 sum of five simple algebras, four of which comprise Hermitian
981 matrices. This method returns a "natural" matrix
982 representation of this element as either a Hermitian matrix or
983 column vector.
984
985 SETUP::
986
987 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
988 ....: QuaternionHermitianEJA)
989
990 EXAMPLES::
991
992 sage: J = ComplexHermitianEJA(3)
993 sage: J.one()
994 e0 + e3 + e8
995 sage: J.one().to_matrix()
996 [1 0 0 0 0 0]
997 [0 1 0 0 0 0]
998 [0 0 1 0 0 0]
999 [0 0 0 1 0 0]
1000 [0 0 0 0 1 0]
1001 [0 0 0 0 0 1]
1002
1003 ::
1004
1005 sage: J = QuaternionHermitianEJA(3)
1006 sage: J.one()
1007 e0 + e5 + e14
1008 sage: J.one().to_matrix()
1009 [1 0 0 0 0 0 0 0 0 0 0 0]
1010 [0 1 0 0 0 0 0 0 0 0 0 0]
1011 [0 0 1 0 0 0 0 0 0 0 0 0]
1012 [0 0 0 1 0 0 0 0 0 0 0 0]
1013 [0 0 0 0 1 0 0 0 0 0 0 0]
1014 [0 0 0 0 0 1 0 0 0 0 0 0]
1015 [0 0 0 0 0 0 1 0 0 0 0 0]
1016 [0 0 0 0 0 0 0 1 0 0 0 0]
1017 [0 0 0 0 0 0 0 0 1 0 0 0]
1018 [0 0 0 0 0 0 0 0 0 1 0 0]
1019 [0 0 0 0 0 0 0 0 0 0 1 0]
1020 [0 0 0 0 0 0 0 0 0 0 0 1]
1021 """
1022 B = self.parent().matrix_basis()
1023 W = self.parent().matrix_space()
1024
1025 # This is just a manual "from_vector()", but of course
1026 # matrix spaces aren't vector spaces in sage, so they
1027 # don't have a from_vector() method.
1028 return W.linear_combination( zip(B, self.to_vector()) )
1029
1030
1031 def norm(self):
1032 """
1033 The norm of this element with respect to :meth:`inner_product`.
1034
1035 SETUP::
1036
1037 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1038 ....: HadamardEJA)
1039
1040 EXAMPLES::
1041
1042 sage: J = HadamardEJA(2)
1043 sage: x = sum(J.gens())
1044 sage: x.norm()
1045 1.414213562373095?
1046
1047 ::
1048
1049 sage: J = JordanSpinEJA(4)
1050 sage: x = sum(J.gens())
1051 sage: x.norm()
1052 2
1053
1054 """
1055 return self.inner_product(self).sqrt()
1056
1057
1058 def operator(self):
1059 """
1060 Return the left-multiplication-by-this-element
1061 operator on the ambient algebra.
1062
1063 SETUP::
1064
1065 sage: from mjo.eja.eja_algebra import random_eja
1066
1067 TESTS::
1068
1069 sage: set_random_seed()
1070 sage: J = random_eja()
1071 sage: x,y = J.random_elements(2)
1072 sage: x.operator()(y) == x*y
1073 True
1074 sage: y.operator()(x) == x*y
1075 True
1076
1077 """
1078 P = self.parent()
1079 left_mult_by_self = lambda y: self*y
1080 L = P.module_morphism(function=left_mult_by_self, codomain=P)
1081 return FiniteDimensionalEuclideanJordanAlgebraOperator(
1082 P,
1083 P,
1084 L.matrix() )
1085
1086
1087 def quadratic_representation(self, other=None):
1088 """
1089 Return the quadratic representation of this element.
1090
1091 SETUP::
1092
1093 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1094 ....: random_eja)
1095
1096 EXAMPLES:
1097
1098 The explicit form in the spin factor algebra is given by
1099 Alizadeh's Example 11.12::
1100
1101 sage: set_random_seed()
1102 sage: x = JordanSpinEJA.random_instance().random_element()
1103 sage: x_vec = x.to_vector()
1104 sage: Q = matrix.identity(x.base_ring(), 0)
1105 sage: n = x_vec.degree()
1106 sage: if n > 0:
1107 ....: x0 = x_vec[0]
1108 ....: x_bar = x_vec[1:]
1109 ....: A = matrix(x.base_ring(), 1, [x_vec.inner_product(x_vec)])
1110 ....: B = 2*x0*x_bar.row()
1111 ....: C = 2*x0*x_bar.column()
1112 ....: D = matrix.identity(x.base_ring(), n-1)
1113 ....: D = (x0^2 - x_bar.inner_product(x_bar))*D
1114 ....: D = D + 2*x_bar.tensor_product(x_bar)
1115 ....: Q = matrix.block(2,2,[A,B,C,D])
1116 sage: Q == x.quadratic_representation().matrix()
1117 True
1118
1119 Test all of the properties from Theorem 11.2 in Alizadeh::
1120
1121 sage: set_random_seed()
1122 sage: J = random_eja()
1123 sage: x,y = J.random_elements(2)
1124 sage: Lx = x.operator()
1125 sage: Lxx = (x*x).operator()
1126 sage: Qx = x.quadratic_representation()
1127 sage: Qy = y.quadratic_representation()
1128 sage: Qxy = x.quadratic_representation(y)
1129 sage: Qex = J.one().quadratic_representation(x)
1130 sage: n = ZZ.random_element(10)
1131 sage: Qxn = (x^n).quadratic_representation()
1132
1133 Property 1:
1134
1135 sage: 2*Qxy == (x+y).quadratic_representation() - Qx - Qy
1136 True
1137
1138 Property 2 (multiply on the right for :trac:`28272`):
1139
1140 sage: alpha = J.base_ring().random_element()
1141 sage: (alpha*x).quadratic_representation() == Qx*(alpha^2)
1142 True
1143
1144 Property 3:
1145
1146 sage: not x.is_invertible() or ( Qx(x.inverse()) == x )
1147 True
1148
1149 sage: not x.is_invertible() or (
1150 ....: ~Qx
1151 ....: ==
1152 ....: x.inverse().quadratic_representation() )
1153 True
1154
1155 sage: Qxy(J.one()) == x*y
1156 True
1157
1158 Property 4:
1159
1160 sage: not x.is_invertible() or (
1161 ....: x.quadratic_representation(x.inverse())*Qx
1162 ....: == Qx*x.quadratic_representation(x.inverse()) )
1163 True
1164
1165 sage: not x.is_invertible() or (
1166 ....: x.quadratic_representation(x.inverse())*Qx
1167 ....: ==
1168 ....: 2*Lx*Qex - Qx )
1169 True
1170
1171 sage: 2*Lx*Qex - Qx == Lxx
1172 True
1173
1174 Property 5:
1175
1176 sage: Qy(x).quadratic_representation() == Qy*Qx*Qy
1177 True
1178
1179 Property 6:
1180
1181 sage: Qxn == (Qx)^n
1182 True
1183
1184 Property 7:
1185
1186 sage: not x.is_invertible() or (
1187 ....: Qx*x.inverse().operator() == Lx )
1188 True
1189
1190 Property 8:
1191
1192 sage: not x.operator_commutes_with(y) or (
1193 ....: Qx(y)^n == Qxn(y^n) )
1194 True
1195
1196 """
1197 if other is None:
1198 other=self
1199 elif not other in self.parent():
1200 raise TypeError("'other' must live in the same algebra")
1201
1202 L = self.operator()
1203 M = other.operator()
1204 return ( L*M + M*L - (self*other).operator() )
1205
1206
1207
1208 def spectral_decomposition(self):
1209 """
1210 Return the unique spectral decomposition of this element.
1211
1212 ALGORITHM:
1213
1214 Following Faraut and Korányi's Theorem III.1.1, we restrict this
1215 element's left-multiplication-by operator to the subalgebra it
1216 generates. We then compute the spectral decomposition of that
1217 operator, and the spectral projectors we get back must be the
1218 left-multiplication-by operators for the idempotents we
1219 seek. Thus applying them to the identity element gives us those
1220 idempotents.
1221
1222 Since the eigenvalues are required to be distinct, we take
1223 the spectral decomposition of the zero element to be zero
1224 times the identity element of the algebra (which is idempotent,
1225 obviously).
1226
1227 SETUP::
1228
1229 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1230
1231 EXAMPLES:
1232
1233 The spectral decomposition of the identity is ``1`` times itself,
1234 and the spectral decomposition of zero is ``0`` times the identity::
1235
1236 sage: J = RealSymmetricEJA(3)
1237 sage: J.one()
1238 e0 + e2 + e5
1239 sage: J.one().spectral_decomposition()
1240 [(1, e0 + e2 + e5)]
1241 sage: J.zero().spectral_decomposition()
1242 [(0, e0 + e2 + e5)]
1243
1244 TESTS::
1245
1246 sage: J = RealSymmetricEJA(4)
1247 sage: x = sum(J.gens())
1248 sage: sd = x.spectral_decomposition()
1249 sage: l0 = sd[0][0]
1250 sage: l1 = sd[1][0]
1251 sage: c0 = sd[0][1]
1252 sage: c1 = sd[1][1]
1253 sage: c0.inner_product(c1) == 0
1254 True
1255 sage: c0.is_idempotent()
1256 True
1257 sage: c1.is_idempotent()
1258 True
1259 sage: c0 + c1 == J.one()
1260 True
1261 sage: l0*c0 + l1*c1 == x
1262 True
1263
1264 The spectral decomposition should work in subalgebras, too::
1265
1266 sage: J = RealSymmetricEJA(4)
1267 sage: (e0, e1, e2, e3, e4, e5, e6, e7, e8, e9) = J.gens()
1268 sage: A = 2*e5 - 2*e8
1269 sage: (lambda1, c1) = A.spectral_decomposition()[1]
1270 sage: (J0, J5, J1) = J.peirce_decomposition(c1)
1271 sage: (f0, f1, f2) = J1.gens()
1272 sage: f0.spectral_decomposition()
1273 [(0, f2), (1, f0)]
1274
1275 """
1276 A = self.subalgebra_generated_by(orthonormalize_basis=True)
1277 result = []
1278 for (evalue, proj) in A(self).operator().spectral_decomposition():
1279 result.append( (evalue, proj(A.one()).superalgebra_element()) )
1280 return result
1281
1282 def subalgebra_generated_by(self, orthonormalize_basis=False):
1283 """
1284 Return the associative subalgebra of the parent EJA generated
1285 by this element.
1286
1287 Since our parent algebra is unital, we want "subalgebra" to mean
1288 "unital subalgebra" as well; thus the subalgebra that an element
1289 generates will itself be a Euclidean Jordan algebra after
1290 restricting the algebra operations appropriately. This is the
1291 subalgebra that Faraut and Korányi work with in section II.2, for
1292 example.
1293
1294 SETUP::
1295
1296 sage: from mjo.eja.eja_algebra import random_eja
1297
1298 TESTS:
1299
1300 This subalgebra, being composed of only powers, is associative::
1301
1302 sage: set_random_seed()
1303 sage: x0 = random_eja().random_element()
1304 sage: A = x0.subalgebra_generated_by()
1305 sage: x,y,z = A.random_elements(3)
1306 sage: (x*y)*z == x*(y*z)
1307 True
1308
1309 Squaring in the subalgebra should work the same as in
1310 the superalgebra::
1311
1312 sage: set_random_seed()
1313 sage: x = random_eja().random_element()
1314 sage: A = x.subalgebra_generated_by()
1315 sage: A(x^2) == A(x)*A(x)
1316 True
1317
1318 By definition, the subalgebra generated by the zero element is
1319 the one-dimensional algebra generated by the identity
1320 element... unless the original algebra was trivial, in which
1321 case the subalgebra is trivial too::
1322
1323 sage: set_random_seed()
1324 sage: A = random_eja().zero().subalgebra_generated_by()
1325 sage: (A.is_trivial() and A.dimension() == 0) or A.dimension() == 1
1326 True
1327
1328 """
1329 from mjo.eja.eja_element_subalgebra import FiniteDimensionalEuclideanJordanElementSubalgebra
1330 return FiniteDimensionalEuclideanJordanElementSubalgebra(self, orthonormalize_basis)
1331
1332
1333 def subalgebra_idempotent(self):
1334 """
1335 Find an idempotent in the associative subalgebra I generate
1336 using Proposition 2.3.5 in Baes.
1337
1338 SETUP::
1339
1340 sage: from mjo.eja.eja_algebra import random_eja
1341
1342 TESTS:
1343
1344 Ensure that we can find an idempotent in a non-trivial algebra
1345 where there are non-nilpotent elements, or that we get the dumb
1346 solution in the trivial algebra::
1347
1348 sage: set_random_seed()
1349 sage: J = random_eja()
1350 sage: x = J.random_element()
1351 sage: while x.is_nilpotent() and not J.is_trivial():
1352 ....: x = J.random_element()
1353 sage: c = x.subalgebra_idempotent()
1354 sage: c^2 == c
1355 True
1356
1357 """
1358 if self.parent().is_trivial():
1359 return self
1360
1361 if self.is_nilpotent():
1362 raise ValueError("this only works with non-nilpotent elements!")
1363
1364 J = self.subalgebra_generated_by()
1365 u = J(self)
1366
1367 # The image of the matrix of left-u^m-multiplication
1368 # will be minimal for some natural number s...
1369 s = 0
1370 minimal_dim = J.dimension()
1371 for i in range(1, minimal_dim):
1372 this_dim = (u**i).operator().matrix().image().dimension()
1373 if this_dim < minimal_dim:
1374 minimal_dim = this_dim
1375 s = i
1376
1377 # Now minimal_matrix should correspond to the smallest
1378 # non-zero subspace in Baes's (or really, Koecher's)
1379 # proposition.
1380 #
1381 # However, we need to restrict the matrix to work on the
1382 # subspace... or do we? Can't we just solve, knowing that
1383 # A(c) = u^(s+1) should have a solution in the big space,
1384 # too?
1385 #
1386 # Beware, solve_right() means that we're using COLUMN vectors.
1387 # Our FiniteDimensionalAlgebraElement superclass uses rows.
1388 u_next = u**(s+1)
1389 A = u_next.operator().matrix()
1390 c = J.from_vector(A.solve_right(u_next.to_vector()))
1391
1392 # Now c is the idempotent we want, but it still lives in the subalgebra.
1393 return c.superalgebra_element()
1394
1395
1396 def trace(self):
1397 """
1398 Return my trace, the sum of my eigenvalues.
1399
1400 In a trivial algebra, however you want to look at it, the trace is
1401 an empty sum for which we declare the result to be zero.
1402
1403 SETUP::
1404
1405 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1406 ....: HadamardEJA,
1407 ....: TrivialEJA,
1408 ....: random_eja)
1409
1410 EXAMPLES::
1411
1412 sage: J = TrivialEJA()
1413 sage: J.zero().trace()
1414 0
1415
1416 ::
1417 sage: J = JordanSpinEJA(3)
1418 sage: x = sum(J.gens())
1419 sage: x.trace()
1420 2
1421
1422 ::
1423
1424 sage: J = HadamardEJA(5)
1425 sage: J.one().trace()
1426 5
1427
1428 TESTS:
1429
1430 The trace of an element is a real number::
1431
1432 sage: set_random_seed()
1433 sage: J = random_eja()
1434 sage: J.random_element().trace() in RLF
1435 True
1436
1437 """
1438 P = self.parent()
1439 r = P.rank()
1440
1441 if r == 0:
1442 # Special case for the trivial algebra where
1443 # the trace is an empty sum.
1444 return P.base_ring().zero()
1445
1446 p = P._charpoly_coefficients()[r-1]
1447 # The _charpoly_coeff function already adds the factor of
1448 # -1 to ensure that _charpoly_coeff(r-1) is really what
1449 # appears in front of t^{r-1} in the charpoly. However,
1450 # we want the negative of THAT for the trace.
1451 return -p(*self.to_vector())
1452
1453
1454 def trace_inner_product(self, other):
1455 """
1456 Return the trace inner product of myself and ``other``.
1457
1458 SETUP::
1459
1460 sage: from mjo.eja.eja_algebra import random_eja
1461
1462 TESTS:
1463
1464 The trace inner product is commutative, bilinear, and associative::
1465
1466 sage: set_random_seed()
1467 sage: J = random_eja()
1468 sage: x,y,z = J.random_elements(3)
1469 sage: # commutative
1470 sage: x.trace_inner_product(y) == y.trace_inner_product(x)
1471 True
1472 sage: # bilinear
1473 sage: a = J.base_ring().random_element();
1474 sage: actual = (a*(x+z)).trace_inner_product(y)
1475 sage: expected = ( a*x.trace_inner_product(y) +
1476 ....: a*z.trace_inner_product(y) )
1477 sage: actual == expected
1478 True
1479 sage: actual = x.trace_inner_product(a*(y+z))
1480 sage: expected = ( a*x.trace_inner_product(y) +
1481 ....: a*x.trace_inner_product(z) )
1482 sage: actual == expected
1483 True
1484 sage: # associative
1485 sage: (x*y).trace_inner_product(z) == y.trace_inner_product(x*z)
1486 True
1487
1488 """
1489 if not other in self.parent():
1490 raise TypeError("'other' must live in the same algebra")
1491
1492 return (self*other).trace()
1493
1494
1495 def trace_norm(self):
1496 """
1497 The norm of this element with respect to :meth:`trace_inner_product`.
1498
1499 SETUP::
1500
1501 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1502 ....: HadamardEJA)
1503
1504 EXAMPLES::
1505
1506 sage: J = HadamardEJA(2)
1507 sage: x = sum(J.gens())
1508 sage: x.trace_norm()
1509 1.414213562373095?
1510
1511 ::
1512
1513 sage: J = JordanSpinEJA(4)
1514 sage: x = sum(J.gens())
1515 sage: x.trace_norm()
1516 2.828427124746190?
1517
1518 """
1519 return self.trace_inner_product(self).sqrt()