]> gitweb.michael.orlitzky.com - sage.d.git/blob - mjo/eja/eja_algebra.py
makefile: add --long to the test command.
[sage.d.git] / mjo / eja / eja_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 SETUP::
9
10 sage: from mjo.eja.eja_algebra import random_eja
11
12 EXAMPLES::
13
14 sage: random_eja()
15 Euclidean Jordan algebra of dimension...
16
17 """
18
19 from itertools import repeat
20
21 from sage.algebras.quatalg.quaternion_algebra import QuaternionAlgebra
22 from sage.categories.magmatic_algebras import MagmaticAlgebras
23 from sage.combinat.free_module import CombinatorialFreeModule
24 from sage.matrix.constructor import matrix
25 from sage.matrix.matrix_space import MatrixSpace
26 from sage.misc.cachefunc import cached_method
27 from sage.misc.table import table
28 from sage.modules.free_module import FreeModule, VectorSpace
29 from sage.rings.all import (ZZ, QQ, AA, QQbar, RR, RLF, CLF,
30 PolynomialRing,
31 QuadraticField)
32 from mjo.eja.eja_element import FiniteDimensionalEuclideanJordanAlgebraElement
33 from mjo.eja.eja_operator import FiniteDimensionalEuclideanJordanAlgebraOperator
34 from mjo.eja.eja_utils import _mat2vec
35
36 class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
37
38 def _coerce_map_from_base_ring(self):
39 """
40 Disable the map from the base ring into the algebra.
41
42 Performing a nonsense conversion like this automatically
43 is counterpedagogical. The fallback is to try the usual
44 element constructor, which should also fail.
45
46 SETUP::
47
48 sage: from mjo.eja.eja_algebra import random_eja
49
50 TESTS::
51
52 sage: set_random_seed()
53 sage: J = random_eja()
54 sage: J(1)
55 Traceback (most recent call last):
56 ...
57 ValueError: not an element of this algebra
58
59 """
60 return None
61
62 def __init__(self,
63 field,
64 mult_table,
65 prefix='e',
66 category=None,
67 matrix_basis=None,
68 check_field=True,
69 check_axioms=True):
70 """
71 SETUP::
72
73 sage: from mjo.eja.eja_algebra import (
74 ....: FiniteDimensionalEuclideanJordanAlgebra,
75 ....: JordanSpinEJA,
76 ....: random_eja)
77
78 EXAMPLES:
79
80 By definition, Jordan multiplication commutes::
81
82 sage: set_random_seed()
83 sage: J = random_eja()
84 sage: x,y = J.random_elements(2)
85 sage: x*y == y*x
86 True
87
88 TESTS:
89
90 The ``field`` we're given must be real with ``check_field=True``::
91
92 sage: JordanSpinEJA(2,QQbar)
93 Traceback (most recent call last):
94 ...
95 ValueError: scalar field is not real
96
97 The multiplication table must be square with ``check_axioms=True``::
98
99 sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,((),()))
100 Traceback (most recent call last):
101 ...
102 ValueError: multiplication table is not square
103
104 """
105 if check_field:
106 if not field.is_subring(RR):
107 # Note: this does return true for the real algebraic
108 # field, the rationals, and any quadratic field where
109 # we've specified a real embedding.
110 raise ValueError("scalar field is not real")
111
112 # The multiplication table had better be square
113 n = len(mult_table)
114 if check_axioms:
115 if not all( len(l) == n for l in mult_table ):
116 raise ValueError("multiplication table is not square")
117
118 self._matrix_basis = matrix_basis
119
120 if category is None:
121 category = MagmaticAlgebras(field).FiniteDimensional()
122 category = category.WithBasis().Unital()
123
124 fda = super(FiniteDimensionalEuclideanJordanAlgebra, self)
125 fda.__init__(field,
126 range(n),
127 prefix=prefix,
128 category=category)
129 self.print_options(bracket='')
130
131 # The multiplication table we're given is necessarily in terms
132 # of vectors, because we don't have an algebra yet for
133 # anything to be an element of. However, it's faster in the
134 # long run to have the multiplication table be in terms of
135 # algebra elements. We do this after calling the superclass
136 # constructor so that from_vector() knows what to do.
137 self._multiplication_table = [
138 list(map(lambda x: self.from_vector(x), ls))
139 for ls in mult_table
140 ]
141
142 if check_axioms:
143 if not self._is_commutative():
144 raise ValueError("algebra is not commutative")
145 if not self._is_jordanian():
146 raise ValueError("Jordan identity does not hold")
147 if not self._inner_product_is_associative():
148 raise ValueError("inner product is not associative")
149
150 def _element_constructor_(self, elt):
151 """
152 Construct an element of this algebra from its vector or matrix
153 representation.
154
155 This gets called only after the parent element _call_ method
156 fails to find a coercion for the argument.
157
158 SETUP::
159
160 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
161 ....: HadamardEJA,
162 ....: RealSymmetricEJA)
163
164 EXAMPLES:
165
166 The identity in `S^n` is converted to the identity in the EJA::
167
168 sage: J = RealSymmetricEJA(3)
169 sage: I = matrix.identity(QQ,3)
170 sage: J(I) == J.one()
171 True
172
173 This skew-symmetric matrix can't be represented in the EJA::
174
175 sage: J = RealSymmetricEJA(3)
176 sage: A = matrix(QQ,3, lambda i,j: i-j)
177 sage: J(A)
178 Traceback (most recent call last):
179 ...
180 ValueError: not an element of this algebra
181
182 TESTS:
183
184 Ensure that we can convert any element of the two non-matrix
185 simple algebras (whose matrix representations are columns)
186 back and forth faithfully::
187
188 sage: set_random_seed()
189 sage: J = HadamardEJA.random_instance()
190 sage: x = J.random_element()
191 sage: J(x.to_vector().column()) == x
192 True
193 sage: J = JordanSpinEJA.random_instance()
194 sage: x = J.random_element()
195 sage: J(x.to_vector().column()) == x
196 True
197 """
198 msg = "not an element of this algebra"
199 if elt == 0:
200 # The superclass implementation of random_element()
201 # needs to be able to coerce "0" into the algebra.
202 return self.zero()
203 elif elt in self.base_ring():
204 # Ensure that no base ring -> algebra coercion is performed
205 # by this method. There's some stupidity in sage that would
206 # otherwise propagate to this method; for example, sage thinks
207 # that the integer 3 belongs to the space of 2-by-2 matrices.
208 raise ValueError(msg)
209
210 if elt not in self.matrix_space():
211 raise ValueError(msg)
212
213 # Thanks for nothing! Matrix spaces aren't vector spaces in
214 # Sage, so we have to figure out its matrix-basis coordinates
215 # ourselves. We use the basis space's ring instead of the
216 # element's ring because the basis space might be an algebraic
217 # closure whereas the base ring of the 3-by-3 identity matrix
218 # could be QQ instead of QQbar.
219 V = VectorSpace(self.base_ring(), elt.nrows()*elt.ncols())
220 W = V.span_of_basis( _mat2vec(s) for s in self.matrix_basis() )
221
222 try:
223 coords = W.coordinate_vector(_mat2vec(elt))
224 except ArithmeticError: # vector is not in free module
225 raise ValueError(msg)
226
227 return self.from_vector(coords)
228
229 def _repr_(self):
230 """
231 Return a string representation of ``self``.
232
233 SETUP::
234
235 sage: from mjo.eja.eja_algebra import JordanSpinEJA
236
237 TESTS:
238
239 Ensure that it says what we think it says::
240
241 sage: JordanSpinEJA(2, field=AA)
242 Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
243 sage: JordanSpinEJA(3, field=RDF)
244 Euclidean Jordan algebra of dimension 3 over Real Double Field
245
246 """
247 fmt = "Euclidean Jordan algebra of dimension {} over {}"
248 return fmt.format(self.dimension(), self.base_ring())
249
250 def product_on_basis(self, i, j):
251 return self._multiplication_table[i][j]
252
253 def _is_commutative(self):
254 r"""
255 Whether or not this algebra's multiplication table is commutative.
256
257 This method should of course always return ``True``, unless
258 this algebra was constructed with ``check_axioms=False`` and
259 passed an invalid multiplication table.
260 """
261 return all( self.product_on_basis(i,j) == self.product_on_basis(i,j)
262 for i in range(self.dimension())
263 for j in range(self.dimension()) )
264
265 def _is_jordanian(self):
266 r"""
267 Whether or not this algebra's multiplication table respects the
268 Jordan identity `(x^{2})(xy) = x(x^{2}y)`.
269
270 We only check one arrangement of `x` and `y`, so for a
271 ``True`` result to be truly true, you should also check
272 :meth:`_is_commutative`. This method should of course always
273 return ``True``, unless this algebra was constructed with
274 ``check_axioms=False`` and passed an invalid multiplication table.
275 """
276 return all( (self.monomial(i)**2)*(self.monomial(i)*self.monomial(j))
277 ==
278 (self.monomial(i))*((self.monomial(i)**2)*self.monomial(j))
279 for i in range(self.dimension())
280 for j in range(self.dimension()) )
281
282 def _inner_product_is_associative(self):
283 r"""
284 Return whether or not this algebra's inner product `B` is
285 associative; that is, whether or not `B(xy,z) = B(x,yz)`.
286
287 This method should of course always return ``True``, unless
288 this algebra was constructed with ``check_axioms=False`` and
289 passed an invalid multiplication table.
290 """
291
292 # Used to check whether or not something is zero in an inexact
293 # ring. This number is sufficient to allow the construction of
294 # QuaternionHermitianEJA(2, RDF) with check_axioms=True.
295 epsilon = 1e-16
296
297 for i in range(self.dimension()):
298 for j in range(self.dimension()):
299 for k in range(self.dimension()):
300 x = self.monomial(i)
301 y = self.monomial(j)
302 z = self.monomial(k)
303 diff = (x*y).inner_product(z) - x.inner_product(y*z)
304
305 if self.base_ring().is_exact():
306 if diff != 0:
307 return False
308 else:
309 if diff.abs() > epsilon:
310 return False
311
312 return True
313
314 @cached_method
315 def characteristic_polynomial_of(self):
316 """
317 Return the algebra's "characteristic polynomial of" function,
318 which is itself a multivariate polynomial that, when evaluated
319 at the coordinates of some algebra element, returns that
320 element's characteristic polynomial.
321
322 The resulting polynomial has `n+1` variables, where `n` is the
323 dimension of this algebra. The first `n` variables correspond to
324 the coordinates of an algebra element: when evaluated at the
325 coordinates of an algebra element with respect to a certain
326 basis, the result is a univariate polynomial (in the one
327 remaining variable ``t``), namely the characteristic polynomial
328 of that element.
329
330 SETUP::
331
332 sage: from mjo.eja.eja_algebra import JordanSpinEJA, TrivialEJA
333
334 EXAMPLES:
335
336 The characteristic polynomial in the spin algebra is given in
337 Alizadeh, Example 11.11::
338
339 sage: J = JordanSpinEJA(3)
340 sage: p = J.characteristic_polynomial_of(); p
341 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
342 sage: xvec = J.one().to_vector()
343 sage: p(*xvec)
344 t^2 - 2*t + 1
345
346 By definition, the characteristic polynomial is a monic
347 degree-zero polynomial in a rank-zero algebra. Note that
348 Cayley-Hamilton is indeed satisfied since the polynomial
349 ``1`` evaluates to the identity element of the algebra on
350 any argument::
351
352 sage: J = TrivialEJA()
353 sage: J.characteristic_polynomial_of()
354 1
355
356 """
357 r = self.rank()
358 n = self.dimension()
359
360 # The list of coefficient polynomials a_0, a_1, a_2, ..., a_(r-1).
361 a = self._charpoly_coefficients()
362
363 # We go to a bit of trouble here to reorder the
364 # indeterminates, so that it's easier to evaluate the
365 # characteristic polynomial at x's coordinates and get back
366 # something in terms of t, which is what we want.
367 S = PolynomialRing(self.base_ring(),'t')
368 t = S.gen(0)
369 if r > 0:
370 R = a[0].parent()
371 S = PolynomialRing(S, R.variable_names())
372 t = S(t)
373
374 return (t**r + sum( a[k]*(t**k) for k in range(r) ))
375
376 def coordinate_polynomial_ring(self):
377 r"""
378 The multivariate polynomial ring in which this algebra's
379 :meth:`characteristic_polynomial_of` lives.
380
381 SETUP::
382
383 sage: from mjo.eja.eja_algebra import (HadamardEJA,
384 ....: RealSymmetricEJA)
385
386 EXAMPLES::
387
388 sage: J = HadamardEJA(2)
389 sage: J.coordinate_polynomial_ring()
390 Multivariate Polynomial Ring in X1, X2...
391 sage: J = RealSymmetricEJA(3,QQ)
392 sage: J.coordinate_polynomial_ring()
393 Multivariate Polynomial Ring in X1, X2, X3, X4, X5, X6...
394
395 """
396 var_names = tuple( "X%d" % z for z in range(1, self.dimension()+1) )
397 return PolynomialRing(self.base_ring(), var_names)
398
399 def inner_product(self, x, y):
400 """
401 The inner product associated with this Euclidean Jordan algebra.
402
403 Defaults to the trace inner product, but can be overridden by
404 subclasses if they are sure that the necessary properties are
405 satisfied.
406
407 SETUP::
408
409 sage: from mjo.eja.eja_algebra import (random_eja,
410 ....: HadamardEJA,
411 ....: BilinearFormEJA)
412
413 EXAMPLES:
414
415 Our inner product is "associative," which means the following for
416 a symmetric bilinear form::
417
418 sage: set_random_seed()
419 sage: J = random_eja()
420 sage: x,y,z = J.random_elements(3)
421 sage: (x*y).inner_product(z) == y.inner_product(x*z)
422 True
423
424 TESTS:
425
426 Ensure that this is the usual inner product for the algebras
427 over `R^n`::
428
429 sage: set_random_seed()
430 sage: J = HadamardEJA.random_instance()
431 sage: x,y = J.random_elements(2)
432 sage: actual = x.inner_product(y)
433 sage: expected = x.to_vector().inner_product(y.to_vector())
434 sage: actual == expected
435 True
436
437 Ensure that this is one-half of the trace inner-product in a
438 BilinearFormEJA that isn't just the reals (when ``n`` isn't
439 one). This is in Faraut and Koranyi, and also my "On the
440 symmetry..." paper::
441
442 sage: set_random_seed()
443 sage: J = BilinearFormEJA.random_instance()
444 sage: n = J.dimension()
445 sage: x = J.random_element()
446 sage: y = J.random_element()
447 sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
448 True
449 """
450 B = self._inner_product_matrix
451 return (B*x.to_vector()).inner_product(y.to_vector())
452
453
454 def is_trivial(self):
455 """
456 Return whether or not this algebra is trivial.
457
458 A trivial algebra contains only the zero element.
459
460 SETUP::
461
462 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
463 ....: TrivialEJA)
464
465 EXAMPLES::
466
467 sage: J = ComplexHermitianEJA(3)
468 sage: J.is_trivial()
469 False
470
471 ::
472
473 sage: J = TrivialEJA()
474 sage: J.is_trivial()
475 True
476
477 """
478 return self.dimension() == 0
479
480
481 def multiplication_table(self):
482 """
483 Return a visual representation of this algebra's multiplication
484 table (on basis elements).
485
486 SETUP::
487
488 sage: from mjo.eja.eja_algebra import JordanSpinEJA
489
490 EXAMPLES::
491
492 sage: J = JordanSpinEJA(4)
493 sage: J.multiplication_table()
494 +----++----+----+----+----+
495 | * || e0 | e1 | e2 | e3 |
496 +====++====+====+====+====+
497 | e0 || e0 | e1 | e2 | e3 |
498 +----++----+----+----+----+
499 | e1 || e1 | e0 | 0 | 0 |
500 +----++----+----+----+----+
501 | e2 || e2 | 0 | e0 | 0 |
502 +----++----+----+----+----+
503 | e3 || e3 | 0 | 0 | e0 |
504 +----++----+----+----+----+
505
506 """
507 M = list(self._multiplication_table) # copy
508 for i in range(len(M)):
509 # M had better be "square"
510 M[i] = [self.monomial(i)] + M[i]
511 M = [["*"] + list(self.gens())] + M
512 return table(M, header_row=True, header_column=True, frame=True)
513
514
515 def matrix_basis(self):
516 """
517 Return an (often more natural) representation of this algebras
518 basis as an ordered tuple of matrices.
519
520 Every finite-dimensional Euclidean Jordan Algebra is a, up to
521 Jordan isomorphism, a direct sum of five simple
522 algebras---four of which comprise Hermitian matrices. And the
523 last type of algebra can of course be thought of as `n`-by-`1`
524 column matrices (ambiguusly called column vectors) to avoid
525 special cases. As a result, matrices (and column vectors) are
526 a natural representation format for Euclidean Jordan algebra
527 elements.
528
529 But, when we construct an algebra from a basis of matrices,
530 those matrix representations are lost in favor of coordinate
531 vectors *with respect to* that basis. We could eventually
532 convert back if we tried hard enough, but having the original
533 representations handy is valuable enough that we simply store
534 them and return them from this method.
535
536 Why implement this for non-matrix algebras? Avoiding special
537 cases for the :class:`BilinearFormEJA` pays with simplicity in
538 its own right. But mainly, we would like to be able to assume
539 that elements of a :class:`DirectSumEJA` can be displayed
540 nicely, without having to have special classes for direct sums
541 one of whose components was a matrix algebra.
542
543 SETUP::
544
545 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
546 ....: RealSymmetricEJA)
547
548 EXAMPLES::
549
550 sage: J = RealSymmetricEJA(2)
551 sage: J.basis()
552 Finite family {0: e0, 1: e1, 2: e2}
553 sage: J.matrix_basis()
554 (
555 [1 0] [ 0 0.7071067811865475?] [0 0]
556 [0 0], [0.7071067811865475? 0], [0 1]
557 )
558
559 ::
560
561 sage: J = JordanSpinEJA(2)
562 sage: J.basis()
563 Finite family {0: e0, 1: e1}
564 sage: J.matrix_basis()
565 (
566 [1] [0]
567 [0], [1]
568 )
569 """
570 if self._matrix_basis is None:
571 M = self.matrix_space()
572 return tuple( M(b.to_vector()) for b in self.basis() )
573 else:
574 return self._matrix_basis
575
576
577 def matrix_space(self):
578 """
579 Return the matrix space in which this algebra's elements live, if
580 we think of them as matrices (including column vectors of the
581 appropriate size).
582
583 Generally this will be an `n`-by-`1` column-vector space,
584 except when the algebra is trivial. There it's `n`-by-`n`
585 (where `n` is zero), to ensure that two elements of the matrix
586 space (empty matrices) can be multiplied.
587
588 Matrix algebras override this with something more useful.
589 """
590 if self.is_trivial():
591 return MatrixSpace(self.base_ring(), 0)
592 elif self._matrix_basis is None or len(self._matrix_basis) == 0:
593 return MatrixSpace(self.base_ring(), self.dimension(), 1)
594 else:
595 return self._matrix_basis[0].matrix_space()
596
597
598 @cached_method
599 def one(self):
600 """
601 Return the unit element of this algebra.
602
603 SETUP::
604
605 sage: from mjo.eja.eja_algebra import (HadamardEJA,
606 ....: random_eja)
607
608 EXAMPLES::
609
610 sage: J = HadamardEJA(5)
611 sage: J.one()
612 e0 + e1 + e2 + e3 + e4
613
614 TESTS:
615
616 The identity element acts like the identity::
617
618 sage: set_random_seed()
619 sage: J = random_eja()
620 sage: x = J.random_element()
621 sage: J.one()*x == x and x*J.one() == x
622 True
623
624 The matrix of the unit element's operator is the identity::
625
626 sage: set_random_seed()
627 sage: J = random_eja()
628 sage: actual = J.one().operator().matrix()
629 sage: expected = matrix.identity(J.base_ring(), J.dimension())
630 sage: actual == expected
631 True
632
633 Ensure that the cached unit element (often precomputed by
634 hand) agrees with the computed one::
635
636 sage: set_random_seed()
637 sage: J = random_eja()
638 sage: cached = J.one()
639 sage: J.one.clear_cache()
640 sage: J.one() == cached
641 True
642
643 """
644 # We can brute-force compute the matrices of the operators
645 # that correspond to the basis elements of this algebra.
646 # If some linear combination of those basis elements is the
647 # algebra identity, then the same linear combination of
648 # their matrices has to be the identity matrix.
649 #
650 # Of course, matrices aren't vectors in sage, so we have to
651 # appeal to the "long vectors" isometry.
652 oper_vecs = [ _mat2vec(g.operator().matrix()) for g in self.gens() ]
653
654 # Now we use basic linear algebra to find the coefficients,
655 # of the matrices-as-vectors-linear-combination, which should
656 # work for the original algebra basis too.
657 A = matrix(self.base_ring(), oper_vecs)
658
659 # We used the isometry on the left-hand side already, but we
660 # still need to do it for the right-hand side. Recall that we
661 # wanted something that summed to the identity matrix.
662 b = _mat2vec( matrix.identity(self.base_ring(), self.dimension()) )
663
664 # Now if there's an identity element in the algebra, this
665 # should work. We solve on the left to avoid having to
666 # transpose the matrix "A".
667 return self.from_vector(A.solve_left(b))
668
669
670 def peirce_decomposition(self, c):
671 """
672 The Peirce decomposition of this algebra relative to the
673 idempotent ``c``.
674
675 In the future, this can be extended to a complete system of
676 orthogonal idempotents.
677
678 INPUT:
679
680 - ``c`` -- an idempotent of this algebra.
681
682 OUTPUT:
683
684 A triple (J0, J5, J1) containing two subalgebras and one subspace
685 of this algebra,
686
687 - ``J0`` -- the algebra on the eigenspace of ``c.operator()``
688 corresponding to the eigenvalue zero.
689
690 - ``J5`` -- the eigenspace (NOT a subalgebra) of ``c.operator()``
691 corresponding to the eigenvalue one-half.
692
693 - ``J1`` -- the algebra on the eigenspace of ``c.operator()``
694 corresponding to the eigenvalue one.
695
696 These are the only possible eigenspaces for that operator, and this
697 algebra is a direct sum of them. The spaces ``J0`` and ``J1`` are
698 orthogonal, and are subalgebras of this algebra with the appropriate
699 restrictions.
700
701 SETUP::
702
703 sage: from mjo.eja.eja_algebra import random_eja, RealSymmetricEJA
704
705 EXAMPLES:
706
707 The canonical example comes from the symmetric matrices, which
708 decompose into diagonal and off-diagonal parts::
709
710 sage: J = RealSymmetricEJA(3)
711 sage: C = matrix(QQ, [ [1,0,0],
712 ....: [0,1,0],
713 ....: [0,0,0] ])
714 sage: c = J(C)
715 sage: J0,J5,J1 = J.peirce_decomposition(c)
716 sage: J0
717 Euclidean Jordan algebra of dimension 1...
718 sage: J5
719 Vector space of degree 6 and dimension 2...
720 sage: J1
721 Euclidean Jordan algebra of dimension 3...
722 sage: J0.one().to_matrix()
723 [0 0 0]
724 [0 0 0]
725 [0 0 1]
726 sage: orig_df = AA.options.display_format
727 sage: AA.options.display_format = 'radical'
728 sage: J.from_vector(J5.basis()[0]).to_matrix()
729 [ 0 0 1/2*sqrt(2)]
730 [ 0 0 0]
731 [1/2*sqrt(2) 0 0]
732 sage: J.from_vector(J5.basis()[1]).to_matrix()
733 [ 0 0 0]
734 [ 0 0 1/2*sqrt(2)]
735 [ 0 1/2*sqrt(2) 0]
736 sage: AA.options.display_format = orig_df
737 sage: J1.one().to_matrix()
738 [1 0 0]
739 [0 1 0]
740 [0 0 0]
741
742 TESTS:
743
744 Every algebra decomposes trivially with respect to its identity
745 element::
746
747 sage: set_random_seed()
748 sage: J = random_eja()
749 sage: J0,J5,J1 = J.peirce_decomposition(J.one())
750 sage: J0.dimension() == 0 and J5.dimension() == 0
751 True
752 sage: J1.superalgebra() == J and J1.dimension() == J.dimension()
753 True
754
755 The decomposition is into eigenspaces, and its components are
756 therefore necessarily orthogonal. Moreover, the identity
757 elements in the two subalgebras are the projections onto their
758 respective subspaces of the superalgebra's identity element::
759
760 sage: set_random_seed()
761 sage: J = random_eja()
762 sage: x = J.random_element()
763 sage: if not J.is_trivial():
764 ....: while x.is_nilpotent():
765 ....: x = J.random_element()
766 sage: c = x.subalgebra_idempotent()
767 sage: J0,J5,J1 = J.peirce_decomposition(c)
768 sage: ipsum = 0
769 sage: for (w,y,z) in zip(J0.basis(), J5.basis(), J1.basis()):
770 ....: w = w.superalgebra_element()
771 ....: y = J.from_vector(y)
772 ....: z = z.superalgebra_element()
773 ....: ipsum += w.inner_product(y).abs()
774 ....: ipsum += w.inner_product(z).abs()
775 ....: ipsum += y.inner_product(z).abs()
776 sage: ipsum
777 0
778 sage: J1(c) == J1.one()
779 True
780 sage: J0(J.one() - c) == J0.one()
781 True
782
783 """
784 if not c.is_idempotent():
785 raise ValueError("element is not idempotent: %s" % c)
786
787 from mjo.eja.eja_subalgebra import FiniteDimensionalEuclideanJordanSubalgebra
788
789 # Default these to what they should be if they turn out to be
790 # trivial, because eigenspaces_left() won't return eigenvalues
791 # corresponding to trivial spaces (e.g. it returns only the
792 # eigenspace corresponding to lambda=1 if you take the
793 # decomposition relative to the identity element).
794 trivial = FiniteDimensionalEuclideanJordanSubalgebra(self, ())
795 J0 = trivial # eigenvalue zero
796 J5 = VectorSpace(self.base_ring(), 0) # eigenvalue one-half
797 J1 = trivial # eigenvalue one
798
799 for (eigval, eigspace) in c.operator().matrix().right_eigenspaces():
800 if eigval == ~(self.base_ring()(2)):
801 J5 = eigspace
802 else:
803 gens = tuple( self.from_vector(b) for b in eigspace.basis() )
804 subalg = FiniteDimensionalEuclideanJordanSubalgebra(self,
805 gens,
806 check_axioms=False)
807 if eigval == 0:
808 J0 = subalg
809 elif eigval == 1:
810 J1 = subalg
811 else:
812 raise ValueError("unexpected eigenvalue: %s" % eigval)
813
814 return (J0, J5, J1)
815
816
817 def random_element(self, thorough=False):
818 r"""
819 Return a random element of this algebra.
820
821 Our algebra superclass method only returns a linear
822 combination of at most two basis elements. We instead
823 want the vector space "random element" method that
824 returns a more diverse selection.
825
826 INPUT:
827
828 - ``thorough`` -- (boolean; default False) whether or not we
829 should generate irrational coefficients for the random
830 element when our base ring is irrational; this slows the
831 algebra operations to a crawl, but any truly random method
832 should include them
833
834 """
835 # For a general base ring... maybe we can trust this to do the
836 # right thing? Unlikely, but.
837 V = self.vector_space()
838 v = V.random_element()
839
840 if self.base_ring() is AA:
841 # The "random element" method of the algebraic reals is
842 # stupid at the moment, and only returns integers between
843 # -2 and 2, inclusive:
844 #
845 # https://trac.sagemath.org/ticket/30875
846 #
847 # Instead, we implement our own "random vector" method,
848 # and then coerce that into the algebra. We use the vector
849 # space degree here instead of the dimension because a
850 # subalgebra could (for example) be spanned by only two
851 # vectors, each with five coordinates. We need to
852 # generate all five coordinates.
853 if thorough:
854 v *= QQbar.random_element().real()
855 else:
856 v *= QQ.random_element()
857
858 return self.from_vector(V.coordinate_vector(v))
859
860 def random_elements(self, count, thorough=False):
861 """
862 Return ``count`` random elements as a tuple.
863
864 INPUT:
865
866 - ``thorough`` -- (boolean; default False) whether or not we
867 should generate irrational coefficients for the random
868 elements when our base ring is irrational; this slows the
869 algebra operations to a crawl, but any truly random method
870 should include them
871
872 SETUP::
873
874 sage: from mjo.eja.eja_algebra import JordanSpinEJA
875
876 EXAMPLES::
877
878 sage: J = JordanSpinEJA(3)
879 sage: x,y,z = J.random_elements(3)
880 sage: all( [ x in J, y in J, z in J ])
881 True
882 sage: len( J.random_elements(10) ) == 10
883 True
884
885 """
886 return tuple( self.random_element(thorough)
887 for idx in range(count) )
888
889
890 @cached_method
891 def _charpoly_coefficients(self):
892 r"""
893 The `r` polynomial coefficients of the "characteristic polynomial
894 of" function.
895 """
896 n = self.dimension()
897 R = self.coordinate_polynomial_ring()
898 vars = R.gens()
899 F = R.fraction_field()
900
901 def L_x_i_j(i,j):
902 # From a result in my book, these are the entries of the
903 # basis representation of L_x.
904 return sum( vars[k]*self.monomial(k).operator().matrix()[i,j]
905 for k in range(n) )
906
907 L_x = matrix(F, n, n, L_x_i_j)
908
909 r = None
910 if self.rank.is_in_cache():
911 r = self.rank()
912 # There's no need to pad the system with redundant
913 # columns if we *know* they'll be redundant.
914 n = r
915
916 # Compute an extra power in case the rank is equal to
917 # the dimension (otherwise, we would stop at x^(r-1)).
918 x_powers = [ (L_x**k)*self.one().to_vector()
919 for k in range(n+1) ]
920 A = matrix.column(F, x_powers[:n])
921 AE = A.extended_echelon_form()
922 E = AE[:,n:]
923 A_rref = AE[:,:n]
924 if r is None:
925 r = A_rref.rank()
926 b = x_powers[r]
927
928 # The theory says that only the first "r" coefficients are
929 # nonzero, and they actually live in the original polynomial
930 # ring and not the fraction field. We negate them because
931 # in the actual characteristic polynomial, they get moved
932 # to the other side where x^r lives.
933 return -A_rref.solve_right(E*b).change_ring(R)[:r]
934
935 @cached_method
936 def rank(self):
937 r"""
938 Return the rank of this EJA.
939
940 This is a cached method because we know the rank a priori for
941 all of the algebras we can construct. Thus we can avoid the
942 expensive ``_charpoly_coefficients()`` call unless we truly
943 need to compute the whole characteristic polynomial.
944
945 SETUP::
946
947 sage: from mjo.eja.eja_algebra import (HadamardEJA,
948 ....: JordanSpinEJA,
949 ....: RealSymmetricEJA,
950 ....: ComplexHermitianEJA,
951 ....: QuaternionHermitianEJA,
952 ....: random_eja)
953
954 EXAMPLES:
955
956 The rank of the Jordan spin algebra is always two::
957
958 sage: JordanSpinEJA(2).rank()
959 2
960 sage: JordanSpinEJA(3).rank()
961 2
962 sage: JordanSpinEJA(4).rank()
963 2
964
965 The rank of the `n`-by-`n` Hermitian real, complex, or
966 quaternion matrices is `n`::
967
968 sage: RealSymmetricEJA(4).rank()
969 4
970 sage: ComplexHermitianEJA(3).rank()
971 3
972 sage: QuaternionHermitianEJA(2).rank()
973 2
974
975 TESTS:
976
977 Ensure that every EJA that we know how to construct has a
978 positive integer rank, unless the algebra is trivial in
979 which case its rank will be zero::
980
981 sage: set_random_seed()
982 sage: J = random_eja()
983 sage: r = J.rank()
984 sage: r in ZZ
985 True
986 sage: r > 0 or (r == 0 and J.is_trivial())
987 True
988
989 Ensure that computing the rank actually works, since the ranks
990 of all simple algebras are known and will be cached by default::
991
992 sage: set_random_seed() # long time
993 sage: J = random_eja() # long time
994 sage: caches = J.rank() # long time
995 sage: J.rank.clear_cache() # long time
996 sage: J.rank() == cached # long time
997 True
998
999 """
1000 return len(self._charpoly_coefficients())
1001
1002
1003 def vector_space(self):
1004 """
1005 Return the vector space that underlies this algebra.
1006
1007 SETUP::
1008
1009 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1010
1011 EXAMPLES::
1012
1013 sage: J = RealSymmetricEJA(2)
1014 sage: J.vector_space()
1015 Vector space of dimension 3 over...
1016
1017 """
1018 return self.zero().to_vector().parent().ambient_vector_space()
1019
1020
1021 Element = FiniteDimensionalEuclideanJordanAlgebraElement
1022
1023
1024 class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
1025 r"""
1026 Algebras whose basis consists of vectors with rational
1027 entries. Equivalently, algebras whose multiplication tables
1028 contain only rational coefficients.
1029
1030 When an EJA has a basis that can be made rational, we can speed up
1031 the computation of its characteristic polynomial by doing it over
1032 ``QQ``. All of the named EJA constructors that we provide fall
1033 into this category.
1034 """
1035 @cached_method
1036 def _charpoly_coefficients(self):
1037 r"""
1038 Override the parent method with something that tries to compute
1039 over a faster (non-extension) field.
1040
1041 SETUP::
1042
1043 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1044
1045 EXAMPLES:
1046
1047 The base ring of the resulting polynomial coefficients is what
1048 it should be, and not the rationals (unless the algebra was
1049 already over the rationals)::
1050
1051 sage: J = JordanSpinEJA(3)
1052 sage: J._charpoly_coefficients()
1053 (X1^2 - X2^2 - X3^2, -2*X1)
1054 sage: a0 = J._charpoly_coefficients()[0]
1055 sage: J.base_ring()
1056 Algebraic Real Field
1057 sage: a0.base_ring()
1058 Algebraic Real Field
1059
1060 """
1061 if self.base_ring() is QQ:
1062 # There's no need to construct *another* algebra over the
1063 # rationals if this one is already over the rationals.
1064 superclass = super(RationalBasisEuclideanJordanAlgebra, self)
1065 return superclass._charpoly_coefficients()
1066
1067 mult_table = tuple(
1068 map(lambda x: x.to_vector(), ls)
1069 for ls in self._multiplication_table
1070 )
1071
1072 # Do the computation over the rationals. The answer will be
1073 # the same, because our basis coordinates are (essentially)
1074 # rational.
1075 J = FiniteDimensionalEuclideanJordanAlgebra(QQ,
1076 mult_table,
1077 check_field=False,
1078 check_axioms=False)
1079 a = J._charpoly_coefficients()
1080 return tuple(map(lambda x: x.change_ring(self.base_ring()), a))
1081
1082
1083 class ConcreteEuclideanJordanAlgebra:
1084 r"""
1085 A class for the Euclidean Jordan algebras that we know by name.
1086
1087 These are the Jordan algebras whose basis, multiplication table,
1088 rank, and so on are known a priori. More to the point, they are
1089 the Euclidean Jordan algebras for which we are able to conjure up
1090 a "random instance."
1091
1092 SETUP::
1093
1094 sage: from mjo.eja.eja_algebra import ConcreteEuclideanJordanAlgebra
1095
1096 TESTS:
1097
1098 Our basis is normalized with respect to the algebra's inner
1099 product, unless we specify otherwise::
1100
1101 sage: set_random_seed()
1102 sage: J = ConcreteEuclideanJordanAlgebra.random_instance()
1103 sage: all( b.norm() == 1 for b in J.gens() )
1104 True
1105
1106 Since our basis is orthonormal with respect to the algebra's inner
1107 product, and since we know that this algebra is an EJA, any
1108 left-multiplication operator's matrix will be symmetric because
1109 natural->EJA basis representation is an isometry and within the
1110 EJA the operator is self-adjoint by the Jordan axiom::
1111
1112 sage: set_random_seed()
1113 sage: J = ConcreteEuclideanJordanAlgebra.random_instance()
1114 sage: x = J.random_element()
1115 sage: x.operator().is_self_adjoint()
1116 True
1117 """
1118
1119 @staticmethod
1120 def _max_random_instance_size():
1121 """
1122 Return an integer "size" that is an upper bound on the size of
1123 this algebra when it is used in a random test
1124 case. Unfortunately, the term "size" is ambiguous -- when
1125 dealing with `R^n` under either the Hadamard or Jordan spin
1126 product, the "size" refers to the dimension `n`. When dealing
1127 with a matrix algebra (real symmetric or complex/quaternion
1128 Hermitian), it refers to the size of the matrix, which is far
1129 less than the dimension of the underlying vector space.
1130
1131 This method must be implemented in each subclass.
1132 """
1133 raise NotImplementedError
1134
1135 @classmethod
1136 def random_instance(cls, field=AA, **kwargs):
1137 """
1138 Return a random instance of this type of algebra.
1139
1140 This method should be implemented in each subclass.
1141 """
1142 from sage.misc.prandom import choice
1143 eja_class = choice(cls.__subclasses__())
1144 return eja_class.random_instance(field)
1145
1146
1147 class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
1148
1149 def __init__(self, field, basis, normalize_basis=True, **kwargs):
1150 """
1151 Compared to the superclass constructor, we take a basis instead of
1152 a multiplication table because the latter can be computed in terms
1153 of the former when the product is known (like it is here).
1154 """
1155 # Used in this class's fast _charpoly_coefficients() override.
1156 self._basis_normalizers = None
1157
1158 # We're going to loop through this a few times, so now's a good
1159 # time to ensure that it isn't a generator expression.
1160 basis = tuple(basis)
1161
1162 algebra_dim = len(basis)
1163 degree = 0 # size of the matrices
1164 if algebra_dim > 0:
1165 degree = basis[0].nrows()
1166
1167 if algebra_dim > 1 and normalize_basis:
1168 # We'll need sqrt(2) to normalize the basis, and this
1169 # winds up in the multiplication table, so the whole
1170 # algebra needs to be over the field extension.
1171 R = PolynomialRing(field, 'z')
1172 z = R.gen()
1173 p = z**2 - 2
1174 if p.is_irreducible():
1175 field = field.extension(p, 'sqrt2', embedding=RLF(2).sqrt())
1176 basis = tuple( s.change_ring(field) for s in basis )
1177 self._basis_normalizers = tuple(
1178 ~(self.matrix_inner_product(s,s).sqrt()) for s in basis )
1179 basis = tuple(s*c for (s,c) in zip(basis,self._basis_normalizers))
1180
1181 # Now compute the multiplication and inner product tables.
1182 # We have to do this *after* normalizing the basis, because
1183 # scaling affects the answers.
1184 V = VectorSpace(field, degree**2)
1185 W = V.span_of_basis( _mat2vec(s) for s in basis )
1186 mult_table = [[W.zero() for j in range(algebra_dim)]
1187 for i in range(algebra_dim)]
1188 ip_table = [[W.zero() for j in range(algebra_dim)]
1189 for i in range(algebra_dim)]
1190 for i in range(algebra_dim):
1191 for j in range(algebra_dim):
1192 mat_entry = (basis[i]*basis[j] + basis[j]*basis[i])/2
1193 mult_table[i][j] = W.coordinate_vector(_mat2vec(mat_entry))
1194
1195 try:
1196 # HACK: ignore the error here if we don't need the
1197 # inner product (as is the case when we construct
1198 # a dummy QQ-algebra for fast charpoly coefficients.
1199 ip_table[i][j] = self.matrix_inner_product(basis[i],
1200 basis[j])
1201 except:
1202 pass
1203
1204 try:
1205 # HACK PART DEUX
1206 self._inner_product_matrix = matrix(field,ip_table)
1207 except:
1208 pass
1209
1210 super(MatrixEuclideanJordanAlgebra, self).__init__(field,
1211 mult_table,
1212 matrix_basis=basis,
1213 **kwargs)
1214
1215 if algebra_dim == 0:
1216 self.one.set_cache(self.zero())
1217 else:
1218 n = basis[0].nrows()
1219 # The identity wrt (A,B) -> (AB + BA)/2 is independent of the
1220 # details of this algebra.
1221 self.one.set_cache(self(matrix.identity(field,n)))
1222
1223
1224 @cached_method
1225 def _charpoly_coefficients(self):
1226 r"""
1227 Override the parent method with something that tries to compute
1228 over a faster (non-extension) field.
1229 """
1230 if self._basis_normalizers is None or self.base_ring() is QQ:
1231 # We didn't normalize, or the basis we started with had
1232 # entries in a nice field already. Just compute the thing.
1233 return super(MatrixEuclideanJordanAlgebra, self)._charpoly_coefficients()
1234
1235 basis = ( (b/n) for (b,n) in zip(self.matrix_basis(),
1236 self._basis_normalizers) )
1237
1238 # Do this over the rationals and convert back at the end.
1239 # Only works because we know the entries of the basis are
1240 # integers. The argument ``check_axioms=False`` is required
1241 # because the trace inner-product method for this
1242 # class is a stub and can't actually be checked.
1243 J = MatrixEuclideanJordanAlgebra(QQ,
1244 basis,
1245 normalize_basis=False,
1246 check_field=False,
1247 check_axioms=False)
1248 a = J._charpoly_coefficients()
1249
1250 # Unfortunately, changing the basis does change the
1251 # coefficients of the characteristic polynomial, but since
1252 # these are really the coefficients of the "characteristic
1253 # polynomial of" function, everything is still nice and
1254 # unevaluated. It's therefore "obvious" how scaling the
1255 # basis affects the coordinate variables X1, X2, et
1256 # cetera. Scaling the first basis vector up by "n" adds a
1257 # factor of 1/n into every "X1" term, for example. So here
1258 # we simply undo the basis_normalizer scaling that we
1259 # performed earlier.
1260 #
1261 # The a[0] access here is safe because trivial algebras
1262 # won't have any basis normalizers and therefore won't
1263 # make it to this "else" branch.
1264 XS = a[0].parent().gens()
1265 subs_dict = { XS[i]: self._basis_normalizers[i]*XS[i]
1266 for i in range(len(XS)) }
1267 return tuple( a_i.subs(subs_dict) for a_i in a )
1268
1269
1270 @staticmethod
1271 def real_embed(M):
1272 """
1273 Embed the matrix ``M`` into a space of real matrices.
1274
1275 The matrix ``M`` can have entries in any field at the moment:
1276 the real numbers, complex numbers, or quaternions. And although
1277 they are not a field, we can probably support octonions at some
1278 point, too. This function returns a real matrix that "acts like"
1279 the original with respect to matrix multiplication; i.e.
1280
1281 real_embed(M*N) = real_embed(M)*real_embed(N)
1282
1283 """
1284 raise NotImplementedError
1285
1286
1287 @staticmethod
1288 def real_unembed(M):
1289 """
1290 The inverse of :meth:`real_embed`.
1291 """
1292 raise NotImplementedError
1293
1294 @classmethod
1295 def matrix_inner_product(cls,X,Y):
1296 Xu = cls.real_unembed(X)
1297 Yu = cls.real_unembed(Y)
1298 tr = (Xu*Yu).trace()
1299
1300 try:
1301 # Works in QQ, AA, RDF, et cetera.
1302 return tr.real()
1303 except AttributeError:
1304 # A quaternion doesn't have a real() method, but does
1305 # have coefficient_tuple() method that returns the
1306 # coefficients of 1, i, j, and k -- in that order.
1307 return tr.coefficient_tuple()[0]
1308
1309
1310 class RealMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
1311 @staticmethod
1312 def real_embed(M):
1313 """
1314 The identity function, for embedding real matrices into real
1315 matrices.
1316 """
1317 return M
1318
1319 @staticmethod
1320 def real_unembed(M):
1321 """
1322 The identity function, for unembedding real matrices from real
1323 matrices.
1324 """
1325 return M
1326
1327
1328 class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra,
1329 ConcreteEuclideanJordanAlgebra):
1330 """
1331 The rank-n simple EJA consisting of real symmetric n-by-n
1332 matrices, the usual symmetric Jordan product, and the trace inner
1333 product. It has dimension `(n^2 + n)/2` over the reals.
1334
1335 SETUP::
1336
1337 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1338
1339 EXAMPLES::
1340
1341 sage: J = RealSymmetricEJA(2)
1342 sage: e0, e1, e2 = J.gens()
1343 sage: e0*e0
1344 e0
1345 sage: e1*e1
1346 1/2*e0 + 1/2*e2
1347 sage: e2*e2
1348 e2
1349
1350 In theory, our "field" can be any subfield of the reals::
1351
1352 sage: RealSymmetricEJA(2, RDF)
1353 Euclidean Jordan algebra of dimension 3 over Real Double Field
1354 sage: RealSymmetricEJA(2, RR)
1355 Euclidean Jordan algebra of dimension 3 over Real Field with
1356 53 bits of precision
1357
1358 TESTS:
1359
1360 The dimension of this algebra is `(n^2 + n) / 2`::
1361
1362 sage: set_random_seed()
1363 sage: n_max = RealSymmetricEJA._max_random_instance_size()
1364 sage: n = ZZ.random_element(1, n_max)
1365 sage: J = RealSymmetricEJA(n)
1366 sage: J.dimension() == (n^2 + n)/2
1367 True
1368
1369 The Jordan multiplication is what we think it is::
1370
1371 sage: set_random_seed()
1372 sage: J = RealSymmetricEJA.random_instance()
1373 sage: x,y = J.random_elements(2)
1374 sage: actual = (x*y).to_matrix()
1375 sage: X = x.to_matrix()
1376 sage: Y = y.to_matrix()
1377 sage: expected = (X*Y + Y*X)/2
1378 sage: actual == expected
1379 True
1380 sage: J(expected) == x*y
1381 True
1382
1383 We can change the generator prefix::
1384
1385 sage: RealSymmetricEJA(3, prefix='q').gens()
1386 (q0, q1, q2, q3, q4, q5)
1387
1388 We can construct the (trivial) algebra of rank zero::
1389
1390 sage: RealSymmetricEJA(0)
1391 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1392
1393 """
1394 @classmethod
1395 def _denormalized_basis(cls, n, field):
1396 """
1397 Return a basis for the space of real symmetric n-by-n matrices.
1398
1399 SETUP::
1400
1401 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1402
1403 TESTS::
1404
1405 sage: set_random_seed()
1406 sage: n = ZZ.random_element(1,5)
1407 sage: B = RealSymmetricEJA._denormalized_basis(n,QQ)
1408 sage: all( M.is_symmetric() for M in B)
1409 True
1410
1411 """
1412 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1413 # coordinates.
1414 S = []
1415 for i in range(n):
1416 for j in range(i+1):
1417 Eij = matrix(field, n, lambda k,l: k==i and l==j)
1418 if i == j:
1419 Sij = Eij
1420 else:
1421 Sij = Eij + Eij.transpose()
1422 S.append(Sij)
1423 return S
1424
1425
1426 @staticmethod
1427 def _max_random_instance_size():
1428 return 4 # Dimension 10
1429
1430 @classmethod
1431 def random_instance(cls, field=AA, **kwargs):
1432 """
1433 Return a random instance of this type of algebra.
1434 """
1435 n = ZZ.random_element(cls._max_random_instance_size() + 1)
1436 return cls(n, field, **kwargs)
1437
1438 def __init__(self, n, field=AA, **kwargs):
1439 basis = self._denormalized_basis(n, field)
1440 super(RealSymmetricEJA, self).__init__(field,
1441 basis,
1442 check_axioms=False,
1443 **kwargs)
1444 self.rank.set_cache(n)
1445
1446
1447 class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
1448 @staticmethod
1449 def real_embed(M):
1450 """
1451 Embed the n-by-n complex matrix ``M`` into the space of real
1452 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1453 bi` to the block matrix ``[[a,b],[-b,a]]``.
1454
1455 SETUP::
1456
1457 sage: from mjo.eja.eja_algebra import \
1458 ....: ComplexMatrixEuclideanJordanAlgebra
1459
1460 EXAMPLES::
1461
1462 sage: F = QuadraticField(-1, 'I')
1463 sage: x1 = F(4 - 2*i)
1464 sage: x2 = F(1 + 2*i)
1465 sage: x3 = F(-i)
1466 sage: x4 = F(6)
1467 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1468 sage: ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1469 [ 4 -2| 1 2]
1470 [ 2 4|-2 1]
1471 [-----+-----]
1472 [ 0 -1| 6 0]
1473 [ 1 0| 0 6]
1474
1475 TESTS:
1476
1477 Embedding is a homomorphism (isomorphism, in fact)::
1478
1479 sage: set_random_seed()
1480 sage: n = ZZ.random_element(3)
1481 sage: F = QuadraticField(-1, 'I')
1482 sage: X = random_matrix(F, n)
1483 sage: Y = random_matrix(F, n)
1484 sage: Xe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X)
1485 sage: Ye = ComplexMatrixEuclideanJordanAlgebra.real_embed(Y)
1486 sage: XYe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1487 sage: Xe*Ye == XYe
1488 True
1489
1490 """
1491 n = M.nrows()
1492 if M.ncols() != n:
1493 raise ValueError("the matrix 'M' must be square")
1494
1495 # We don't need any adjoined elements...
1496 field = M.base_ring().base_ring()
1497
1498 blocks = []
1499 for z in M.list():
1500 a = z.list()[0] # real part, I guess
1501 b = z.list()[1] # imag part, I guess
1502 blocks.append(matrix(field, 2, [[a,b],[-b,a]]))
1503
1504 return matrix.block(field, n, blocks)
1505
1506
1507 @staticmethod
1508 def real_unembed(M):
1509 """
1510 The inverse of _embed_complex_matrix().
1511
1512 SETUP::
1513
1514 sage: from mjo.eja.eja_algebra import \
1515 ....: ComplexMatrixEuclideanJordanAlgebra
1516
1517 EXAMPLES::
1518
1519 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1520 ....: [-2, 1, -4, 3],
1521 ....: [ 9, 10, 11, 12],
1522 ....: [-10, 9, -12, 11] ])
1523 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(A)
1524 [ 2*I + 1 4*I + 3]
1525 [ 10*I + 9 12*I + 11]
1526
1527 TESTS:
1528
1529 Unembedding is the inverse of embedding::
1530
1531 sage: set_random_seed()
1532 sage: F = QuadraticField(-1, 'I')
1533 sage: M = random_matrix(F, 3)
1534 sage: Me = ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1535 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1536 True
1537
1538 """
1539 n = ZZ(M.nrows())
1540 if M.ncols() != n:
1541 raise ValueError("the matrix 'M' must be square")
1542 if not n.mod(2).is_zero():
1543 raise ValueError("the matrix 'M' must be a complex embedding")
1544
1545 # If "M" was normalized, its base ring might have roots
1546 # adjoined and they can stick around after unembedding.
1547 field = M.base_ring()
1548 R = PolynomialRing(field, 'z')
1549 z = R.gen()
1550 if field is AA:
1551 # Sage doesn't know how to embed AA into QQbar, i.e. how
1552 # to adjoin sqrt(-1) to AA.
1553 F = QQbar
1554 else:
1555 F = field.extension(z**2 + 1, 'I', embedding=CLF(-1).sqrt())
1556 i = F.gen()
1557
1558 # Go top-left to bottom-right (reading order), converting every
1559 # 2-by-2 block we see to a single complex element.
1560 elements = []
1561 for k in range(n/2):
1562 for j in range(n/2):
1563 submat = M[2*k:2*k+2,2*j:2*j+2]
1564 if submat[0,0] != submat[1,1]:
1565 raise ValueError('bad on-diagonal submatrix')
1566 if submat[0,1] != -submat[1,0]:
1567 raise ValueError('bad off-diagonal submatrix')
1568 z = submat[0,0] + submat[0,1]*i
1569 elements.append(z)
1570
1571 return matrix(F, n/2, elements)
1572
1573
1574 @classmethod
1575 def matrix_inner_product(cls,X,Y):
1576 """
1577 Compute a matrix inner product in this algebra directly from
1578 its real embedding.
1579
1580 SETUP::
1581
1582 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1583
1584 TESTS:
1585
1586 This gives the same answer as the slow, default method implemented
1587 in :class:`MatrixEuclideanJordanAlgebra`::
1588
1589 sage: set_random_seed()
1590 sage: J = ComplexHermitianEJA.random_instance()
1591 sage: x,y = J.random_elements(2)
1592 sage: Xe = x.to_matrix()
1593 sage: Ye = y.to_matrix()
1594 sage: X = ComplexHermitianEJA.real_unembed(Xe)
1595 sage: Y = ComplexHermitianEJA.real_unembed(Ye)
1596 sage: expected = (X*Y).trace().real()
1597 sage: actual = ComplexHermitianEJA.matrix_inner_product(Xe,Ye)
1598 sage: actual == expected
1599 True
1600
1601 """
1602 return RealMatrixEuclideanJordanAlgebra.matrix_inner_product(X,Y)/2
1603
1604
1605 class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra,
1606 ConcreteEuclideanJordanAlgebra):
1607 """
1608 The rank-n simple EJA consisting of complex Hermitian n-by-n
1609 matrices over the real numbers, the usual symmetric Jordan product,
1610 and the real-part-of-trace inner product. It has dimension `n^2` over
1611 the reals.
1612
1613 SETUP::
1614
1615 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1616
1617 EXAMPLES:
1618
1619 In theory, our "field" can be any subfield of the reals::
1620
1621 sage: ComplexHermitianEJA(2, RDF)
1622 Euclidean Jordan algebra of dimension 4 over Real Double Field
1623 sage: ComplexHermitianEJA(2, RR)
1624 Euclidean Jordan algebra of dimension 4 over Real Field with
1625 53 bits of precision
1626
1627 TESTS:
1628
1629 The dimension of this algebra is `n^2`::
1630
1631 sage: set_random_seed()
1632 sage: n_max = ComplexHermitianEJA._max_random_instance_size()
1633 sage: n = ZZ.random_element(1, n_max)
1634 sage: J = ComplexHermitianEJA(n)
1635 sage: J.dimension() == n^2
1636 True
1637
1638 The Jordan multiplication is what we think it is::
1639
1640 sage: set_random_seed()
1641 sage: J = ComplexHermitianEJA.random_instance()
1642 sage: x,y = J.random_elements(2)
1643 sage: actual = (x*y).to_matrix()
1644 sage: X = x.to_matrix()
1645 sage: Y = y.to_matrix()
1646 sage: expected = (X*Y + Y*X)/2
1647 sage: actual == expected
1648 True
1649 sage: J(expected) == x*y
1650 True
1651
1652 We can change the generator prefix::
1653
1654 sage: ComplexHermitianEJA(2, prefix='z').gens()
1655 (z0, z1, z2, z3)
1656
1657 We can construct the (trivial) algebra of rank zero::
1658
1659 sage: ComplexHermitianEJA(0)
1660 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1661
1662 """
1663
1664 @classmethod
1665 def _denormalized_basis(cls, n, field):
1666 """
1667 Returns a basis for the space of complex Hermitian n-by-n matrices.
1668
1669 Why do we embed these? Basically, because all of numerical linear
1670 algebra assumes that you're working with vectors consisting of `n`
1671 entries from a field and scalars from the same field. There's no way
1672 to tell SageMath that (for example) the vectors contain complex
1673 numbers, while the scalar field is real.
1674
1675 SETUP::
1676
1677 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1678
1679 TESTS::
1680
1681 sage: set_random_seed()
1682 sage: n = ZZ.random_element(1,5)
1683 sage: field = QuadraticField(2, 'sqrt2')
1684 sage: B = ComplexHermitianEJA._denormalized_basis(n, field)
1685 sage: all( M.is_symmetric() for M in B)
1686 True
1687
1688 """
1689 R = PolynomialRing(field, 'z')
1690 z = R.gen()
1691 F = field.extension(z**2 + 1, 'I')
1692 I = F.gen()
1693
1694 # This is like the symmetric case, but we need to be careful:
1695 #
1696 # * We want conjugate-symmetry, not just symmetry.
1697 # * The diagonal will (as a result) be real.
1698 #
1699 S = []
1700 for i in range(n):
1701 for j in range(i+1):
1702 Eij = matrix(F, n, lambda k,l: k==i and l==j)
1703 if i == j:
1704 Sij = cls.real_embed(Eij)
1705 S.append(Sij)
1706 else:
1707 # The second one has a minus because it's conjugated.
1708 Sij_real = cls.real_embed(Eij + Eij.transpose())
1709 S.append(Sij_real)
1710 Sij_imag = cls.real_embed(I*Eij - I*Eij.transpose())
1711 S.append(Sij_imag)
1712
1713 # Since we embedded these, we can drop back to the "field" that we
1714 # started with instead of the complex extension "F".
1715 return ( s.change_ring(field) for s in S )
1716
1717
1718 def __init__(self, n, field=AA, **kwargs):
1719 basis = self._denormalized_basis(n,field)
1720 super(ComplexHermitianEJA,self).__init__(field,
1721 basis,
1722 check_axioms=False,
1723 **kwargs)
1724 self.rank.set_cache(n)
1725
1726 @staticmethod
1727 def _max_random_instance_size():
1728 return 3 # Dimension 9
1729
1730 @classmethod
1731 def random_instance(cls, field=AA, **kwargs):
1732 """
1733 Return a random instance of this type of algebra.
1734 """
1735 n = ZZ.random_element(cls._max_random_instance_size() + 1)
1736 return cls(n, field, **kwargs)
1737
1738 class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
1739 @staticmethod
1740 def real_embed(M):
1741 """
1742 Embed the n-by-n quaternion matrix ``M`` into the space of real
1743 matrices of size 4n-by-4n by first sending each quaternion entry `z
1744 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
1745 c+di],[-c + di, a-bi]]`, and then embedding those into a real
1746 matrix.
1747
1748 SETUP::
1749
1750 sage: from mjo.eja.eja_algebra import \
1751 ....: QuaternionMatrixEuclideanJordanAlgebra
1752
1753 EXAMPLES::
1754
1755 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1756 sage: i,j,k = Q.gens()
1757 sage: x = 1 + 2*i + 3*j + 4*k
1758 sage: M = matrix(Q, 1, [[x]])
1759 sage: QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1760 [ 1 2 3 4]
1761 [-2 1 -4 3]
1762 [-3 4 1 -2]
1763 [-4 -3 2 1]
1764
1765 Embedding is a homomorphism (isomorphism, in fact)::
1766
1767 sage: set_random_seed()
1768 sage: n = ZZ.random_element(2)
1769 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1770 sage: X = random_matrix(Q, n)
1771 sage: Y = random_matrix(Q, n)
1772 sage: Xe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X)
1773 sage: Ye = QuaternionMatrixEuclideanJordanAlgebra.real_embed(Y)
1774 sage: XYe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1775 sage: Xe*Ye == XYe
1776 True
1777
1778 """
1779 quaternions = M.base_ring()
1780 n = M.nrows()
1781 if M.ncols() != n:
1782 raise ValueError("the matrix 'M' must be square")
1783
1784 F = QuadraticField(-1, 'I')
1785 i = F.gen()
1786
1787 blocks = []
1788 for z in M.list():
1789 t = z.coefficient_tuple()
1790 a = t[0]
1791 b = t[1]
1792 c = t[2]
1793 d = t[3]
1794 cplxM = matrix(F, 2, [[ a + b*i, c + d*i],
1795 [-c + d*i, a - b*i]])
1796 realM = ComplexMatrixEuclideanJordanAlgebra.real_embed(cplxM)
1797 blocks.append(realM)
1798
1799 # We should have real entries by now, so use the realest field
1800 # we've got for the return value.
1801 return matrix.block(quaternions.base_ring(), n, blocks)
1802
1803
1804
1805 @staticmethod
1806 def real_unembed(M):
1807 """
1808 The inverse of _embed_quaternion_matrix().
1809
1810 SETUP::
1811
1812 sage: from mjo.eja.eja_algebra import \
1813 ....: QuaternionMatrixEuclideanJordanAlgebra
1814
1815 EXAMPLES::
1816
1817 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
1818 ....: [-2, 1, -4, 3],
1819 ....: [-3, 4, 1, -2],
1820 ....: [-4, -3, 2, 1]])
1821 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(M)
1822 [1 + 2*i + 3*j + 4*k]
1823
1824 TESTS:
1825
1826 Unembedding is the inverse of embedding::
1827
1828 sage: set_random_seed()
1829 sage: Q = QuaternionAlgebra(QQ, -1, -1)
1830 sage: M = random_matrix(Q, 3)
1831 sage: Me = QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1832 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1833 True
1834
1835 """
1836 n = ZZ(M.nrows())
1837 if M.ncols() != n:
1838 raise ValueError("the matrix 'M' must be square")
1839 if not n.mod(4).is_zero():
1840 raise ValueError("the matrix 'M' must be a quaternion embedding")
1841
1842 # Use the base ring of the matrix to ensure that its entries can be
1843 # multiplied by elements of the quaternion algebra.
1844 field = M.base_ring()
1845 Q = QuaternionAlgebra(field,-1,-1)
1846 i,j,k = Q.gens()
1847
1848 # Go top-left to bottom-right (reading order), converting every
1849 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
1850 # quaternion block.
1851 elements = []
1852 for l in range(n/4):
1853 for m in range(n/4):
1854 submat = ComplexMatrixEuclideanJordanAlgebra.real_unembed(
1855 M[4*l:4*l+4,4*m:4*m+4] )
1856 if submat[0,0] != submat[1,1].conjugate():
1857 raise ValueError('bad on-diagonal submatrix')
1858 if submat[0,1] != -submat[1,0].conjugate():
1859 raise ValueError('bad off-diagonal submatrix')
1860 z = submat[0,0].real()
1861 z += submat[0,0].imag()*i
1862 z += submat[0,1].real()*j
1863 z += submat[0,1].imag()*k
1864 elements.append(z)
1865
1866 return matrix(Q, n/4, elements)
1867
1868
1869 @classmethod
1870 def matrix_inner_product(cls,X,Y):
1871 """
1872 Compute a matrix inner product in this algebra directly from
1873 its real embedding.
1874
1875 SETUP::
1876
1877 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1878
1879 TESTS:
1880
1881 This gives the same answer as the slow, default method implemented
1882 in :class:`MatrixEuclideanJordanAlgebra`::
1883
1884 sage: set_random_seed()
1885 sage: J = QuaternionHermitianEJA.random_instance()
1886 sage: x,y = J.random_elements(2)
1887 sage: Xe = x.to_matrix()
1888 sage: Ye = y.to_matrix()
1889 sage: X = QuaternionHermitianEJA.real_unembed(Xe)
1890 sage: Y = QuaternionHermitianEJA.real_unembed(Ye)
1891 sage: expected = (X*Y).trace().coefficient_tuple()[0]
1892 sage: actual = QuaternionHermitianEJA.matrix_inner_product(Xe,Ye)
1893 sage: actual == expected
1894 True
1895
1896 """
1897 return RealMatrixEuclideanJordanAlgebra.matrix_inner_product(X,Y)/4
1898
1899
1900 class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra,
1901 ConcreteEuclideanJordanAlgebra):
1902 r"""
1903 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
1904 matrices, the usual symmetric Jordan product, and the
1905 real-part-of-trace inner product. It has dimension `2n^2 - n` over
1906 the reals.
1907
1908 SETUP::
1909
1910 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1911
1912 EXAMPLES:
1913
1914 In theory, our "field" can be any subfield of the reals::
1915
1916 sage: QuaternionHermitianEJA(2, RDF)
1917 Euclidean Jordan algebra of dimension 6 over Real Double Field
1918 sage: QuaternionHermitianEJA(2, RR)
1919 Euclidean Jordan algebra of dimension 6 over Real Field with
1920 53 bits of precision
1921
1922 TESTS:
1923
1924 The dimension of this algebra is `2*n^2 - n`::
1925
1926 sage: set_random_seed()
1927 sage: n_max = QuaternionHermitianEJA._max_random_instance_size()
1928 sage: n = ZZ.random_element(1, n_max)
1929 sage: J = QuaternionHermitianEJA(n)
1930 sage: J.dimension() == 2*(n^2) - n
1931 True
1932
1933 The Jordan multiplication is what we think it is::
1934
1935 sage: set_random_seed()
1936 sage: J = QuaternionHermitianEJA.random_instance()
1937 sage: x,y = J.random_elements(2)
1938 sage: actual = (x*y).to_matrix()
1939 sage: X = x.to_matrix()
1940 sage: Y = y.to_matrix()
1941 sage: expected = (X*Y + Y*X)/2
1942 sage: actual == expected
1943 True
1944 sage: J(expected) == x*y
1945 True
1946
1947 We can change the generator prefix::
1948
1949 sage: QuaternionHermitianEJA(2, prefix='a').gens()
1950 (a0, a1, a2, a3, a4, a5)
1951
1952 We can construct the (trivial) algebra of rank zero::
1953
1954 sage: QuaternionHermitianEJA(0)
1955 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1956
1957 """
1958 @classmethod
1959 def _denormalized_basis(cls, n, field):
1960 """
1961 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
1962
1963 Why do we embed these? Basically, because all of numerical
1964 linear algebra assumes that you're working with vectors consisting
1965 of `n` entries from a field and scalars from the same field. There's
1966 no way to tell SageMath that (for example) the vectors contain
1967 complex numbers, while the scalar field is real.
1968
1969 SETUP::
1970
1971 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1972
1973 TESTS::
1974
1975 sage: set_random_seed()
1976 sage: n = ZZ.random_element(1,5)
1977 sage: B = QuaternionHermitianEJA._denormalized_basis(n,QQ)
1978 sage: all( M.is_symmetric() for M in B )
1979 True
1980
1981 """
1982 Q = QuaternionAlgebra(QQ,-1,-1)
1983 I,J,K = Q.gens()
1984
1985 # This is like the symmetric case, but we need to be careful:
1986 #
1987 # * We want conjugate-symmetry, not just symmetry.
1988 # * The diagonal will (as a result) be real.
1989 #
1990 S = []
1991 for i in range(n):
1992 for j in range(i+1):
1993 Eij = matrix(Q, n, lambda k,l: k==i and l==j)
1994 if i == j:
1995 Sij = cls.real_embed(Eij)
1996 S.append(Sij)
1997 else:
1998 # The second, third, and fourth ones have a minus
1999 # because they're conjugated.
2000 Sij_real = cls.real_embed(Eij + Eij.transpose())
2001 S.append(Sij_real)
2002 Sij_I = cls.real_embed(I*Eij - I*Eij.transpose())
2003 S.append(Sij_I)
2004 Sij_J = cls.real_embed(J*Eij - J*Eij.transpose())
2005 S.append(Sij_J)
2006 Sij_K = cls.real_embed(K*Eij - K*Eij.transpose())
2007 S.append(Sij_K)
2008
2009 # Since we embedded these, we can drop back to the "field" that we
2010 # started with instead of the quaternion algebra "Q".
2011 return ( s.change_ring(field) for s in S )
2012
2013
2014 def __init__(self, n, field=AA, **kwargs):
2015 basis = self._denormalized_basis(n,field)
2016 super(QuaternionHermitianEJA,self).__init__(field,
2017 basis,
2018 check_axioms=False,
2019 **kwargs)
2020 self.rank.set_cache(n)
2021
2022 @staticmethod
2023 def _max_random_instance_size():
2024 r"""
2025 The maximum rank of a random QuaternionHermitianEJA.
2026 """
2027 return 2 # Dimension 6
2028
2029 @classmethod
2030 def random_instance(cls, field=AA, **kwargs):
2031 """
2032 Return a random instance of this type of algebra.
2033 """
2034 n = ZZ.random_element(cls._max_random_instance_size() + 1)
2035 return cls(n, field, **kwargs)
2036
2037
2038 class HadamardEJA(RationalBasisEuclideanJordanAlgebra,
2039 ConcreteEuclideanJordanAlgebra):
2040 """
2041 Return the Euclidean Jordan Algebra corresponding to the set
2042 `R^n` under the Hadamard product.
2043
2044 Note: this is nothing more than the Cartesian product of ``n``
2045 copies of the spin algebra. Once Cartesian product algebras
2046 are implemented, this can go.
2047
2048 SETUP::
2049
2050 sage: from mjo.eja.eja_algebra import HadamardEJA
2051
2052 EXAMPLES:
2053
2054 This multiplication table can be verified by hand::
2055
2056 sage: J = HadamardEJA(3)
2057 sage: e0,e1,e2 = J.gens()
2058 sage: e0*e0
2059 e0
2060 sage: e0*e1
2061 0
2062 sage: e0*e2
2063 0
2064 sage: e1*e1
2065 e1
2066 sage: e1*e2
2067 0
2068 sage: e2*e2
2069 e2
2070
2071 TESTS:
2072
2073 We can change the generator prefix::
2074
2075 sage: HadamardEJA(3, prefix='r').gens()
2076 (r0, r1, r2)
2077
2078 """
2079 def __init__(self, n, field=AA, **kwargs):
2080 V = VectorSpace(field, n)
2081 mult_table = [ [ V.gen(i)*(i == j) for j in range(n) ]
2082 for i in range(n) ]
2083
2084 # Inner products are real numbers and not algebra
2085 # elements, so once we turn the algebra element
2086 # into a vector in inner_product(), we never go
2087 # back. As a result -- contrary to what we do with
2088 # self._multiplication_table -- we store the inner
2089 # product table as a plain old matrix and not as
2090 # an algebra operator.
2091 ip_table = matrix.identity(field,n)
2092 self._inner_product_matrix = ip_table
2093
2094 super(HadamardEJA, self).__init__(field,
2095 mult_table,
2096 check_axioms=False,
2097 **kwargs)
2098 self.rank.set_cache(n)
2099
2100 if n == 0:
2101 self.one.set_cache( self.zero() )
2102 else:
2103 self.one.set_cache( sum(self.gens()) )
2104
2105 @staticmethod
2106 def _max_random_instance_size():
2107 r"""
2108 The maximum dimension of a random HadamardEJA.
2109 """
2110 return 5
2111
2112 @classmethod
2113 def random_instance(cls, field=AA, **kwargs):
2114 """
2115 Return a random instance of this type of algebra.
2116 """
2117 n = ZZ.random_element(cls._max_random_instance_size() + 1)
2118 return cls(n, field, **kwargs)
2119
2120
2121 class BilinearFormEJA(RationalBasisEuclideanJordanAlgebra,
2122 ConcreteEuclideanJordanAlgebra):
2123 r"""
2124 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2125 with the half-trace inner product and jordan product ``x*y =
2126 (<Bx,y>,y_bar>, x0*y_bar + y0*x_bar)`` where `B = 1 \times B22` is
2127 a symmetric positive-definite "bilinear form" matrix. Its
2128 dimension is the size of `B`, and it has rank two in dimensions
2129 larger than two. It reduces to the ``JordanSpinEJA`` when `B` is
2130 the identity matrix of order ``n``.
2131
2132 We insist that the one-by-one upper-left identity block of `B` be
2133 passed in as well so that we can be passed a matrix of size zero
2134 to construct a trivial algebra.
2135
2136 SETUP::
2137
2138 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
2139 ....: JordanSpinEJA)
2140
2141 EXAMPLES:
2142
2143 When no bilinear form is specified, the identity matrix is used,
2144 and the resulting algebra is the Jordan spin algebra::
2145
2146 sage: B = matrix.identity(AA,3)
2147 sage: J0 = BilinearFormEJA(B)
2148 sage: J1 = JordanSpinEJA(3)
2149 sage: J0.multiplication_table() == J0.multiplication_table()
2150 True
2151
2152 An error is raised if the matrix `B` does not correspond to a
2153 positive-definite bilinear form::
2154
2155 sage: B = matrix.random(QQ,2,3)
2156 sage: J = BilinearFormEJA(B)
2157 Traceback (most recent call last):
2158 ...
2159 ValueError: bilinear form is not positive-definite
2160 sage: B = matrix.zero(QQ,3)
2161 sage: J = BilinearFormEJA(B)
2162 Traceback (most recent call last):
2163 ...
2164 ValueError: bilinear form is not positive-definite
2165
2166 TESTS:
2167
2168 We can create a zero-dimensional algebra::
2169
2170 sage: B = matrix.identity(AA,0)
2171 sage: J = BilinearFormEJA(B)
2172 sage: J.basis()
2173 Finite family {}
2174
2175 We can check the multiplication condition given in the Jordan, von
2176 Neumann, and Wigner paper (and also discussed on my "On the
2177 symmetry..." paper). Note that this relies heavily on the standard
2178 choice of basis, as does anything utilizing the bilinear form matrix::
2179
2180 sage: set_random_seed()
2181 sage: n = ZZ.random_element(5)
2182 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2183 sage: B11 = matrix.identity(QQ,1)
2184 sage: B22 = M.transpose()*M
2185 sage: B = block_matrix(2,2,[ [B11,0 ],
2186 ....: [0, B22 ] ])
2187 sage: J = BilinearFormEJA(B)
2188 sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
2189 sage: V = J.vector_space()
2190 sage: sis = [ J.from_vector(V([0] + (M.inverse()*ei).list()))
2191 ....: for ei in eis ]
2192 sage: actual = [ sis[i]*sis[j]
2193 ....: for i in range(n-1)
2194 ....: for j in range(n-1) ]
2195 sage: expected = [ J.one() if i == j else J.zero()
2196 ....: for i in range(n-1)
2197 ....: for j in range(n-1) ]
2198 sage: actual == expected
2199 True
2200 """
2201 def __init__(self, B, field=AA, **kwargs):
2202 n = B.nrows()
2203
2204 if not B.is_positive_definite():
2205 raise ValueError("bilinear form is not positive-definite")
2206
2207 V = VectorSpace(field, n)
2208 mult_table = [[V.zero() for j in range(n)] for i in range(n)]
2209 for i in range(n):
2210 for j in range(n):
2211 x = V.gen(i)
2212 y = V.gen(j)
2213 x0 = x[0]
2214 xbar = x[1:]
2215 y0 = y[0]
2216 ybar = y[1:]
2217 z0 = (B*x).inner_product(y)
2218 zbar = y0*xbar + x0*ybar
2219 z = V([z0] + zbar.list())
2220 mult_table[i][j] = z
2221
2222 # Inner products are real numbers and not algebra
2223 # elements, so once we turn the algebra element
2224 # into a vector in inner_product(), we never go
2225 # back. As a result -- contrary to what we do with
2226 # self._multiplication_table -- we store the inner
2227 # product table as a plain old matrix and not as
2228 # an algebra operator.
2229 ip_table = B
2230 self._inner_product_matrix = ip_table
2231
2232 super(BilinearFormEJA, self).__init__(field,
2233 mult_table,
2234 check_axioms=False,
2235 **kwargs)
2236
2237 # The rank of this algebra is two, unless we're in a
2238 # one-dimensional ambient space (because the rank is bounded
2239 # by the ambient dimension).
2240 self.rank.set_cache(min(n,2))
2241
2242 if n == 0:
2243 self.one.set_cache( self.zero() )
2244 else:
2245 self.one.set_cache( self.monomial(0) )
2246
2247 @staticmethod
2248 def _max_random_instance_size():
2249 r"""
2250 The maximum dimension of a random BilinearFormEJA.
2251 """
2252 return 5
2253
2254 @classmethod
2255 def random_instance(cls, field=AA, **kwargs):
2256 """
2257 Return a random instance of this algebra.
2258 """
2259 n = ZZ.random_element(cls._max_random_instance_size() + 1)
2260 if n.is_zero():
2261 B = matrix.identity(field, n)
2262 return cls(B, field, **kwargs)
2263
2264 B11 = matrix.identity(field,1)
2265 M = matrix.random(field, n-1)
2266 I = matrix.identity(field, n-1)
2267 alpha = field.zero()
2268 while alpha.is_zero():
2269 alpha = field.random_element().abs()
2270 B22 = M.transpose()*M + alpha*I
2271
2272 from sage.matrix.special import block_matrix
2273 B = block_matrix(2,2, [ [B11, ZZ(0) ],
2274 [ZZ(0), B22 ] ])
2275
2276 return cls(B, field, **kwargs)
2277
2278
2279 class JordanSpinEJA(BilinearFormEJA):
2280 """
2281 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2282 with the usual inner product and jordan product ``x*y =
2283 (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2284 the reals.
2285
2286 SETUP::
2287
2288 sage: from mjo.eja.eja_algebra import JordanSpinEJA
2289
2290 EXAMPLES:
2291
2292 This multiplication table can be verified by hand::
2293
2294 sage: J = JordanSpinEJA(4)
2295 sage: e0,e1,e2,e3 = J.gens()
2296 sage: e0*e0
2297 e0
2298 sage: e0*e1
2299 e1
2300 sage: e0*e2
2301 e2
2302 sage: e0*e3
2303 e3
2304 sage: e1*e2
2305 0
2306 sage: e1*e3
2307 0
2308 sage: e2*e3
2309 0
2310
2311 We can change the generator prefix::
2312
2313 sage: JordanSpinEJA(2, prefix='B').gens()
2314 (B0, B1)
2315
2316 TESTS:
2317
2318 Ensure that we have the usual inner product on `R^n`::
2319
2320 sage: set_random_seed()
2321 sage: J = JordanSpinEJA.random_instance()
2322 sage: x,y = J.random_elements(2)
2323 sage: actual = x.inner_product(y)
2324 sage: expected = x.to_vector().inner_product(y.to_vector())
2325 sage: actual == expected
2326 True
2327
2328 """
2329 def __init__(self, n, field=AA, **kwargs):
2330 # This is a special case of the BilinearFormEJA with the identity
2331 # matrix as its bilinear form.
2332 B = matrix.identity(field, n)
2333 super(JordanSpinEJA, self).__init__(B, field, **kwargs)
2334
2335 @staticmethod
2336 def _max_random_instance_size():
2337 r"""
2338 The maximum dimension of a random JordanSpinEJA.
2339 """
2340 return 5
2341
2342 @classmethod
2343 def random_instance(cls, field=AA, **kwargs):
2344 """
2345 Return a random instance of this type of algebra.
2346
2347 Needed here to override the implementation for ``BilinearFormEJA``.
2348 """
2349 n = ZZ.random_element(cls._max_random_instance_size() + 1)
2350 return cls(n, field, **kwargs)
2351
2352
2353 class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra,
2354 ConcreteEuclideanJordanAlgebra):
2355 """
2356 The trivial Euclidean Jordan algebra consisting of only a zero element.
2357
2358 SETUP::
2359
2360 sage: from mjo.eja.eja_algebra import TrivialEJA
2361
2362 EXAMPLES::
2363
2364 sage: J = TrivialEJA()
2365 sage: J.dimension()
2366 0
2367 sage: J.zero()
2368 0
2369 sage: J.one()
2370 0
2371 sage: 7*J.one()*12*J.one()
2372 0
2373 sage: J.one().inner_product(J.one())
2374 0
2375 sage: J.one().norm()
2376 0
2377 sage: J.one().subalgebra_generated_by()
2378 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2379 sage: J.rank()
2380 0
2381
2382 """
2383 def __init__(self, field=AA, **kwargs):
2384 mult_table = []
2385 self._inner_product_matrix = matrix(field,0)
2386 super(TrivialEJA, self).__init__(field,
2387 mult_table,
2388 check_axioms=False,
2389 **kwargs)
2390 # The rank is zero using my definition, namely the dimension of the
2391 # largest subalgebra generated by any element.
2392 self.rank.set_cache(0)
2393 self.one.set_cache( self.zero() )
2394
2395 @classmethod
2396 def random_instance(cls, field=AA, **kwargs):
2397 # We don't take a "size" argument so the superclass method is
2398 # inappropriate for us.
2399 return cls(field, **kwargs)
2400
2401 class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra):
2402 r"""
2403 The external (orthogonal) direct sum of two other Euclidean Jordan
2404 algebras. Essentially the Cartesian product of its two factors.
2405 Every Euclidean Jordan algebra decomposes into an orthogonal
2406 direct sum of simple Euclidean Jordan algebras, so no generality
2407 is lost by providing only this construction.
2408
2409 SETUP::
2410
2411 sage: from mjo.eja.eja_algebra import (random_eja,
2412 ....: HadamardEJA,
2413 ....: RealSymmetricEJA,
2414 ....: DirectSumEJA)
2415
2416 EXAMPLES::
2417
2418 sage: J1 = HadamardEJA(2)
2419 sage: J2 = RealSymmetricEJA(3)
2420 sage: J = DirectSumEJA(J1,J2)
2421 sage: J.dimension()
2422 8
2423 sage: J.rank()
2424 5
2425
2426 TESTS:
2427
2428 The external direct sum construction is only valid when the two factors
2429 have the same base ring; an error is raised otherwise::
2430
2431 sage: set_random_seed()
2432 sage: J1 = random_eja(AA)
2433 sage: J2 = random_eja(QQ)
2434 sage: J = DirectSumEJA(J1,J2)
2435 Traceback (most recent call last):
2436 ...
2437 ValueError: algebras must share the same base field
2438
2439 """
2440 def __init__(self, J1, J2, **kwargs):
2441 if J1.base_ring() != J2.base_ring():
2442 raise ValueError("algebras must share the same base field")
2443 field = J1.base_ring()
2444
2445 self._factors = (J1, J2)
2446 n1 = J1.dimension()
2447 n2 = J2.dimension()
2448 n = n1+n2
2449 V = VectorSpace(field, n)
2450 mult_table = [ [ V.zero() for j in range(n) ]
2451 for i in range(n) ]
2452 for i in range(n1):
2453 for j in range(n1):
2454 p = (J1.monomial(i)*J1.monomial(j)).to_vector()
2455 mult_table[i][j] = V(p.list() + [field.zero()]*n2)
2456
2457 for i in range(n2):
2458 for j in range(n2):
2459 p = (J2.monomial(i)*J2.monomial(j)).to_vector()
2460 mult_table[n1+i][n1+j] = V([field.zero()]*n1 + p.list())
2461
2462 super(DirectSumEJA, self).__init__(field,
2463 mult_table,
2464 check_axioms=False,
2465 **kwargs)
2466 self.rank.set_cache(J1.rank() + J2.rank())
2467
2468
2469 def factors(self):
2470 r"""
2471 Return the pair of this algebra's factors.
2472
2473 SETUP::
2474
2475 sage: from mjo.eja.eja_algebra import (HadamardEJA,
2476 ....: JordanSpinEJA,
2477 ....: DirectSumEJA)
2478
2479 EXAMPLES::
2480
2481 sage: J1 = HadamardEJA(2,QQ)
2482 sage: J2 = JordanSpinEJA(3,QQ)
2483 sage: J = DirectSumEJA(J1,J2)
2484 sage: J.factors()
2485 (Euclidean Jordan algebra of dimension 2 over Rational Field,
2486 Euclidean Jordan algebra of dimension 3 over Rational Field)
2487
2488 """
2489 return self._factors
2490
2491 def projections(self):
2492 r"""
2493 Return a pair of projections onto this algebra's factors.
2494
2495 SETUP::
2496
2497 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
2498 ....: ComplexHermitianEJA,
2499 ....: DirectSumEJA)
2500
2501 EXAMPLES::
2502
2503 sage: J1 = JordanSpinEJA(2)
2504 sage: J2 = ComplexHermitianEJA(2)
2505 sage: J = DirectSumEJA(J1,J2)
2506 sage: (pi_left, pi_right) = J.projections()
2507 sage: J.one().to_vector()
2508 (1, 0, 1, 0, 0, 1)
2509 sage: pi_left(J.one()).to_vector()
2510 (1, 0)
2511 sage: pi_right(J.one()).to_vector()
2512 (1, 0, 0, 1)
2513
2514 """
2515 (J1,J2) = self.factors()
2516 m = J1.dimension()
2517 n = J2.dimension()
2518 V_basis = self.vector_space().basis()
2519 # Need to specify the dimensions explicitly so that we don't
2520 # wind up with a zero-by-zero matrix when we want e.g. a
2521 # zero-by-two matrix (important for composing things).
2522 P1 = matrix(self.base_ring(), m, m+n, V_basis[:m])
2523 P2 = matrix(self.base_ring(), n, m+n, V_basis[m:])
2524 pi_left = FiniteDimensionalEuclideanJordanAlgebraOperator(self,J1,P1)
2525 pi_right = FiniteDimensionalEuclideanJordanAlgebraOperator(self,J2,P2)
2526 return (pi_left, pi_right)
2527
2528 def inclusions(self):
2529 r"""
2530 Return the pair of inclusion maps from our factors into us.
2531
2532 SETUP::
2533
2534 sage: from mjo.eja.eja_algebra import (random_eja,
2535 ....: JordanSpinEJA,
2536 ....: RealSymmetricEJA,
2537 ....: DirectSumEJA)
2538
2539 EXAMPLES::
2540
2541 sage: J1 = JordanSpinEJA(3)
2542 sage: J2 = RealSymmetricEJA(2)
2543 sage: J = DirectSumEJA(J1,J2)
2544 sage: (iota_left, iota_right) = J.inclusions()
2545 sage: iota_left(J1.zero()) == J.zero()
2546 True
2547 sage: iota_right(J2.zero()) == J.zero()
2548 True
2549 sage: J1.one().to_vector()
2550 (1, 0, 0)
2551 sage: iota_left(J1.one()).to_vector()
2552 (1, 0, 0, 0, 0, 0)
2553 sage: J2.one().to_vector()
2554 (1, 0, 1)
2555 sage: iota_right(J2.one()).to_vector()
2556 (0, 0, 0, 1, 0, 1)
2557 sage: J.one().to_vector()
2558 (1, 0, 0, 1, 0, 1)
2559
2560 TESTS:
2561
2562 Composing a projection with the corresponding inclusion should
2563 produce the identity map, and mismatching them should produce
2564 the zero map::
2565
2566 sage: set_random_seed()
2567 sage: J1 = random_eja()
2568 sage: J2 = random_eja()
2569 sage: J = DirectSumEJA(J1,J2)
2570 sage: (iota_left, iota_right) = J.inclusions()
2571 sage: (pi_left, pi_right) = J.projections()
2572 sage: pi_left*iota_left == J1.one().operator()
2573 True
2574 sage: pi_right*iota_right == J2.one().operator()
2575 True
2576 sage: (pi_left*iota_right).is_zero()
2577 True
2578 sage: (pi_right*iota_left).is_zero()
2579 True
2580
2581 """
2582 (J1,J2) = self.factors()
2583 m = J1.dimension()
2584 n = J2.dimension()
2585 V_basis = self.vector_space().basis()
2586 # Need to specify the dimensions explicitly so that we don't
2587 # wind up with a zero-by-zero matrix when we want e.g. a
2588 # two-by-zero matrix (important for composing things).
2589 I1 = matrix.column(self.base_ring(), m, m+n, V_basis[:m])
2590 I2 = matrix.column(self.base_ring(), n, m+n, V_basis[m:])
2591 iota_left = FiniteDimensionalEuclideanJordanAlgebraOperator(J1,self,I1)
2592 iota_right = FiniteDimensionalEuclideanJordanAlgebraOperator(J2,self,I2)
2593 return (iota_left, iota_right)
2594
2595 def inner_product(self, x, y):
2596 r"""
2597 The standard Cartesian inner-product.
2598
2599 We project ``x`` and ``y`` onto our factors, and add up the
2600 inner-products from the subalgebras.
2601
2602 SETUP::
2603
2604
2605 sage: from mjo.eja.eja_algebra import (HadamardEJA,
2606 ....: QuaternionHermitianEJA,
2607 ....: DirectSumEJA)
2608
2609 EXAMPLE::
2610
2611 sage: J1 = HadamardEJA(3,QQ)
2612 sage: J2 = QuaternionHermitianEJA(2,QQ,normalize_basis=False)
2613 sage: J = DirectSumEJA(J1,J2)
2614 sage: x1 = J1.one()
2615 sage: x2 = x1
2616 sage: y1 = J2.one()
2617 sage: y2 = y1
2618 sage: x1.inner_product(x2)
2619 3
2620 sage: y1.inner_product(y2)
2621 2
2622 sage: J.one().inner_product(J.one())
2623 5
2624
2625 """
2626 (pi_left, pi_right) = self.projections()
2627 x1 = pi_left(x)
2628 x2 = pi_right(x)
2629 y1 = pi_left(y)
2630 y2 = pi_right(y)
2631
2632 return (x1.inner_product(y1) + x2.inner_product(y2))
2633
2634
2635
2636 random_eja = ConcreteEuclideanJordanAlgebra.random_instance