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