]> gitweb.michael.orlitzky.com - sage.d.git/blob - mjo/eja/eja_algebra.py
eja: remove more completed TODO items.
[sage.d.git] / mjo / eja / eja_algebra.py
1 """
2 Representations and constructions for Euclidean Jordan algebras.
3
4 A Euclidean Jordan algebra is a Jordan algebra that has some
5 additional properties:
6
7 1. It is finite-dimensional.
8 2. Its scalar field is the real numbers.
9 3a. An inner product is defined on it, and...
10 3b. That inner product is compatible with the Jordan product
11 in the sense that `<x*y,z> = <y,x*z>` for all elements
12 `x,y,z` in the algebra.
13
14 Every Euclidean Jordan algebra is formally-real: for any two elements
15 `x` and `y` in the algebra, `x^{2} + y^{2} = 0` implies that `x = y =
16 0`. Conversely, every finite-dimensional formally-real Jordan algebra
17 can be made into a Euclidean Jordan algebra with an appropriate choice
18 of inner-product.
19
20 Formally-real Jordan algebras were originally studied as a framework
21 for quantum mechanics. Today, Euclidean Jordan algebras are crucial in
22 symmetric cone optimization, since every symmetric cone arises as the
23 cone of squares in some Euclidean Jordan algebra.
24
25 It is known that every Euclidean Jordan algebra decomposes into an
26 orthogonal direct sum (essentially, a Cartesian product) of simple
27 algebras, and that moreover, up to Jordan-algebra isomorphism, there
28 are only five families of simple algebras. We provide constructions
29 for these simple algebras:
30
31 * :class:`BilinearFormEJA`
32 * :class:`RealSymmetricEJA`
33 * :class:`ComplexHermitianEJA`
34 * :class:`QuaternionHermitianEJA`
35 * :class:`OctonionHermitianEJA`
36
37 In addition to these, we provide two other example constructions,
38
39 * :class:`JordanSpinEJA`
40 * :class:`HadamardEJA`
41 * :class:`AlbertEJA`
42 * :class:`TrivialEJA`
43
44 The Jordan spin algebra is a bilinear form algebra where the bilinear
45 form is the identity. The Hadamard EJA is simply a Cartesian product
46 of one-dimensional spin algebras. The Albert EJA is simply a special
47 case of the :class:`OctonionHermitianEJA` where the matrices are
48 three-by-three and the resulting space has dimension 27. And
49 last/least, the trivial EJA is exactly what you think it is; it could
50 also be obtained by constructing a dimension-zero instance of any of
51 the other algebras. Cartesian products of these are also supported
52 using the usual ``cartesian_product()`` function; as a result, we
53 support (up to isomorphism) all Euclidean Jordan algebras.
54
55 SETUP::
56
57 sage: from mjo.eja.eja_algebra import random_eja
58
59 EXAMPLES::
60
61 sage: random_eja()
62 Euclidean Jordan algebra of dimension...
63 """
64
65 from sage.algebras.quatalg.quaternion_algebra import QuaternionAlgebra
66 from sage.categories.magmatic_algebras import MagmaticAlgebras
67 from sage.categories.sets_cat import cartesian_product
68 from sage.combinat.free_module import CombinatorialFreeModule
69 from sage.matrix.constructor import matrix
70 from sage.matrix.matrix_space import MatrixSpace
71 from sage.misc.cachefunc import cached_method
72 from sage.misc.table import table
73 from sage.modules.free_module import FreeModule, VectorSpace
74 from sage.rings.all import (ZZ, QQ, AA, QQbar, RR, RLF, CLF,
75 PolynomialRing,
76 QuadraticField)
77 from mjo.eja.eja_element import FiniteDimensionalEJAElement
78 from mjo.eja.eja_operator import FiniteDimensionalEJAOperator
79 from mjo.eja.eja_utils import _all2list, _mat2vec
80
81 class FiniteDimensionalEJA(CombinatorialFreeModule):
82 r"""
83 A finite-dimensional Euclidean Jordan algebra.
84
85 INPUT:
86
87 - ``basis`` -- a tuple; a tuple of basis elements in "matrix
88 form," which must be the same form as the arguments to
89 ``jordan_product`` and ``inner_product``. In reality, "matrix
90 form" can be either vectors, matrices, or a Cartesian product
91 (ordered tuple) of vectors or matrices. All of these would
92 ideally be vector spaces in sage with no special-casing
93 needed; but in reality we turn vectors into column-matrices
94 and Cartesian products `(a,b)` into column matrices
95 `(a,b)^{T}` after converting `a` and `b` themselves.
96
97 - ``jordan_product`` -- a function; afunction of two ``basis``
98 elements (in matrix form) that returns their jordan product,
99 also in matrix form; this will be applied to ``basis`` to
100 compute a multiplication table for the algebra.
101
102 - ``inner_product`` -- a function; a function of two ``basis``
103 elements (in matrix form) that returns their inner
104 product. This will be applied to ``basis`` to compute an
105 inner-product table (basically a matrix) for this algebra.
106
107 - ``matrix_space`` -- the space that your matrix basis lives in,
108 or ``None`` (the default). So long as your basis does not have
109 length zero you can omit this. But in trivial algebras, it is
110 required.
111
112 - ``field`` -- a subfield of the reals (default: ``AA``); the scalar
113 field for the algebra.
114
115 - ``orthonormalize`` -- boolean (default: ``True``); whether or
116 not to orthonormalize the basis. Doing so is expensive and
117 generally rules out using the rationals as your ``field``, but
118 is required for spectral decompositions.
119
120 SETUP::
121
122 sage: from mjo.eja.eja_algebra import random_eja
123
124 TESTS:
125
126 We should compute that an element subalgebra is associative even
127 if we circumvent the element method::
128
129 sage: set_random_seed()
130 sage: J = random_eja(field=QQ,orthonormalize=False)
131 sage: x = J.random_element()
132 sage: A = x.subalgebra_generated_by(orthonormalize=False)
133 sage: basis = tuple(b.superalgebra_element() for b in A.basis())
134 sage: J.subalgebra(basis, orthonormalize=False).is_associative()
135 True
136 """
137 Element = FiniteDimensionalEJAElement
138
139 def __init__(self,
140 basis,
141 jordan_product,
142 inner_product,
143 field=AA,
144 matrix_space=None,
145 orthonormalize=True,
146 associative=None,
147 cartesian_product=False,
148 check_field=True,
149 check_axioms=True,
150 prefix="b"):
151
152 n = len(basis)
153
154 if check_field:
155 if not field.is_subring(RR):
156 # Note: this does return true for the real algebraic
157 # field, the rationals, and any quadratic field where
158 # we've specified a real embedding.
159 raise ValueError("scalar field is not real")
160
161 if check_axioms:
162 # Check commutativity of the Jordan and inner-products.
163 # This has to be done before we build the multiplication
164 # and inner-product tables/matrices, because we take
165 # advantage of symmetry in the process.
166 if not all( jordan_product(bi,bj) == jordan_product(bj,bi)
167 for bi in basis
168 for bj in basis ):
169 raise ValueError("Jordan product is not commutative")
170
171 if not all( inner_product(bi,bj) == inner_product(bj,bi)
172 for bi in basis
173 for bj in basis ):
174 raise ValueError("inner-product is not commutative")
175
176
177 category = MagmaticAlgebras(field).FiniteDimensional()
178 category = category.WithBasis().Unital().Commutative()
179
180 if n <= 1:
181 # All zero- and one-dimensional algebras are just the real
182 # numbers with (some positive multiples of) the usual
183 # multiplication as its Jordan and inner-product.
184 associative = True
185 if associative is None:
186 # We should figure it out. As with check_axioms, we have to do
187 # this without the help of the _jordan_product_is_associative()
188 # method because we need to know the category before we
189 # initialize the algebra.
190 associative = all( jordan_product(jordan_product(bi,bj),bk)
191 ==
192 jordan_product(bi,jordan_product(bj,bk))
193 for bi in basis
194 for bj in basis
195 for bk in basis)
196
197 if associative:
198 # Element subalgebras can take advantage of this.
199 category = category.Associative()
200 if cartesian_product:
201 # Use join() here because otherwise we only get the
202 # "Cartesian product of..." and not the things themselves.
203 category = category.join([category,
204 category.CartesianProducts()])
205
206 # Call the superclass constructor so that we can use its from_vector()
207 # method to build our multiplication table.
208 CombinatorialFreeModule.__init__(self,
209 field,
210 range(n),
211 prefix=prefix,
212 category=category,
213 bracket=False)
214
215 # Now comes all of the hard work. We'll be constructing an
216 # ambient vector space V that our (vectorized) basis lives in,
217 # as well as a subspace W of V spanned by those (vectorized)
218 # basis elements. The W-coordinates are the coefficients that
219 # we see in things like x = 1*b1 + 2*b2.
220 vector_basis = basis
221
222 degree = 0
223 if n > 0:
224 degree = len(_all2list(basis[0]))
225
226 # Build an ambient space that fits our matrix basis when
227 # written out as "long vectors."
228 V = VectorSpace(field, degree)
229
230 # The matrix that will hole the orthonormal -> unorthonormal
231 # coordinate transformation.
232 self._deortho_matrix = None
233
234 if orthonormalize:
235 # Save a copy of the un-orthonormalized basis for later.
236 # Convert it to ambient V (vector) coordinates while we're
237 # at it, because we'd have to do it later anyway.
238 deortho_vector_basis = tuple( V(_all2list(b)) for b in basis )
239
240 from mjo.eja.eja_utils import gram_schmidt
241 basis = tuple(gram_schmidt(basis, inner_product))
242
243 # Save the (possibly orthonormalized) matrix basis for
244 # later, as well as the space that its elements live in.
245 # In most cases we can deduce the matrix space, but when
246 # n == 0 (that is, there are no basis elements) we cannot.
247 self._matrix_basis = basis
248 if matrix_space is None:
249 self._matrix_space = self._matrix_basis[0].parent()
250 else:
251 self._matrix_space = matrix_space
252
253 # Now create the vector space for the algebra, which will have
254 # its own set of non-ambient coordinates (in terms of the
255 # supplied basis).
256 vector_basis = tuple( V(_all2list(b)) for b in basis )
257 W = V.span_of_basis( vector_basis, check=check_axioms)
258
259 if orthonormalize:
260 # Now "W" is the vector space of our algebra coordinates. The
261 # variables "X1", "X2",... refer to the entries of vectors in
262 # W. Thus to convert back and forth between the orthonormal
263 # coordinates and the given ones, we need to stick the original
264 # basis in W.
265 U = V.span_of_basis( deortho_vector_basis, check=check_axioms)
266 self._deortho_matrix = matrix( U.coordinate_vector(q)
267 for q in vector_basis )
268
269
270 # Now we actually compute the multiplication and inner-product
271 # tables/matrices using the possibly-orthonormalized basis.
272 self._inner_product_matrix = matrix.identity(field, n)
273 self._multiplication_table = [ [0 for j in range(i+1)]
274 for i in range(n) ]
275
276 # Note: the Jordan and inner-products are defined in terms
277 # of the ambient basis. It's important that their arguments
278 # are in ambient coordinates as well.
279 for i in range(n):
280 for j in range(i+1):
281 # ortho basis w.r.t. ambient coords
282 q_i = basis[i]
283 q_j = basis[j]
284
285 # The jordan product returns a matrixy answer, so we
286 # have to convert it to the algebra coordinates.
287 elt = jordan_product(q_i, q_j)
288 elt = W.coordinate_vector(V(_all2list(elt)))
289 self._multiplication_table[i][j] = self.from_vector(elt)
290
291 if not orthonormalize:
292 # If we're orthonormalizing the basis with respect
293 # to an inner-product, then the inner-product
294 # matrix with respect to the resulting basis is
295 # just going to be the identity.
296 ip = inner_product(q_i, q_j)
297 self._inner_product_matrix[i,j] = ip
298 self._inner_product_matrix[j,i] = ip
299
300 self._inner_product_matrix._cache = {'hermitian': True}
301 self._inner_product_matrix.set_immutable()
302
303 if check_axioms:
304 if not self._is_jordanian():
305 raise ValueError("Jordan identity does not hold")
306 if not self._inner_product_is_associative():
307 raise ValueError("inner product is not associative")
308
309
310 def _coerce_map_from_base_ring(self):
311 """
312 Disable the map from the base ring into the algebra.
313
314 Performing a nonsense conversion like this automatically
315 is counterpedagogical. The fallback is to try the usual
316 element constructor, which should also fail.
317
318 SETUP::
319
320 sage: from mjo.eja.eja_algebra import random_eja
321
322 TESTS::
323
324 sage: set_random_seed()
325 sage: J = random_eja()
326 sage: J(1)
327 Traceback (most recent call last):
328 ...
329 ValueError: not an element of this algebra
330
331 """
332 return None
333
334
335 def product_on_basis(self, i, j):
336 r"""
337 Returns the Jordan product of the `i` and `j`th basis elements.
338
339 This completely defines the Jordan product on the algebra, and
340 is used direclty by our superclass machinery to implement
341 :meth:`product`.
342
343 SETUP::
344
345 sage: from mjo.eja.eja_algebra import random_eja
346
347 TESTS::
348
349 sage: set_random_seed()
350 sage: J = random_eja()
351 sage: n = J.dimension()
352 sage: bi = J.zero()
353 sage: bj = J.zero()
354 sage: bi_bj = J.zero()*J.zero()
355 sage: if n > 0:
356 ....: i = ZZ.random_element(n)
357 ....: j = ZZ.random_element(n)
358 ....: bi = J.monomial(i)
359 ....: bj = J.monomial(j)
360 ....: bi_bj = J.product_on_basis(i,j)
361 sage: bi*bj == bi_bj
362 True
363
364 """
365 # We only stored the lower-triangular portion of the
366 # multiplication table.
367 if j <= i:
368 return self._multiplication_table[i][j]
369 else:
370 return self._multiplication_table[j][i]
371
372 def inner_product(self, x, y):
373 """
374 The inner product associated with this Euclidean Jordan algebra.
375
376 Defaults to the trace inner product, but can be overridden by
377 subclasses if they are sure that the necessary properties are
378 satisfied.
379
380 SETUP::
381
382 sage: from mjo.eja.eja_algebra import (random_eja,
383 ....: HadamardEJA,
384 ....: BilinearFormEJA)
385
386 EXAMPLES:
387
388 Our inner product is "associative," which means the following for
389 a symmetric bilinear form::
390
391 sage: set_random_seed()
392 sage: J = random_eja()
393 sage: x,y,z = J.random_elements(3)
394 sage: (x*y).inner_product(z) == y.inner_product(x*z)
395 True
396
397 TESTS:
398
399 Ensure that this is the usual inner product for the algebras
400 over `R^n`::
401
402 sage: set_random_seed()
403 sage: J = HadamardEJA.random_instance()
404 sage: x,y = J.random_elements(2)
405 sage: actual = x.inner_product(y)
406 sage: expected = x.to_vector().inner_product(y.to_vector())
407 sage: actual == expected
408 True
409
410 Ensure that this is one-half of the trace inner-product in a
411 BilinearFormEJA that isn't just the reals (when ``n`` isn't
412 one). This is in Faraut and Koranyi, and also my "On the
413 symmetry..." paper::
414
415 sage: set_random_seed()
416 sage: J = BilinearFormEJA.random_instance()
417 sage: n = J.dimension()
418 sage: x = J.random_element()
419 sage: y = J.random_element()
420 sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
421 True
422
423 """
424 B = self._inner_product_matrix
425 return (B*x.to_vector()).inner_product(y.to_vector())
426
427
428 def is_associative(self):
429 r"""
430 Return whether or not this algebra's Jordan product is associative.
431
432 SETUP::
433
434 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
435
436 EXAMPLES::
437
438 sage: J = ComplexHermitianEJA(3, field=QQ, orthonormalize=False)
439 sage: J.is_associative()
440 False
441 sage: x = sum(J.gens())
442 sage: A = x.subalgebra_generated_by(orthonormalize=False)
443 sage: A.is_associative()
444 True
445
446 """
447 return "Associative" in self.category().axioms()
448
449 def _is_commutative(self):
450 r"""
451 Whether or not this algebra's multiplication table is commutative.
452
453 This method should of course always return ``True``, unless
454 this algebra was constructed with ``check_axioms=False`` and
455 passed an invalid multiplication table.
456 """
457 return all( x*y == y*x for x in self.gens() for y in self.gens() )
458
459 def _is_jordanian(self):
460 r"""
461 Whether or not this algebra's multiplication table respects the
462 Jordan identity `(x^{2})(xy) = x(x^{2}y)`.
463
464 We only check one arrangement of `x` and `y`, so for a
465 ``True`` result to be truly true, you should also check
466 :meth:`_is_commutative`. This method should of course always
467 return ``True``, unless this algebra was constructed with
468 ``check_axioms=False`` and passed an invalid multiplication table.
469 """
470 return all( (self.monomial(i)**2)*(self.monomial(i)*self.monomial(j))
471 ==
472 (self.monomial(i))*((self.monomial(i)**2)*self.monomial(j))
473 for i in range(self.dimension())
474 for j in range(self.dimension()) )
475
476 def _jordan_product_is_associative(self):
477 r"""
478 Return whether or not this algebra's Jordan product is
479 associative; that is, whether or not `x*(y*z) = (x*y)*z`
480 for all `x,y,x`.
481
482 This method should agree with :meth:`is_associative` unless
483 you lied about the value of the ``associative`` parameter
484 when you constructed the algebra.
485
486 SETUP::
487
488 sage: from mjo.eja.eja_algebra import (random_eja,
489 ....: RealSymmetricEJA,
490 ....: ComplexHermitianEJA,
491 ....: QuaternionHermitianEJA)
492
493 EXAMPLES::
494
495 sage: J = RealSymmetricEJA(4, orthonormalize=False)
496 sage: J._jordan_product_is_associative()
497 False
498 sage: x = sum(J.gens())
499 sage: A = x.subalgebra_generated_by()
500 sage: A._jordan_product_is_associative()
501 True
502
503 ::
504
505 sage: J = ComplexHermitianEJA(2,field=QQ,orthonormalize=False)
506 sage: J._jordan_product_is_associative()
507 False
508 sage: x = sum(J.gens())
509 sage: A = x.subalgebra_generated_by(orthonormalize=False)
510 sage: A._jordan_product_is_associative()
511 True
512
513 ::
514
515 sage: J = QuaternionHermitianEJA(2)
516 sage: J._jordan_product_is_associative()
517 False
518 sage: x = sum(J.gens())
519 sage: A = x.subalgebra_generated_by()
520 sage: A._jordan_product_is_associative()
521 True
522
523 TESTS:
524
525 The values we've presupplied to the constructors agree with
526 the computation::
527
528 sage: set_random_seed()
529 sage: J = random_eja()
530 sage: J.is_associative() == J._jordan_product_is_associative()
531 True
532
533 """
534 R = self.base_ring()
535
536 # Used to check whether or not something is zero.
537 epsilon = R.zero()
538 if not R.is_exact():
539 # I don't know of any examples that make this magnitude
540 # necessary because I don't know how to make an
541 # associative algebra when the element subalgebra
542 # construction is unreliable (as it is over RDF; we can't
543 # find the degree of an element because we can't compute
544 # the rank of a matrix). But even multiplication of floats
545 # is non-associative, so *some* epsilon is needed... let's
546 # just take the one from _inner_product_is_associative?
547 epsilon = 1e-15
548
549 for i in range(self.dimension()):
550 for j in range(self.dimension()):
551 for k in range(self.dimension()):
552 x = self.monomial(i)
553 y = self.monomial(j)
554 z = self.monomial(k)
555 diff = (x*y)*z - x*(y*z)
556
557 if diff.norm() > epsilon:
558 return False
559
560 return True
561
562 def _inner_product_is_associative(self):
563 r"""
564 Return whether or not this algebra's inner product `B` is
565 associative; that is, whether or not `B(xy,z) = B(x,yz)`.
566
567 This method should of course always return ``True``, unless
568 this algebra was constructed with ``check_axioms=False`` and
569 passed an invalid Jordan or inner-product.
570 """
571 R = self.base_ring()
572
573 # Used to check whether or not something is zero.
574 epsilon = R.zero()
575 if not R.is_exact():
576 # This choice is sufficient to allow the construction of
577 # QuaternionHermitianEJA(2, field=RDF) with check_axioms=True.
578 epsilon = 1e-15
579
580 for i in range(self.dimension()):
581 for j in range(self.dimension()):
582 for k in range(self.dimension()):
583 x = self.monomial(i)
584 y = self.monomial(j)
585 z = self.monomial(k)
586 diff = (x*y).inner_product(z) - x.inner_product(y*z)
587
588 if diff.abs() > epsilon:
589 return False
590
591 return True
592
593 def _element_constructor_(self, elt):
594 """
595 Construct an element of this algebra from its vector or matrix
596 representation.
597
598 This gets called only after the parent element _call_ method
599 fails to find a coercion for the argument.
600
601 SETUP::
602
603 sage: from mjo.eja.eja_algebra import (random_eja,
604 ....: JordanSpinEJA,
605 ....: HadamardEJA,
606 ....: RealSymmetricEJA)
607
608 EXAMPLES:
609
610 The identity in `S^n` is converted to the identity in the EJA::
611
612 sage: J = RealSymmetricEJA(3)
613 sage: I = matrix.identity(QQ,3)
614 sage: J(I) == J.one()
615 True
616
617 This skew-symmetric matrix can't be represented in the EJA::
618
619 sage: J = RealSymmetricEJA(3)
620 sage: A = matrix(QQ,3, lambda i,j: i-j)
621 sage: J(A)
622 Traceback (most recent call last):
623 ...
624 ValueError: not an element of this algebra
625
626 Tuples work as well, provided that the matrix basis for the
627 algebra consists of them::
628
629 sage: J1 = HadamardEJA(3)
630 sage: J2 = RealSymmetricEJA(2)
631 sage: J = cartesian_product([J1,J2])
632 sage: J( (J1.matrix_basis()[1], J2.matrix_basis()[2]) )
633 b1 + b5
634
635 TESTS:
636
637 Ensure that we can convert any element back and forth
638 faithfully between its matrix and algebra representations::
639
640 sage: set_random_seed()
641 sage: J = random_eja()
642 sage: x = J.random_element()
643 sage: J(x.to_matrix()) == x
644 True
645
646 We cannot coerce elements between algebras just because their
647 matrix representations are compatible::
648
649 sage: J1 = HadamardEJA(3)
650 sage: J2 = JordanSpinEJA(3)
651 sage: J2(J1.one())
652 Traceback (most recent call last):
653 ...
654 ValueError: not an element of this algebra
655 sage: J1(J2.zero())
656 Traceback (most recent call last):
657 ...
658 ValueError: not an element of this algebra
659 """
660 msg = "not an element of this algebra"
661 if elt in self.base_ring():
662 # Ensure that no base ring -> algebra coercion is performed
663 # by this method. There's some stupidity in sage that would
664 # otherwise propagate to this method; for example, sage thinks
665 # that the integer 3 belongs to the space of 2-by-2 matrices.
666 raise ValueError(msg)
667
668 try:
669 # Try to convert a vector into a column-matrix...
670 elt = elt.column()
671 except (AttributeError, TypeError):
672 # and ignore failure, because we weren't really expecting
673 # a vector as an argument anyway.
674 pass
675
676 if elt not in self.matrix_space():
677 raise ValueError(msg)
678
679 # Thanks for nothing! Matrix spaces aren't vector spaces in
680 # Sage, so we have to figure out its matrix-basis coordinates
681 # ourselves. We use the basis space's ring instead of the
682 # element's ring because the basis space might be an algebraic
683 # closure whereas the base ring of the 3-by-3 identity matrix
684 # could be QQ instead of QQbar.
685 #
686 # And, we also have to handle Cartesian product bases (when
687 # the matrix basis consists of tuples) here. The "good news"
688 # is that we're already converting everything to long vectors,
689 # and that strategy works for tuples as well.
690 #
691 # We pass check=False because the matrix basis is "guaranteed"
692 # to be linearly independent... right? Ha ha.
693 elt = _all2list(elt)
694 V = VectorSpace(self.base_ring(), len(elt))
695 W = V.span_of_basis( (V(_all2list(s)) for s in self.matrix_basis()),
696 check=False)
697
698 try:
699 coords = W.coordinate_vector(V(elt))
700 except ArithmeticError: # vector is not in free module
701 raise ValueError(msg)
702
703 return self.from_vector(coords)
704
705 def _repr_(self):
706 """
707 Return a string representation of ``self``.
708
709 SETUP::
710
711 sage: from mjo.eja.eja_algebra import JordanSpinEJA
712
713 TESTS:
714
715 Ensure that it says what we think it says::
716
717 sage: JordanSpinEJA(2, field=AA)
718 Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
719 sage: JordanSpinEJA(3, field=RDF)
720 Euclidean Jordan algebra of dimension 3 over Real Double Field
721
722 """
723 fmt = "Euclidean Jordan algebra of dimension {} over {}"
724 return fmt.format(self.dimension(), self.base_ring())
725
726
727 @cached_method
728 def characteristic_polynomial_of(self):
729 """
730 Return the algebra's "characteristic polynomial of" function,
731 which is itself a multivariate polynomial that, when evaluated
732 at the coordinates of some algebra element, returns that
733 element's characteristic polynomial.
734
735 The resulting polynomial has `n+1` variables, where `n` is the
736 dimension of this algebra. The first `n` variables correspond to
737 the coordinates of an algebra element: when evaluated at the
738 coordinates of an algebra element with respect to a certain
739 basis, the result is a univariate polynomial (in the one
740 remaining variable ``t``), namely the characteristic polynomial
741 of that element.
742
743 SETUP::
744
745 sage: from mjo.eja.eja_algebra import JordanSpinEJA, TrivialEJA
746
747 EXAMPLES:
748
749 The characteristic polynomial in the spin algebra is given in
750 Alizadeh, Example 11.11::
751
752 sage: J = JordanSpinEJA(3)
753 sage: p = J.characteristic_polynomial_of(); p
754 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
755 sage: xvec = J.one().to_vector()
756 sage: p(*xvec)
757 t^2 - 2*t + 1
758
759 By definition, the characteristic polynomial is a monic
760 degree-zero polynomial in a rank-zero algebra. Note that
761 Cayley-Hamilton is indeed satisfied since the polynomial
762 ``1`` evaluates to the identity element of the algebra on
763 any argument::
764
765 sage: J = TrivialEJA()
766 sage: J.characteristic_polynomial_of()
767 1
768
769 """
770 r = self.rank()
771 n = self.dimension()
772
773 # The list of coefficient polynomials a_0, a_1, a_2, ..., a_(r-1).
774 a = self._charpoly_coefficients()
775
776 # We go to a bit of trouble here to reorder the
777 # indeterminates, so that it's easier to evaluate the
778 # characteristic polynomial at x's coordinates and get back
779 # something in terms of t, which is what we want.
780 S = PolynomialRing(self.base_ring(),'t')
781 t = S.gen(0)
782 if r > 0:
783 R = a[0].parent()
784 S = PolynomialRing(S, R.variable_names())
785 t = S(t)
786
787 return (t**r + sum( a[k]*(t**k) for k in range(r) ))
788
789 def coordinate_polynomial_ring(self):
790 r"""
791 The multivariate polynomial ring in which this algebra's
792 :meth:`characteristic_polynomial_of` lives.
793
794 SETUP::
795
796 sage: from mjo.eja.eja_algebra import (HadamardEJA,
797 ....: RealSymmetricEJA)
798
799 EXAMPLES::
800
801 sage: J = HadamardEJA(2)
802 sage: J.coordinate_polynomial_ring()
803 Multivariate Polynomial Ring in X1, X2...
804 sage: J = RealSymmetricEJA(3,field=QQ,orthonormalize=False)
805 sage: J.coordinate_polynomial_ring()
806 Multivariate Polynomial Ring in X1, X2, X3, X4, X5, X6...
807
808 """
809 var_names = tuple( "X%d" % z for z in range(1, self.dimension()+1) )
810 return PolynomialRing(self.base_ring(), var_names)
811
812 def inner_product(self, x, y):
813 """
814 The inner product associated with this Euclidean Jordan algebra.
815
816 Defaults to the trace inner product, but can be overridden by
817 subclasses if they are sure that the necessary properties are
818 satisfied.
819
820 SETUP::
821
822 sage: from mjo.eja.eja_algebra import (random_eja,
823 ....: HadamardEJA,
824 ....: BilinearFormEJA)
825
826 EXAMPLES:
827
828 Our inner product is "associative," which means the following for
829 a symmetric bilinear form::
830
831 sage: set_random_seed()
832 sage: J = random_eja()
833 sage: x,y,z = J.random_elements(3)
834 sage: (x*y).inner_product(z) == y.inner_product(x*z)
835 True
836
837 TESTS:
838
839 Ensure that this is the usual inner product for the algebras
840 over `R^n`::
841
842 sage: set_random_seed()
843 sage: J = HadamardEJA.random_instance()
844 sage: x,y = J.random_elements(2)
845 sage: actual = x.inner_product(y)
846 sage: expected = x.to_vector().inner_product(y.to_vector())
847 sage: actual == expected
848 True
849
850 Ensure that this is one-half of the trace inner-product in a
851 BilinearFormEJA that isn't just the reals (when ``n`` isn't
852 one). This is in Faraut and Koranyi, and also my "On the
853 symmetry..." paper::
854
855 sage: set_random_seed()
856 sage: J = BilinearFormEJA.random_instance()
857 sage: n = J.dimension()
858 sage: x = J.random_element()
859 sage: y = J.random_element()
860 sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
861 True
862 """
863 B = self._inner_product_matrix
864 return (B*x.to_vector()).inner_product(y.to_vector())
865
866
867 def is_trivial(self):
868 """
869 Return whether or not this algebra is trivial.
870
871 A trivial algebra contains only the zero element.
872
873 SETUP::
874
875 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
876 ....: TrivialEJA)
877
878 EXAMPLES::
879
880 sage: J = ComplexHermitianEJA(3)
881 sage: J.is_trivial()
882 False
883
884 ::
885
886 sage: J = TrivialEJA()
887 sage: J.is_trivial()
888 True
889
890 """
891 return self.dimension() == 0
892
893
894 def multiplication_table(self):
895 """
896 Return a visual representation of this algebra's multiplication
897 table (on basis elements).
898
899 SETUP::
900
901 sage: from mjo.eja.eja_algebra import JordanSpinEJA
902
903 EXAMPLES::
904
905 sage: J = JordanSpinEJA(4)
906 sage: J.multiplication_table()
907 +----++----+----+----+----+
908 | * || b0 | b1 | b2 | b3 |
909 +====++====+====+====+====+
910 | b0 || b0 | b1 | b2 | b3 |
911 +----++----+----+----+----+
912 | b1 || b1 | b0 | 0 | 0 |
913 +----++----+----+----+----+
914 | b2 || b2 | 0 | b0 | 0 |
915 +----++----+----+----+----+
916 | b3 || b3 | 0 | 0 | b0 |
917 +----++----+----+----+----+
918
919 """
920 n = self.dimension()
921 # Prepend the header row.
922 M = [["*"] + list(self.gens())]
923
924 # And to each subsequent row, prepend an entry that belongs to
925 # the left-side "header column."
926 M += [ [self.monomial(i)] + [ self.monomial(i)*self.monomial(j)
927 for j in range(n) ]
928 for i in range(n) ]
929
930 return table(M, header_row=True, header_column=True, frame=True)
931
932
933 def matrix_basis(self):
934 """
935 Return an (often more natural) representation of this algebras
936 basis as an ordered tuple of matrices.
937
938 Every finite-dimensional Euclidean Jordan Algebra is a, up to
939 Jordan isomorphism, a direct sum of five simple
940 algebras---four of which comprise Hermitian matrices. And the
941 last type of algebra can of course be thought of as `n`-by-`1`
942 column matrices (ambiguusly called column vectors) to avoid
943 special cases. As a result, matrices (and column vectors) are
944 a natural representation format for Euclidean Jordan algebra
945 elements.
946
947 But, when we construct an algebra from a basis of matrices,
948 those matrix representations are lost in favor of coordinate
949 vectors *with respect to* that basis. We could eventually
950 convert back if we tried hard enough, but having the original
951 representations handy is valuable enough that we simply store
952 them and return them from this method.
953
954 Why implement this for non-matrix algebras? Avoiding special
955 cases for the :class:`BilinearFormEJA` pays with simplicity in
956 its own right. But mainly, we would like to be able to assume
957 that elements of a :class:`CartesianProductEJA` can be displayed
958 nicely, without having to have special classes for direct sums
959 one of whose components was a matrix algebra.
960
961 SETUP::
962
963 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
964 ....: RealSymmetricEJA)
965
966 EXAMPLES::
967
968 sage: J = RealSymmetricEJA(2)
969 sage: J.basis()
970 Finite family {0: b0, 1: b1, 2: b2}
971 sage: J.matrix_basis()
972 (
973 [1 0] [ 0 0.7071067811865475?] [0 0]
974 [0 0], [0.7071067811865475? 0], [0 1]
975 )
976
977 ::
978
979 sage: J = JordanSpinEJA(2)
980 sage: J.basis()
981 Finite family {0: b0, 1: b1}
982 sage: J.matrix_basis()
983 (
984 [1] [0]
985 [0], [1]
986 )
987 """
988 return self._matrix_basis
989
990
991 def matrix_space(self):
992 """
993 Return the matrix space in which this algebra's elements live, if
994 we think of them as matrices (including column vectors of the
995 appropriate size).
996
997 "By default" this will be an `n`-by-`1` column-matrix space,
998 except when the algebra is trivial. There it's `n`-by-`n`
999 (where `n` is zero), to ensure that two elements of the matrix
1000 space (empty matrices) can be multiplied. For algebras of
1001 matrices, this returns the space in which their
1002 real embeddings live.
1003
1004 SETUP::
1005
1006 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
1007 ....: JordanSpinEJA,
1008 ....: QuaternionHermitianEJA,
1009 ....: TrivialEJA)
1010
1011 EXAMPLES:
1012
1013 By default, the matrix representation is just a column-matrix
1014 equivalent to the vector representation::
1015
1016 sage: J = JordanSpinEJA(3)
1017 sage: J.matrix_space()
1018 Full MatrixSpace of 3 by 1 dense matrices over Algebraic
1019 Real Field
1020
1021 The matrix representation in the trivial algebra is
1022 zero-by-zero instead of the usual `n`-by-one::
1023
1024 sage: J = TrivialEJA()
1025 sage: J.matrix_space()
1026 Full MatrixSpace of 0 by 0 dense matrices over Algebraic
1027 Real Field
1028
1029 The matrix space for complex/quaternion Hermitian matrix EJA
1030 is the space in which their real-embeddings live, not the
1031 original complex/quaternion matrix space::
1032
1033 sage: J = ComplexHermitianEJA(2,field=QQ,orthonormalize=False)
1034 sage: J.matrix_space()
1035 Module of 2 by 2 matrices with entries in Algebraic Field over
1036 the scalar ring Rational Field
1037 sage: J = QuaternionHermitianEJA(1,field=QQ,orthonormalize=False)
1038 sage: J.matrix_space()
1039 Module of 1 by 1 matrices with entries in Quaternion
1040 Algebra (-1, -1) with base ring Rational Field over
1041 the scalar ring Rational Field
1042
1043 """
1044 return self._matrix_space
1045
1046
1047 @cached_method
1048 def one(self):
1049 """
1050 Return the unit element of this algebra.
1051
1052 SETUP::
1053
1054 sage: from mjo.eja.eja_algebra import (HadamardEJA,
1055 ....: random_eja)
1056
1057 EXAMPLES:
1058
1059 We can compute unit element in the Hadamard EJA::
1060
1061 sage: J = HadamardEJA(5)
1062 sage: J.one()
1063 b0 + b1 + b2 + b3 + b4
1064
1065 The unit element in the Hadamard EJA is inherited in the
1066 subalgebras generated by its elements::
1067
1068 sage: J = HadamardEJA(5)
1069 sage: J.one()
1070 b0 + b1 + b2 + b3 + b4
1071 sage: x = sum(J.gens())
1072 sage: A = x.subalgebra_generated_by(orthonormalize=False)
1073 sage: A.one()
1074 c0
1075 sage: A.one().superalgebra_element()
1076 b0 + b1 + b2 + b3 + b4
1077
1078 TESTS:
1079
1080 The identity element acts like the identity, regardless of
1081 whether or not we orthonormalize::
1082
1083 sage: set_random_seed()
1084 sage: J = random_eja()
1085 sage: x = J.random_element()
1086 sage: J.one()*x == x and x*J.one() == x
1087 True
1088 sage: A = x.subalgebra_generated_by()
1089 sage: y = A.random_element()
1090 sage: A.one()*y == y and y*A.one() == y
1091 True
1092
1093 ::
1094
1095 sage: set_random_seed()
1096 sage: J = random_eja(field=QQ, orthonormalize=False)
1097 sage: x = J.random_element()
1098 sage: J.one()*x == x and x*J.one() == x
1099 True
1100 sage: A = x.subalgebra_generated_by(orthonormalize=False)
1101 sage: y = A.random_element()
1102 sage: A.one()*y == y and y*A.one() == y
1103 True
1104
1105 The matrix of the unit element's operator is the identity,
1106 regardless of the base field and whether or not we
1107 orthonormalize::
1108
1109 sage: set_random_seed()
1110 sage: J = random_eja()
1111 sage: actual = J.one().operator().matrix()
1112 sage: expected = matrix.identity(J.base_ring(), J.dimension())
1113 sage: actual == expected
1114 True
1115 sage: x = J.random_element()
1116 sage: A = x.subalgebra_generated_by()
1117 sage: actual = A.one().operator().matrix()
1118 sage: expected = matrix.identity(A.base_ring(), A.dimension())
1119 sage: actual == expected
1120 True
1121
1122 ::
1123
1124 sage: set_random_seed()
1125 sage: J = random_eja(field=QQ, orthonormalize=False)
1126 sage: actual = J.one().operator().matrix()
1127 sage: expected = matrix.identity(J.base_ring(), J.dimension())
1128 sage: actual == expected
1129 True
1130 sage: x = J.random_element()
1131 sage: A = x.subalgebra_generated_by(orthonormalize=False)
1132 sage: actual = A.one().operator().matrix()
1133 sage: expected = matrix.identity(A.base_ring(), A.dimension())
1134 sage: actual == expected
1135 True
1136
1137 Ensure that the cached unit element (often precomputed by
1138 hand) agrees with the computed one::
1139
1140 sage: set_random_seed()
1141 sage: J = random_eja()
1142 sage: cached = J.one()
1143 sage: J.one.clear_cache()
1144 sage: J.one() == cached
1145 True
1146
1147 ::
1148
1149 sage: set_random_seed()
1150 sage: J = random_eja(field=QQ, orthonormalize=False)
1151 sage: cached = J.one()
1152 sage: J.one.clear_cache()
1153 sage: J.one() == cached
1154 True
1155
1156 """
1157 # We can brute-force compute the matrices of the operators
1158 # that correspond to the basis elements of this algebra.
1159 # If some linear combination of those basis elements is the
1160 # algebra identity, then the same linear combination of
1161 # their matrices has to be the identity matrix.
1162 #
1163 # Of course, matrices aren't vectors in sage, so we have to
1164 # appeal to the "long vectors" isometry.
1165 oper_vecs = [ _mat2vec(g.operator().matrix()) for g in self.gens() ]
1166
1167 # Now we use basic linear algebra to find the coefficients,
1168 # of the matrices-as-vectors-linear-combination, which should
1169 # work for the original algebra basis too.
1170 A = matrix(self.base_ring(), oper_vecs)
1171
1172 # We used the isometry on the left-hand side already, but we
1173 # still need to do it for the right-hand side. Recall that we
1174 # wanted something that summed to the identity matrix.
1175 b = _mat2vec( matrix.identity(self.base_ring(), self.dimension()) )
1176
1177 # Now if there's an identity element in the algebra, this
1178 # should work. We solve on the left to avoid having to
1179 # transpose the matrix "A".
1180 return self.from_vector(A.solve_left(b))
1181
1182
1183 def peirce_decomposition(self, c):
1184 """
1185 The Peirce decomposition of this algebra relative to the
1186 idempotent ``c``.
1187
1188 In the future, this can be extended to a complete system of
1189 orthogonal idempotents.
1190
1191 INPUT:
1192
1193 - ``c`` -- an idempotent of this algebra.
1194
1195 OUTPUT:
1196
1197 A triple (J0, J5, J1) containing two subalgebras and one subspace
1198 of this algebra,
1199
1200 - ``J0`` -- the algebra on the eigenspace of ``c.operator()``
1201 corresponding to the eigenvalue zero.
1202
1203 - ``J5`` -- the eigenspace (NOT a subalgebra) of ``c.operator()``
1204 corresponding to the eigenvalue one-half.
1205
1206 - ``J1`` -- the algebra on the eigenspace of ``c.operator()``
1207 corresponding to the eigenvalue one.
1208
1209 These are the only possible eigenspaces for that operator, and this
1210 algebra is a direct sum of them. The spaces ``J0`` and ``J1`` are
1211 orthogonal, and are subalgebras of this algebra with the appropriate
1212 restrictions.
1213
1214 SETUP::
1215
1216 sage: from mjo.eja.eja_algebra import random_eja, RealSymmetricEJA
1217
1218 EXAMPLES:
1219
1220 The canonical example comes from the symmetric matrices, which
1221 decompose into diagonal and off-diagonal parts::
1222
1223 sage: J = RealSymmetricEJA(3)
1224 sage: C = matrix(QQ, [ [1,0,0],
1225 ....: [0,1,0],
1226 ....: [0,0,0] ])
1227 sage: c = J(C)
1228 sage: J0,J5,J1 = J.peirce_decomposition(c)
1229 sage: J0
1230 Euclidean Jordan algebra of dimension 1...
1231 sage: J5
1232 Vector space of degree 6 and dimension 2...
1233 sage: J1
1234 Euclidean Jordan algebra of dimension 3...
1235 sage: J0.one().to_matrix()
1236 [0 0 0]
1237 [0 0 0]
1238 [0 0 1]
1239 sage: orig_df = AA.options.display_format
1240 sage: AA.options.display_format = 'radical'
1241 sage: J.from_vector(J5.basis()[0]).to_matrix()
1242 [ 0 0 1/2*sqrt(2)]
1243 [ 0 0 0]
1244 [1/2*sqrt(2) 0 0]
1245 sage: J.from_vector(J5.basis()[1]).to_matrix()
1246 [ 0 0 0]
1247 [ 0 0 1/2*sqrt(2)]
1248 [ 0 1/2*sqrt(2) 0]
1249 sage: AA.options.display_format = orig_df
1250 sage: J1.one().to_matrix()
1251 [1 0 0]
1252 [0 1 0]
1253 [0 0 0]
1254
1255 TESTS:
1256
1257 Every algebra decomposes trivially with respect to its identity
1258 element::
1259
1260 sage: set_random_seed()
1261 sage: J = random_eja()
1262 sage: J0,J5,J1 = J.peirce_decomposition(J.one())
1263 sage: J0.dimension() == 0 and J5.dimension() == 0
1264 True
1265 sage: J1.superalgebra() == J and J1.dimension() == J.dimension()
1266 True
1267
1268 The decomposition is into eigenspaces, and its components are
1269 therefore necessarily orthogonal. Moreover, the identity
1270 elements in the two subalgebras are the projections onto their
1271 respective subspaces of the superalgebra's identity element::
1272
1273 sage: set_random_seed()
1274 sage: J = random_eja()
1275 sage: x = J.random_element()
1276 sage: if not J.is_trivial():
1277 ....: while x.is_nilpotent():
1278 ....: x = J.random_element()
1279 sage: c = x.subalgebra_idempotent()
1280 sage: J0,J5,J1 = J.peirce_decomposition(c)
1281 sage: ipsum = 0
1282 sage: for (w,y,z) in zip(J0.basis(), J5.basis(), J1.basis()):
1283 ....: w = w.superalgebra_element()
1284 ....: y = J.from_vector(y)
1285 ....: z = z.superalgebra_element()
1286 ....: ipsum += w.inner_product(y).abs()
1287 ....: ipsum += w.inner_product(z).abs()
1288 ....: ipsum += y.inner_product(z).abs()
1289 sage: ipsum
1290 0
1291 sage: J1(c) == J1.one()
1292 True
1293 sage: J0(J.one() - c) == J0.one()
1294 True
1295
1296 """
1297 if not c.is_idempotent():
1298 raise ValueError("element is not idempotent: %s" % c)
1299
1300 # Default these to what they should be if they turn out to be
1301 # trivial, because eigenspaces_left() won't return eigenvalues
1302 # corresponding to trivial spaces (e.g. it returns only the
1303 # eigenspace corresponding to lambda=1 if you take the
1304 # decomposition relative to the identity element).
1305 trivial = self.subalgebra(())
1306 J0 = trivial # eigenvalue zero
1307 J5 = VectorSpace(self.base_ring(), 0) # eigenvalue one-half
1308 J1 = trivial # eigenvalue one
1309
1310 for (eigval, eigspace) in c.operator().matrix().right_eigenspaces():
1311 if eigval == ~(self.base_ring()(2)):
1312 J5 = eigspace
1313 else:
1314 gens = tuple( self.from_vector(b) for b in eigspace.basis() )
1315 subalg = self.subalgebra(gens, check_axioms=False)
1316 if eigval == 0:
1317 J0 = subalg
1318 elif eigval == 1:
1319 J1 = subalg
1320 else:
1321 raise ValueError("unexpected eigenvalue: %s" % eigval)
1322
1323 return (J0, J5, J1)
1324
1325
1326 def random_element(self, thorough=False):
1327 r"""
1328 Return a random element of this algebra.
1329
1330 Our algebra superclass method only returns a linear
1331 combination of at most two basis elements. We instead
1332 want the vector space "random element" method that
1333 returns a more diverse selection.
1334
1335 INPUT:
1336
1337 - ``thorough`` -- (boolean; default False) whether or not we
1338 should generate irrational coefficients for the random
1339 element when our base ring is irrational; this slows the
1340 algebra operations to a crawl, but any truly random method
1341 should include them
1342
1343 """
1344 # For a general base ring... maybe we can trust this to do the
1345 # right thing? Unlikely, but.
1346 V = self.vector_space()
1347 v = V.random_element()
1348
1349 if self.base_ring() is AA:
1350 # The "random element" method of the algebraic reals is
1351 # stupid at the moment, and only returns integers between
1352 # -2 and 2, inclusive:
1353 #
1354 # https://trac.sagemath.org/ticket/30875
1355 #
1356 # Instead, we implement our own "random vector" method,
1357 # and then coerce that into the algebra. We use the vector
1358 # space degree here instead of the dimension because a
1359 # subalgebra could (for example) be spanned by only two
1360 # vectors, each with five coordinates. We need to
1361 # generate all five coordinates.
1362 if thorough:
1363 v *= QQbar.random_element().real()
1364 else:
1365 v *= QQ.random_element()
1366
1367 return self.from_vector(V.coordinate_vector(v))
1368
1369 def random_elements(self, count, thorough=False):
1370 """
1371 Return ``count`` random elements as a tuple.
1372
1373 INPUT:
1374
1375 - ``thorough`` -- (boolean; default False) whether or not we
1376 should generate irrational coefficients for the random
1377 elements when our base ring is irrational; this slows the
1378 algebra operations to a crawl, but any truly random method
1379 should include them
1380
1381 SETUP::
1382
1383 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1384
1385 EXAMPLES::
1386
1387 sage: J = JordanSpinEJA(3)
1388 sage: x,y,z = J.random_elements(3)
1389 sage: all( [ x in J, y in J, z in J ])
1390 True
1391 sage: len( J.random_elements(10) ) == 10
1392 True
1393
1394 """
1395 return tuple( self.random_element(thorough)
1396 for idx in range(count) )
1397
1398
1399 @cached_method
1400 def _charpoly_coefficients(self):
1401 r"""
1402 The `r` polynomial coefficients of the "characteristic polynomial
1403 of" function.
1404
1405 SETUP::
1406
1407 sage: from mjo.eja.eja_algebra import random_eja
1408
1409 TESTS:
1410
1411 The theory shows that these are all homogeneous polynomials of
1412 a known degree::
1413
1414 sage: set_random_seed()
1415 sage: J = random_eja()
1416 sage: all(p.is_homogeneous() for p in J._charpoly_coefficients())
1417 True
1418
1419 """
1420 n = self.dimension()
1421 R = self.coordinate_polynomial_ring()
1422 vars = R.gens()
1423 F = R.fraction_field()
1424
1425 def L_x_i_j(i,j):
1426 # From a result in my book, these are the entries of the
1427 # basis representation of L_x.
1428 return sum( vars[k]*self.monomial(k).operator().matrix()[i,j]
1429 for k in range(n) )
1430
1431 L_x = matrix(F, n, n, L_x_i_j)
1432
1433 r = None
1434 if self.rank.is_in_cache():
1435 r = self.rank()
1436 # There's no need to pad the system with redundant
1437 # columns if we *know* they'll be redundant.
1438 n = r
1439
1440 # Compute an extra power in case the rank is equal to
1441 # the dimension (otherwise, we would stop at x^(r-1)).
1442 x_powers = [ (L_x**k)*self.one().to_vector()
1443 for k in range(n+1) ]
1444 A = matrix.column(F, x_powers[:n])
1445 AE = A.extended_echelon_form()
1446 E = AE[:,n:]
1447 A_rref = AE[:,:n]
1448 if r is None:
1449 r = A_rref.rank()
1450 b = x_powers[r]
1451
1452 # The theory says that only the first "r" coefficients are
1453 # nonzero, and they actually live in the original polynomial
1454 # ring and not the fraction field. We negate them because in
1455 # the actual characteristic polynomial, they get moved to the
1456 # other side where x^r lives. We don't bother to trim A_rref
1457 # down to a square matrix and solve the resulting system,
1458 # because the upper-left r-by-r portion of A_rref is
1459 # guaranteed to be the identity matrix, so e.g.
1460 #
1461 # A_rref.solve_right(Y)
1462 #
1463 # would just be returning Y.
1464 return (-E*b)[:r].change_ring(R)
1465
1466 @cached_method
1467 def rank(self):
1468 r"""
1469 Return the rank of this EJA.
1470
1471 This is a cached method because we know the rank a priori for
1472 all of the algebras we can construct. Thus we can avoid the
1473 expensive ``_charpoly_coefficients()`` call unless we truly
1474 need to compute the whole characteristic polynomial.
1475
1476 SETUP::
1477
1478 sage: from mjo.eja.eja_algebra import (HadamardEJA,
1479 ....: JordanSpinEJA,
1480 ....: RealSymmetricEJA,
1481 ....: ComplexHermitianEJA,
1482 ....: QuaternionHermitianEJA,
1483 ....: random_eja)
1484
1485 EXAMPLES:
1486
1487 The rank of the Jordan spin algebra is always two::
1488
1489 sage: JordanSpinEJA(2).rank()
1490 2
1491 sage: JordanSpinEJA(3).rank()
1492 2
1493 sage: JordanSpinEJA(4).rank()
1494 2
1495
1496 The rank of the `n`-by-`n` Hermitian real, complex, or
1497 quaternion matrices is `n`::
1498
1499 sage: RealSymmetricEJA(4).rank()
1500 4
1501 sage: ComplexHermitianEJA(3).rank()
1502 3
1503 sage: QuaternionHermitianEJA(2).rank()
1504 2
1505
1506 TESTS:
1507
1508 Ensure that every EJA that we know how to construct has a
1509 positive integer rank, unless the algebra is trivial in
1510 which case its rank will be zero::
1511
1512 sage: set_random_seed()
1513 sage: J = random_eja()
1514 sage: r = J.rank()
1515 sage: r in ZZ
1516 True
1517 sage: r > 0 or (r == 0 and J.is_trivial())
1518 True
1519
1520 Ensure that computing the rank actually works, since the ranks
1521 of all simple algebras are known and will be cached by default::
1522
1523 sage: set_random_seed() # long time
1524 sage: J = random_eja() # long time
1525 sage: cached = J.rank() # long time
1526 sage: J.rank.clear_cache() # long time
1527 sage: J.rank() == cached # long time
1528 True
1529
1530 """
1531 return len(self._charpoly_coefficients())
1532
1533
1534 def subalgebra(self, basis, **kwargs):
1535 r"""
1536 Create a subalgebra of this algebra from the given basis.
1537 """
1538 from mjo.eja.eja_subalgebra import FiniteDimensionalEJASubalgebra
1539 return FiniteDimensionalEJASubalgebra(self, basis, **kwargs)
1540
1541
1542 def vector_space(self):
1543 """
1544 Return the vector space that underlies this algebra.
1545
1546 SETUP::
1547
1548 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1549
1550 EXAMPLES::
1551
1552 sage: J = RealSymmetricEJA(2)
1553 sage: J.vector_space()
1554 Vector space of dimension 3 over...
1555
1556 """
1557 return self.zero().to_vector().parent().ambient_vector_space()
1558
1559
1560
1561 class RationalBasisEJA(FiniteDimensionalEJA):
1562 r"""
1563 Algebras whose supplied basis elements have all rational entries.
1564
1565 SETUP::
1566
1567 sage: from mjo.eja.eja_algebra import BilinearFormEJA
1568
1569 EXAMPLES:
1570
1571 The supplied basis is orthonormalized by default::
1572
1573 sage: B = matrix(QQ, [[1, 0, 0], [0, 25, -32], [0, -32, 41]])
1574 sage: J = BilinearFormEJA(B)
1575 sage: J.matrix_basis()
1576 (
1577 [1] [ 0] [ 0]
1578 [0] [1/5] [32/5]
1579 [0], [ 0], [ 5]
1580 )
1581
1582 """
1583 def __init__(self,
1584 basis,
1585 jordan_product,
1586 inner_product,
1587 field=AA,
1588 check_field=True,
1589 **kwargs):
1590
1591 if check_field:
1592 # Abuse the check_field parameter to check that the entries of
1593 # out basis (in ambient coordinates) are in the field QQ.
1594 # Use _all2list to get the vector coordinates of octonion
1595 # entries and not the octonions themselves (which are not
1596 # rational).
1597 if not all( all(b_i in QQ for b_i in _all2list(b))
1598 for b in basis ):
1599 raise TypeError("basis not rational")
1600
1601 super().__init__(basis,
1602 jordan_product,
1603 inner_product,
1604 field=field,
1605 check_field=check_field,
1606 **kwargs)
1607
1608 self._rational_algebra = None
1609 if field is not QQ:
1610 # There's no point in constructing the extra algebra if this
1611 # one is already rational.
1612 #
1613 # Note: the same Jordan and inner-products work here,
1614 # because they are necessarily defined with respect to
1615 # ambient coordinates and not any particular basis.
1616 self._rational_algebra = FiniteDimensionalEJA(
1617 basis,
1618 jordan_product,
1619 inner_product,
1620 field=QQ,
1621 matrix_space=self.matrix_space(),
1622 associative=self.is_associative(),
1623 orthonormalize=False,
1624 check_field=False,
1625 check_axioms=False)
1626
1627 @cached_method
1628 def _charpoly_coefficients(self):
1629 r"""
1630 SETUP::
1631
1632 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
1633 ....: JordanSpinEJA)
1634
1635 EXAMPLES:
1636
1637 The base ring of the resulting polynomial coefficients is what
1638 it should be, and not the rationals (unless the algebra was
1639 already over the rationals)::
1640
1641 sage: J = JordanSpinEJA(3)
1642 sage: J._charpoly_coefficients()
1643 (X1^2 - X2^2 - X3^2, -2*X1)
1644 sage: a0 = J._charpoly_coefficients()[0]
1645 sage: J.base_ring()
1646 Algebraic Real Field
1647 sage: a0.base_ring()
1648 Algebraic Real Field
1649
1650 """
1651 if self._rational_algebra is None:
1652 # There's no need to construct *another* algebra over the
1653 # rationals if this one is already over the
1654 # rationals. Likewise, if we never orthonormalized our
1655 # basis, we might as well just use the given one.
1656 return super()._charpoly_coefficients()
1657
1658 # Do the computation over the rationals. The answer will be
1659 # the same, because all we've done is a change of basis.
1660 # Then, change back from QQ to our real base ring
1661 a = ( a_i.change_ring(self.base_ring())
1662 for a_i in self._rational_algebra._charpoly_coefficients() )
1663
1664 if self._deortho_matrix is None:
1665 # This can happen if our base ring was, say, AA and we
1666 # chose not to (or didn't need to) orthonormalize. It's
1667 # still faster to do the computations over QQ even if
1668 # the numbers in the boxes stay the same.
1669 return tuple(a)
1670
1671 # Otherwise, convert the coordinate variables back to the
1672 # deorthonormalized ones.
1673 R = self.coordinate_polynomial_ring()
1674 from sage.modules.free_module_element import vector
1675 X = vector(R, R.gens())
1676 BX = self._deortho_matrix*X
1677
1678 subs_dict = { X[i]: BX[i] for i in range(len(X)) }
1679 return tuple( a_i.subs(subs_dict) for a_i in a )
1680
1681 class ConcreteEJA(FiniteDimensionalEJA):
1682 r"""
1683 A class for the Euclidean Jordan algebras that we know by name.
1684
1685 These are the Jordan algebras whose basis, multiplication table,
1686 rank, and so on are known a priori. More to the point, they are
1687 the Euclidean Jordan algebras for which we are able to conjure up
1688 a "random instance."
1689
1690 SETUP::
1691
1692 sage: from mjo.eja.eja_algebra import ConcreteEJA
1693
1694 TESTS:
1695
1696 Our basis is normalized with respect to the algebra's inner
1697 product, unless we specify otherwise::
1698
1699 sage: set_random_seed()
1700 sage: J = ConcreteEJA.random_instance()
1701 sage: all( b.norm() == 1 for b in J.gens() )
1702 True
1703
1704 Since our basis is orthonormal with respect to the algebra's inner
1705 product, and since we know that this algebra is an EJA, any
1706 left-multiplication operator's matrix will be symmetric because
1707 natural->EJA basis representation is an isometry and within the
1708 EJA the operator is self-adjoint by the Jordan axiom::
1709
1710 sage: set_random_seed()
1711 sage: J = ConcreteEJA.random_instance()
1712 sage: x = J.random_element()
1713 sage: x.operator().is_self_adjoint()
1714 True
1715 """
1716
1717 @staticmethod
1718 def _max_random_instance_dimension():
1719 r"""
1720 The maximum dimension of any random instance. Ten dimensions seems
1721 to be about the point where everything takes a turn for the
1722 worse. And dimension ten (but not nine) allows the 4-by-4 real
1723 Hermitian matrices, the 2-by-2 quaternion Hermitian matrices,
1724 and the 2-by-2 octonion Hermitian matrices.
1725 """
1726 return 10
1727
1728 @staticmethod
1729 def _max_random_instance_size(max_dimension):
1730 """
1731 Return an integer "size" that is an upper bound on the size of
1732 this algebra when it is used in a random test case. This size
1733 (which can be passed to the algebra's constructor) is itself
1734 based on the ``max_dimension`` parameter.
1735
1736 This method must be implemented in each subclass.
1737 """
1738 raise NotImplementedError
1739
1740 @classmethod
1741 def random_instance(cls, max_dimension=None, *args, **kwargs):
1742 """
1743 Return a random instance of this type of algebra whose dimension
1744 is less than or equal to ``max_dimension``. If the dimension bound
1745 is omitted, then the ``_max_random_instance_dimension()`` is used
1746 to get a suitable bound.
1747
1748 This method should be implemented in each subclass.
1749 """
1750 from sage.misc.prandom import choice
1751 eja_class = choice(cls.__subclasses__())
1752
1753 # These all bubble up to the RationalBasisEJA superclass
1754 # constructor, so any (kw)args valid there are also valid
1755 # here.
1756 return eja_class.random_instance(*args, **kwargs)
1757
1758
1759 class MatrixEJA(FiniteDimensionalEJA):
1760 @staticmethod
1761 def _denormalized_basis(A):
1762 """
1763 Returns a basis for the space of complex Hermitian n-by-n matrices.
1764
1765 Why do we embed these? Basically, because all of numerical linear
1766 algebra assumes that you're working with vectors consisting of `n`
1767 entries from a field and scalars from the same field. There's no way
1768 to tell SageMath that (for example) the vectors contain complex
1769 numbers, while the scalar field is real.
1770
1771 SETUP::
1772
1773 sage: from mjo.hurwitz import (ComplexMatrixAlgebra,
1774 ....: QuaternionMatrixAlgebra,
1775 ....: OctonionMatrixAlgebra)
1776 sage: from mjo.eja.eja_algebra import MatrixEJA
1777
1778 TESTS::
1779
1780 sage: set_random_seed()
1781 sage: n = ZZ.random_element(1,5)
1782 sage: A = MatrixSpace(QQ, n)
1783 sage: B = MatrixEJA._denormalized_basis(A)
1784 sage: all( M.is_hermitian() for M in B)
1785 True
1786
1787 ::
1788
1789 sage: set_random_seed()
1790 sage: n = ZZ.random_element(1,5)
1791 sage: A = ComplexMatrixAlgebra(n, scalars=QQ)
1792 sage: B = MatrixEJA._denormalized_basis(A)
1793 sage: all( M.is_hermitian() for M in B)
1794 True
1795
1796 ::
1797
1798 sage: set_random_seed()
1799 sage: n = ZZ.random_element(1,5)
1800 sage: A = QuaternionMatrixAlgebra(n, scalars=QQ)
1801 sage: B = MatrixEJA._denormalized_basis(A)
1802 sage: all( M.is_hermitian() for M in B )
1803 True
1804
1805 ::
1806
1807 sage: set_random_seed()
1808 sage: n = ZZ.random_element(1,5)
1809 sage: A = OctonionMatrixAlgebra(n, scalars=QQ)
1810 sage: B = MatrixEJA._denormalized_basis(A)
1811 sage: all( M.is_hermitian() for M in B )
1812 True
1813
1814 """
1815 # These work for real MatrixSpace, whose monomials only have
1816 # two coordinates (because the last one would always be "1").
1817 es = A.base_ring().gens()
1818 gen = lambda A,m: A.monomial(m[:2])
1819
1820 if hasattr(A, 'entry_algebra_gens'):
1821 # We've got a MatrixAlgebra, and its monomials will have
1822 # three coordinates.
1823 es = A.entry_algebra_gens()
1824 gen = lambda A,m: A.monomial(m)
1825
1826 basis = []
1827 for i in range(A.nrows()):
1828 for j in range(i+1):
1829 if i == j:
1830 E_ii = gen(A, (i,j,es[0]))
1831 basis.append(E_ii)
1832 else:
1833 for e in es:
1834 E_ij = gen(A, (i,j,e))
1835 E_ij += E_ij.conjugate_transpose()
1836 basis.append(E_ij)
1837
1838 return tuple( basis )
1839
1840 @staticmethod
1841 def jordan_product(X,Y):
1842 return (X*Y + Y*X)/2
1843
1844 @staticmethod
1845 def trace_inner_product(X,Y):
1846 r"""
1847 A trace inner-product for matrices that aren't embedded in the
1848 reals. It takes MATRICES as arguments, not EJA elements.
1849
1850 SETUP::
1851
1852 sage: from mjo.eja.eja_algebra import (RealSymmetricEJA,
1853 ....: ComplexHermitianEJA,
1854 ....: QuaternionHermitianEJA,
1855 ....: OctonionHermitianEJA)
1856
1857 EXAMPLES::
1858
1859 sage: J = RealSymmetricEJA(2,field=QQ,orthonormalize=False)
1860 sage: I = J.one().to_matrix()
1861 sage: J.trace_inner_product(I, -I)
1862 -2
1863
1864 ::
1865
1866 sage: J = ComplexHermitianEJA(2,field=QQ,orthonormalize=False)
1867 sage: I = J.one().to_matrix()
1868 sage: J.trace_inner_product(I, -I)
1869 -2
1870
1871 ::
1872
1873 sage: J = QuaternionHermitianEJA(2,field=QQ,orthonormalize=False)
1874 sage: I = J.one().to_matrix()
1875 sage: J.trace_inner_product(I, -I)
1876 -2
1877
1878 ::
1879
1880 sage: J = OctonionHermitianEJA(2,field=QQ,orthonormalize=False)
1881 sage: I = J.one().to_matrix()
1882 sage: J.trace_inner_product(I, -I)
1883 -2
1884
1885 """
1886 tr = (X*Y).trace()
1887 if hasattr(tr, 'coefficient'):
1888 # Works for octonions, and has to come first because they
1889 # also have a "real()" method that doesn't return an
1890 # element of the scalar ring.
1891 return tr.coefficient(0)
1892 elif hasattr(tr, 'coefficient_tuple'):
1893 # Works for quaternions.
1894 return tr.coefficient_tuple()[0]
1895
1896 # Works for real and complex numbers.
1897 return tr.real()
1898
1899
1900 def __init__(self, matrix_space, **kwargs):
1901 # We know this is a valid EJA, but will double-check
1902 # if the user passes check_axioms=True.
1903 if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
1904
1905
1906 super().__init__(self._denormalized_basis(matrix_space),
1907 self.jordan_product,
1908 self.trace_inner_product,
1909 field=matrix_space.base_ring(),
1910 matrix_space=matrix_space,
1911 **kwargs)
1912
1913 self.rank.set_cache(matrix_space.nrows())
1914 self.one.set_cache( self(matrix_space.one()) )
1915
1916 class RealSymmetricEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA):
1917 """
1918 The rank-n simple EJA consisting of real symmetric n-by-n
1919 matrices, the usual symmetric Jordan product, and the trace inner
1920 product. It has dimension `(n^2 + n)/2` over the reals.
1921
1922 SETUP::
1923
1924 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1925
1926 EXAMPLES::
1927
1928 sage: J = RealSymmetricEJA(2)
1929 sage: b0, b1, b2 = J.gens()
1930 sage: b0*b0
1931 b0
1932 sage: b1*b1
1933 1/2*b0 + 1/2*b2
1934 sage: b2*b2
1935 b2
1936
1937 In theory, our "field" can be any subfield of the reals::
1938
1939 sage: RealSymmetricEJA(2, field=RDF, check_axioms=True)
1940 Euclidean Jordan algebra of dimension 3 over Real Double Field
1941 sage: RealSymmetricEJA(2, field=RR, check_axioms=True)
1942 Euclidean Jordan algebra of dimension 3 over Real Field with
1943 53 bits of precision
1944
1945 TESTS:
1946
1947 The dimension of this algebra is `(n^2 + n) / 2`::
1948
1949 sage: set_random_seed()
1950 sage: d = RealSymmetricEJA._max_random_instance_dimension()
1951 sage: n = RealSymmetricEJA._max_random_instance_size(d)
1952 sage: J = RealSymmetricEJA(n)
1953 sage: J.dimension() == (n^2 + n)/2
1954 True
1955
1956 The Jordan multiplication is what we think it is::
1957
1958 sage: set_random_seed()
1959 sage: J = RealSymmetricEJA.random_instance()
1960 sage: x,y = J.random_elements(2)
1961 sage: actual = (x*y).to_matrix()
1962 sage: X = x.to_matrix()
1963 sage: Y = y.to_matrix()
1964 sage: expected = (X*Y + Y*X)/2
1965 sage: actual == expected
1966 True
1967 sage: J(expected) == x*y
1968 True
1969
1970 We can change the generator prefix::
1971
1972 sage: RealSymmetricEJA(3, prefix='q').gens()
1973 (q0, q1, q2, q3, q4, q5)
1974
1975 We can construct the (trivial) algebra of rank zero::
1976
1977 sage: RealSymmetricEJA(0)
1978 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1979
1980 """
1981 @staticmethod
1982 def _max_random_instance_size(max_dimension):
1983 # Obtained by solving d = (n^2 + n)/2.
1984 # The ZZ-int-ZZ thing is just "floor."
1985 return ZZ(int(ZZ(8*max_dimension + 1).sqrt()/2 - 1/2))
1986
1987 @classmethod
1988 def random_instance(cls, max_dimension=None, **kwargs):
1989 """
1990 Return a random instance of this type of algebra.
1991 """
1992 if max_dimension is None:
1993 max_dimension = cls._max_random_instance_dimension()
1994 max_size = cls._max_random_instance_size(max_dimension) + 1
1995 n = ZZ.random_element(max_size)
1996 return cls(n, **kwargs)
1997
1998 def __init__(self, n, field=AA, **kwargs):
1999 # We know this is a valid EJA, but will double-check
2000 # if the user passes check_axioms=True.
2001 if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
2002
2003 A = MatrixSpace(field, n)
2004 super().__init__(A, **kwargs)
2005
2006 from mjo.eja.eja_cache import real_symmetric_eja_coeffs
2007 a = real_symmetric_eja_coeffs(self)
2008 if a is not None:
2009 if self._rational_algebra is None:
2010 self._charpoly_coefficients.set_cache(a)
2011 else:
2012 self._rational_algebra._charpoly_coefficients.set_cache(a)
2013
2014
2015
2016 class ComplexHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA):
2017 """
2018 The rank-n simple EJA consisting of complex Hermitian n-by-n
2019 matrices over the real numbers, the usual symmetric Jordan product,
2020 and the real-part-of-trace inner product. It has dimension `n^2` over
2021 the reals.
2022
2023 SETUP::
2024
2025 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
2026
2027 EXAMPLES:
2028
2029 In theory, our "field" can be any subfield of the reals, but we
2030 can't use inexact real fields at the moment because SageMath
2031 doesn't know how to convert their elements into complex numbers,
2032 or even into algebraic reals::
2033
2034 sage: QQbar(RDF(1))
2035 Traceback (most recent call last):
2036 ...
2037 TypeError: Illegal initializer for algebraic number
2038 sage: AA(RR(1))
2039 Traceback (most recent call last):
2040 ...
2041 TypeError: Illegal initializer for algebraic number
2042
2043 This causes the following error when we try to scale a matrix of
2044 complex numbers by an inexact real number::
2045
2046 sage: ComplexHermitianEJA(2,field=RR)
2047 Traceback (most recent call last):
2048 ...
2049 TypeError: Unable to coerce entries (=(1.00000000000000,
2050 -0.000000000000000)) to coefficients in Algebraic Real Field
2051
2052 TESTS:
2053
2054 The dimension of this algebra is `n^2`::
2055
2056 sage: set_random_seed()
2057 sage: d = ComplexHermitianEJA._max_random_instance_dimension()
2058 sage: n = ComplexHermitianEJA._max_random_instance_size(d)
2059 sage: J = ComplexHermitianEJA(n)
2060 sage: J.dimension() == n^2
2061 True
2062
2063 The Jordan multiplication is what we think it is::
2064
2065 sage: set_random_seed()
2066 sage: J = ComplexHermitianEJA.random_instance()
2067 sage: x,y = J.random_elements(2)
2068 sage: actual = (x*y).to_matrix()
2069 sage: X = x.to_matrix()
2070 sage: Y = y.to_matrix()
2071 sage: expected = (X*Y + Y*X)/2
2072 sage: actual == expected
2073 True
2074 sage: J(expected) == x*y
2075 True
2076
2077 We can change the generator prefix::
2078
2079 sage: ComplexHermitianEJA(2, prefix='z').gens()
2080 (z0, z1, z2, z3)
2081
2082 We can construct the (trivial) algebra of rank zero::
2083
2084 sage: ComplexHermitianEJA(0)
2085 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2086
2087 """
2088 def __init__(self, n, field=AA, **kwargs):
2089 # We know this is a valid EJA, but will double-check
2090 # if the user passes check_axioms=True.
2091 if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
2092
2093 from mjo.hurwitz import ComplexMatrixAlgebra
2094 A = ComplexMatrixAlgebra(n, scalars=field)
2095 super().__init__(A, **kwargs)
2096
2097 from mjo.eja.eja_cache import complex_hermitian_eja_coeffs
2098 a = complex_hermitian_eja_coeffs(self)
2099 if a is not None:
2100 if self._rational_algebra is None:
2101 self._charpoly_coefficients.set_cache(a)
2102 else:
2103 self._rational_algebra._charpoly_coefficients.set_cache(a)
2104
2105 @staticmethod
2106 def _max_random_instance_size(max_dimension):
2107 # Obtained by solving d = n^2.
2108 # The ZZ-int-ZZ thing is just "floor."
2109 return ZZ(int(ZZ(max_dimension).sqrt()))
2110
2111 @classmethod
2112 def random_instance(cls, max_dimension=None, **kwargs):
2113 """
2114 Return a random instance of this type of algebra.
2115 """
2116 if max_dimension is None:
2117 max_dimension = cls._max_random_instance_dimension()
2118 n = ZZ.random_element(cls._max_random_instance_size(max_dimension) + 1)
2119 return cls(n, **kwargs)
2120
2121
2122 class QuaternionHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA):
2123 r"""
2124 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
2125 matrices, the usual symmetric Jordan product, and the
2126 real-part-of-trace inner product. It has dimension `2n^2 - n` over
2127 the reals.
2128
2129 SETUP::
2130
2131 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2132
2133 EXAMPLES:
2134
2135 In theory, our "field" can be any subfield of the reals::
2136
2137 sage: QuaternionHermitianEJA(2, field=RDF, check_axioms=True)
2138 Euclidean Jordan algebra of dimension 6 over Real Double Field
2139 sage: QuaternionHermitianEJA(2, field=RR, check_axioms=True)
2140 Euclidean Jordan algebra of dimension 6 over Real Field with
2141 53 bits of precision
2142
2143 TESTS:
2144
2145 The dimension of this algebra is `2*n^2 - n`::
2146
2147 sage: set_random_seed()
2148 sage: d = QuaternionHermitianEJA._max_random_instance_dimension()
2149 sage: n = QuaternionHermitianEJA._max_random_instance_size(d)
2150 sage: J = QuaternionHermitianEJA(n)
2151 sage: J.dimension() == 2*(n^2) - n
2152 True
2153
2154 The Jordan multiplication is what we think it is::
2155
2156 sage: set_random_seed()
2157 sage: J = QuaternionHermitianEJA.random_instance()
2158 sage: x,y = J.random_elements(2)
2159 sage: actual = (x*y).to_matrix()
2160 sage: X = x.to_matrix()
2161 sage: Y = y.to_matrix()
2162 sage: expected = (X*Y + Y*X)/2
2163 sage: actual == expected
2164 True
2165 sage: J(expected) == x*y
2166 True
2167
2168 We can change the generator prefix::
2169
2170 sage: QuaternionHermitianEJA(2, prefix='a').gens()
2171 (a0, a1, a2, a3, a4, a5)
2172
2173 We can construct the (trivial) algebra of rank zero::
2174
2175 sage: QuaternionHermitianEJA(0)
2176 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2177
2178 """
2179 def __init__(self, n, field=AA, **kwargs):
2180 # We know this is a valid EJA, but will double-check
2181 # if the user passes check_axioms=True.
2182 if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
2183
2184 from mjo.hurwitz import QuaternionMatrixAlgebra
2185 A = QuaternionMatrixAlgebra(n, scalars=field)
2186 super().__init__(A, **kwargs)
2187
2188 from mjo.eja.eja_cache import quaternion_hermitian_eja_coeffs
2189 a = quaternion_hermitian_eja_coeffs(self)
2190 if a is not None:
2191 if self._rational_algebra is None:
2192 self._charpoly_coefficients.set_cache(a)
2193 else:
2194 self._rational_algebra._charpoly_coefficients.set_cache(a)
2195
2196
2197
2198 @staticmethod
2199 def _max_random_instance_size(max_dimension):
2200 r"""
2201 The maximum rank of a random QuaternionHermitianEJA.
2202 """
2203 # Obtained by solving d = 2n^2 - n.
2204 # The ZZ-int-ZZ thing is just "floor."
2205 return ZZ(int(ZZ(8*max_dimension + 1).sqrt()/4 + 1/4))
2206
2207 @classmethod
2208 def random_instance(cls, max_dimension=None, **kwargs):
2209 """
2210 Return a random instance of this type of algebra.
2211 """
2212 if max_dimension is None:
2213 max_dimension = cls._max_random_instance_dimension()
2214 n = ZZ.random_element(cls._max_random_instance_size(max_dimension) + 1)
2215 return cls(n, **kwargs)
2216
2217 class OctonionHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA):
2218 r"""
2219 SETUP::
2220
2221 sage: from mjo.eja.eja_algebra import (FiniteDimensionalEJA,
2222 ....: OctonionHermitianEJA)
2223 sage: from mjo.hurwitz import Octonions, OctonionMatrixAlgebra
2224
2225 EXAMPLES:
2226
2227 The 3-by-3 algebra satisfies the axioms of an EJA::
2228
2229 sage: OctonionHermitianEJA(3, # long time
2230 ....: field=QQ, # long time
2231 ....: orthonormalize=False, # long time
2232 ....: check_axioms=True) # long time
2233 Euclidean Jordan algebra of dimension 27 over Rational Field
2234
2235 After a change-of-basis, the 2-by-2 algebra has the same
2236 multiplication table as the ten-dimensional Jordan spin algebra::
2237
2238 sage: A = OctonionMatrixAlgebra(2,Octonions(QQ),QQ)
2239 sage: b = OctonionHermitianEJA._denormalized_basis(A)
2240 sage: basis = (b[0] + b[9],) + b[1:9] + (b[0] - b[9],)
2241 sage: jp = OctonionHermitianEJA.jordan_product
2242 sage: ip = OctonionHermitianEJA.trace_inner_product
2243 sage: J = FiniteDimensionalEJA(basis,
2244 ....: jp,
2245 ....: ip,
2246 ....: field=QQ,
2247 ....: orthonormalize=False)
2248 sage: J.multiplication_table()
2249 +----++----+----+----+----+----+----+----+----+----+----+
2250 | * || b0 | b1 | b2 | b3 | b4 | b5 | b6 | b7 | b8 | b9 |
2251 +====++====+====+====+====+====+====+====+====+====+====+
2252 | b0 || b0 | b1 | b2 | b3 | b4 | b5 | b6 | b7 | b8 | b9 |
2253 +----++----+----+----+----+----+----+----+----+----+----+
2254 | b1 || b1 | b0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
2255 +----++----+----+----+----+----+----+----+----+----+----+
2256 | b2 || b2 | 0 | b0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
2257 +----++----+----+----+----+----+----+----+----+----+----+
2258 | b3 || b3 | 0 | 0 | b0 | 0 | 0 | 0 | 0 | 0 | 0 |
2259 +----++----+----+----+----+----+----+----+----+----+----+
2260 | b4 || b4 | 0 | 0 | 0 | b0 | 0 | 0 | 0 | 0 | 0 |
2261 +----++----+----+----+----+----+----+----+----+----+----+
2262 | b5 || b5 | 0 | 0 | 0 | 0 | b0 | 0 | 0 | 0 | 0 |
2263 +----++----+----+----+----+----+----+----+----+----+----+
2264 | b6 || b6 | 0 | 0 | 0 | 0 | 0 | b0 | 0 | 0 | 0 |
2265 +----++----+----+----+----+----+----+----+----+----+----+
2266 | b7 || b7 | 0 | 0 | 0 | 0 | 0 | 0 | b0 | 0 | 0 |
2267 +----++----+----+----+----+----+----+----+----+----+----+
2268 | b8 || b8 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | b0 | 0 |
2269 +----++----+----+----+----+----+----+----+----+----+----+
2270 | b9 || b9 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | b0 |
2271 +----++----+----+----+----+----+----+----+----+----+----+
2272
2273 TESTS:
2274
2275 We can actually construct the 27-dimensional Albert algebra,
2276 and we get the right unit element if we recompute it::
2277
2278 sage: J = OctonionHermitianEJA(3, # long time
2279 ....: field=QQ, # long time
2280 ....: orthonormalize=False) # long time
2281 sage: J.one.clear_cache() # long time
2282 sage: J.one() # long time
2283 b0 + b9 + b26
2284 sage: J.one().to_matrix() # long time
2285 +----+----+----+
2286 | e0 | 0 | 0 |
2287 +----+----+----+
2288 | 0 | e0 | 0 |
2289 +----+----+----+
2290 | 0 | 0 | e0 |
2291 +----+----+----+
2292
2293 The 2-by-2 algebra is isomorphic to the ten-dimensional Jordan
2294 spin algebra, but just to be sure, we recompute its rank::
2295
2296 sage: J = OctonionHermitianEJA(2, # long time
2297 ....: field=QQ, # long time
2298 ....: orthonormalize=False) # long time
2299 sage: J.rank.clear_cache() # long time
2300 sage: J.rank() # long time
2301 2
2302
2303 """
2304 @staticmethod
2305 def _max_random_instance_size(max_dimension):
2306 r"""
2307 The maximum rank of a random QuaternionHermitianEJA.
2308 """
2309 # There's certainly a formula for this, but with only four
2310 # cases to worry about, I'm not that motivated to derive it.
2311 if max_dimension >= 27:
2312 return 3
2313 elif max_dimension >= 10:
2314 return 2
2315 elif max_dimension >= 1:
2316 return 1
2317 else:
2318 return 0
2319
2320 @classmethod
2321 def random_instance(cls, max_dimension=None, **kwargs):
2322 """
2323 Return a random instance of this type of algebra.
2324 """
2325 if max_dimension is None:
2326 max_dimension = cls._max_random_instance_dimension()
2327 n = ZZ.random_element(cls._max_random_instance_size(max_dimension) + 1)
2328 return cls(n, **kwargs)
2329
2330 def __init__(self, n, field=AA, **kwargs):
2331 if n > 3:
2332 # Otherwise we don't get an EJA.
2333 raise ValueError("n cannot exceed 3")
2334
2335 # We know this is a valid EJA, but will double-check
2336 # if the user passes check_axioms=True.
2337 if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
2338
2339 from mjo.hurwitz import OctonionMatrixAlgebra
2340 A = OctonionMatrixAlgebra(n, scalars=field)
2341 super().__init__(A, **kwargs)
2342
2343 from mjo.eja.eja_cache import octonion_hermitian_eja_coeffs
2344 a = octonion_hermitian_eja_coeffs(self)
2345 if a is not None:
2346 if self._rational_algebra is None:
2347 self._charpoly_coefficients.set_cache(a)
2348 else:
2349 self._rational_algebra._charpoly_coefficients.set_cache(a)
2350
2351
2352 class AlbertEJA(OctonionHermitianEJA):
2353 r"""
2354 The Albert algebra is the algebra of three-by-three Hermitian
2355 matrices whose entries are octonions.
2356
2357 SETUP::
2358
2359 sage: from mjo.eja.eja_algebra import AlbertEJA
2360
2361 EXAMPLES::
2362
2363 sage: AlbertEJA(field=QQ, orthonormalize=False)
2364 Euclidean Jordan algebra of dimension 27 over Rational Field
2365 sage: AlbertEJA() # long time
2366 Euclidean Jordan algebra of dimension 27 over Algebraic Real Field
2367
2368 """
2369 def __init__(self, *args, **kwargs):
2370 super().__init__(3, *args, **kwargs)
2371
2372
2373 class HadamardEJA(RationalBasisEJA, ConcreteEJA):
2374 """
2375 Return the Euclidean Jordan algebra on `R^n` with the Hadamard
2376 (pointwise real-number multiplication) Jordan product and the
2377 usual inner-product.
2378
2379 This is nothing more than the Cartesian product of ``n`` copies of
2380 the one-dimensional Jordan spin algebra, and is the most common
2381 example of a non-simple Euclidean Jordan algebra.
2382
2383 SETUP::
2384
2385 sage: from mjo.eja.eja_algebra import HadamardEJA
2386
2387 EXAMPLES:
2388
2389 This multiplication table can be verified by hand::
2390
2391 sage: J = HadamardEJA(3)
2392 sage: b0,b1,b2 = J.gens()
2393 sage: b0*b0
2394 b0
2395 sage: b0*b1
2396 0
2397 sage: b0*b2
2398 0
2399 sage: b1*b1
2400 b1
2401 sage: b1*b2
2402 0
2403 sage: b2*b2
2404 b2
2405
2406 TESTS:
2407
2408 We can change the generator prefix::
2409
2410 sage: HadamardEJA(3, prefix='r').gens()
2411 (r0, r1, r2)
2412 """
2413 def __init__(self, n, field=AA, **kwargs):
2414 MS = MatrixSpace(field, n, 1)
2415
2416 if n == 0:
2417 jordan_product = lambda x,y: x
2418 inner_product = lambda x,y: x
2419 else:
2420 def jordan_product(x,y):
2421 return MS( xi*yi for (xi,yi) in zip(x,y) )
2422
2423 def inner_product(x,y):
2424 return (x.T*y)[0,0]
2425
2426 # New defaults for keyword arguments. Don't orthonormalize
2427 # because our basis is already orthonormal with respect to our
2428 # inner-product. Don't check the axioms, because we know this
2429 # is a valid EJA... but do double-check if the user passes
2430 # check_axioms=True. Note: we DON'T override the "check_field"
2431 # default here, because the user can pass in a field!
2432 if "orthonormalize" not in kwargs: kwargs["orthonormalize"] = False
2433 if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
2434
2435 column_basis = tuple( MS(b) for b in FreeModule(field, n).basis() )
2436 super().__init__(column_basis,
2437 jordan_product,
2438 inner_product,
2439 field=field,
2440 matrix_space=MS,
2441 associative=True,
2442 **kwargs)
2443 self.rank.set_cache(n)
2444
2445 self.one.set_cache( self.sum(self.gens()) )
2446
2447 @staticmethod
2448 def _max_random_instance_dimension():
2449 r"""
2450 There's no reason to go higher than five here. That's
2451 enough to get the point across.
2452 """
2453 return 5
2454
2455 @staticmethod
2456 def _max_random_instance_size(max_dimension):
2457 r"""
2458 The maximum size (=dimension) of a random HadamardEJA.
2459 """
2460 return max_dimension
2461
2462 @classmethod
2463 def random_instance(cls, max_dimension=None, **kwargs):
2464 """
2465 Return a random instance of this type of algebra.
2466 """
2467 if max_dimension is None:
2468 max_dimension = cls._max_random_instance_dimension()
2469 n = ZZ.random_element(cls._max_random_instance_size(max_dimension) + 1)
2470 return cls(n, **kwargs)
2471
2472
2473 class BilinearFormEJA(RationalBasisEJA, ConcreteEJA):
2474 r"""
2475 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2476 with the half-trace inner product and jordan product ``x*y =
2477 (<Bx,y>,y_bar>, x0*y_bar + y0*x_bar)`` where `B = 1 \times B22` is
2478 a symmetric positive-definite "bilinear form" matrix. Its
2479 dimension is the size of `B`, and it has rank two in dimensions
2480 larger than two. It reduces to the ``JordanSpinEJA`` when `B` is
2481 the identity matrix of order ``n``.
2482
2483 We insist that the one-by-one upper-left identity block of `B` be
2484 passed in as well so that we can be passed a matrix of size zero
2485 to construct a trivial algebra.
2486
2487 SETUP::
2488
2489 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
2490 ....: JordanSpinEJA)
2491
2492 EXAMPLES:
2493
2494 When no bilinear form is specified, the identity matrix is used,
2495 and the resulting algebra is the Jordan spin algebra::
2496
2497 sage: B = matrix.identity(AA,3)
2498 sage: J0 = BilinearFormEJA(B)
2499 sage: J1 = JordanSpinEJA(3)
2500 sage: J0.multiplication_table() == J0.multiplication_table()
2501 True
2502
2503 An error is raised if the matrix `B` does not correspond to a
2504 positive-definite bilinear form::
2505
2506 sage: B = matrix.random(QQ,2,3)
2507 sage: J = BilinearFormEJA(B)
2508 Traceback (most recent call last):
2509 ...
2510 ValueError: bilinear form is not positive-definite
2511 sage: B = matrix.zero(QQ,3)
2512 sage: J = BilinearFormEJA(B)
2513 Traceback (most recent call last):
2514 ...
2515 ValueError: bilinear form is not positive-definite
2516
2517 TESTS:
2518
2519 We can create a zero-dimensional algebra::
2520
2521 sage: B = matrix.identity(AA,0)
2522 sage: J = BilinearFormEJA(B)
2523 sage: J.basis()
2524 Finite family {}
2525
2526 We can check the multiplication condition given in the Jordan, von
2527 Neumann, and Wigner paper (and also discussed on my "On the
2528 symmetry..." paper). Note that this relies heavily on the standard
2529 choice of basis, as does anything utilizing the bilinear form
2530 matrix. We opt not to orthonormalize the basis, because if we
2531 did, we would have to normalize the `s_{i}` in a similar manner::
2532
2533 sage: set_random_seed()
2534 sage: n = ZZ.random_element(5)
2535 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2536 sage: B11 = matrix.identity(QQ,1)
2537 sage: B22 = M.transpose()*M
2538 sage: B = block_matrix(2,2,[ [B11,0 ],
2539 ....: [0, B22 ] ])
2540 sage: J = BilinearFormEJA(B, orthonormalize=False)
2541 sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
2542 sage: V = J.vector_space()
2543 sage: sis = [ J( V([0] + (M.inverse()*ei).list()).column() )
2544 ....: for ei in eis ]
2545 sage: actual = [ sis[i]*sis[j]
2546 ....: for i in range(n-1)
2547 ....: for j in range(n-1) ]
2548 sage: expected = [ J.one() if i == j else J.zero()
2549 ....: for i in range(n-1)
2550 ....: for j in range(n-1) ]
2551 sage: actual == expected
2552 True
2553
2554 """
2555 def __init__(self, B, field=AA, **kwargs):
2556 # The matrix "B" is supplied by the user in most cases,
2557 # so it makes sense to check whether or not its positive-
2558 # definite unless we are specifically asked not to...
2559 if ("check_axioms" not in kwargs) or kwargs["check_axioms"]:
2560 if not B.is_positive_definite():
2561 raise ValueError("bilinear form is not positive-definite")
2562
2563 # However, all of the other data for this EJA is computed
2564 # by us in manner that guarantees the axioms are
2565 # satisfied. So, again, unless we are specifically asked to
2566 # verify things, we'll skip the rest of the checks.
2567 if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
2568
2569 n = B.nrows()
2570 MS = MatrixSpace(field, n, 1)
2571
2572 def inner_product(x,y):
2573 return (y.T*B*x)[0,0]
2574
2575 def jordan_product(x,y):
2576 x0 = x[0,0]
2577 xbar = x[1:,0]
2578 y0 = y[0,0]
2579 ybar = y[1:,0]
2580 z0 = inner_product(y,x)
2581 zbar = y0*xbar + x0*ybar
2582 return MS([z0] + zbar.list())
2583
2584 column_basis = tuple( MS(b) for b in FreeModule(field, n).basis() )
2585
2586 # TODO: I haven't actually checked this, but it seems legit.
2587 associative = False
2588 if n <= 2:
2589 associative = True
2590
2591 super().__init__(column_basis,
2592 jordan_product,
2593 inner_product,
2594 field=field,
2595 matrix_space=MS,
2596 associative=associative,
2597 **kwargs)
2598
2599 # The rank of this algebra is two, unless we're in a
2600 # one-dimensional ambient space (because the rank is bounded
2601 # by the ambient dimension).
2602 self.rank.set_cache(min(n,2))
2603 if n == 0:
2604 self.one.set_cache( self.zero() )
2605 else:
2606 self.one.set_cache( self.monomial(0) )
2607
2608 @staticmethod
2609 def _max_random_instance_dimension():
2610 r"""
2611 There's no reason to go higher than five here. That's
2612 enough to get the point across.
2613 """
2614 return 5
2615
2616 @staticmethod
2617 def _max_random_instance_size(max_dimension):
2618 r"""
2619 The maximum size (=dimension) of a random BilinearFormEJA.
2620 """
2621 return max_dimension
2622
2623 @classmethod
2624 def random_instance(cls, max_dimension=None, **kwargs):
2625 """
2626 Return a random instance of this algebra.
2627 """
2628 if max_dimension is None:
2629 max_dimension = cls._max_random_instance_dimension()
2630 n = ZZ.random_element(cls._max_random_instance_size(max_dimension) + 1)
2631 if n.is_zero():
2632 B = matrix.identity(ZZ, n)
2633 return cls(B, **kwargs)
2634
2635 B11 = matrix.identity(ZZ, 1)
2636 M = matrix.random(ZZ, n-1)
2637 I = matrix.identity(ZZ, n-1)
2638 alpha = ZZ.zero()
2639 while alpha.is_zero():
2640 alpha = ZZ.random_element().abs()
2641 B22 = M.transpose()*M + alpha*I
2642
2643 from sage.matrix.special import block_matrix
2644 B = block_matrix(2,2, [ [B11, ZZ(0) ],
2645 [ZZ(0), B22 ] ])
2646
2647 return cls(B, **kwargs)
2648
2649
2650 class JordanSpinEJA(BilinearFormEJA):
2651 """
2652 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2653 with the usual inner product and jordan product ``x*y =
2654 (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2655 the reals.
2656
2657 SETUP::
2658
2659 sage: from mjo.eja.eja_algebra import JordanSpinEJA
2660
2661 EXAMPLES:
2662
2663 This multiplication table can be verified by hand::
2664
2665 sage: J = JordanSpinEJA(4)
2666 sage: b0,b1,b2,b3 = J.gens()
2667 sage: b0*b0
2668 b0
2669 sage: b0*b1
2670 b1
2671 sage: b0*b2
2672 b2
2673 sage: b0*b3
2674 b3
2675 sage: b1*b2
2676 0
2677 sage: b1*b3
2678 0
2679 sage: b2*b3
2680 0
2681
2682 We can change the generator prefix::
2683
2684 sage: JordanSpinEJA(2, prefix='B').gens()
2685 (B0, B1)
2686
2687 TESTS:
2688
2689 Ensure that we have the usual inner product on `R^n`::
2690
2691 sage: set_random_seed()
2692 sage: J = JordanSpinEJA.random_instance()
2693 sage: x,y = J.random_elements(2)
2694 sage: actual = x.inner_product(y)
2695 sage: expected = x.to_vector().inner_product(y.to_vector())
2696 sage: actual == expected
2697 True
2698
2699 """
2700 def __init__(self, n, *args, **kwargs):
2701 # This is a special case of the BilinearFormEJA with the
2702 # identity matrix as its bilinear form.
2703 B = matrix.identity(ZZ, n)
2704
2705 # Don't orthonormalize because our basis is already
2706 # orthonormal with respect to our inner-product.
2707 if "orthonormalize" not in kwargs: kwargs["orthonormalize"] = False
2708
2709 # But also don't pass check_field=False here, because the user
2710 # can pass in a field!
2711 super().__init__(B, *args, **kwargs)
2712
2713 @classmethod
2714 def random_instance(cls, max_dimension=None, **kwargs):
2715 """
2716 Return a random instance of this type of algebra.
2717
2718 Needed here to override the implementation for ``BilinearFormEJA``.
2719 """
2720 if max_dimension is None:
2721 max_dimension = cls._max_random_instance_dimension()
2722 n = ZZ.random_element(cls._max_random_instance_size(max_dimension) + 1)
2723 return cls(n, **kwargs)
2724
2725
2726 class TrivialEJA(RationalBasisEJA, ConcreteEJA):
2727 """
2728 The trivial Euclidean Jordan algebra consisting of only a zero element.
2729
2730 SETUP::
2731
2732 sage: from mjo.eja.eja_algebra import TrivialEJA
2733
2734 EXAMPLES::
2735
2736 sage: J = TrivialEJA()
2737 sage: J.dimension()
2738 0
2739 sage: J.zero()
2740 0
2741 sage: J.one()
2742 0
2743 sage: 7*J.one()*12*J.one()
2744 0
2745 sage: J.one().inner_product(J.one())
2746 0
2747 sage: J.one().norm()
2748 0
2749 sage: J.one().subalgebra_generated_by()
2750 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2751 sage: J.rank()
2752 0
2753
2754 """
2755 def __init__(self, field=AA, **kwargs):
2756 jordan_product = lambda x,y: x
2757 inner_product = lambda x,y: field.zero()
2758 basis = ()
2759 MS = MatrixSpace(field,0)
2760
2761 # New defaults for keyword arguments
2762 if "orthonormalize" not in kwargs: kwargs["orthonormalize"] = False
2763 if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
2764
2765 super().__init__(basis,
2766 jordan_product,
2767 inner_product,
2768 associative=True,
2769 field=field,
2770 matrix_space=MS,
2771 **kwargs)
2772
2773 # The rank is zero using my definition, namely the dimension of the
2774 # largest subalgebra generated by any element.
2775 self.rank.set_cache(0)
2776 self.one.set_cache( self.zero() )
2777
2778 @classmethod
2779 def random_instance(cls, **kwargs):
2780 # We don't take a "size" argument so the superclass method is
2781 # inappropriate for us.
2782 return cls(**kwargs)
2783
2784
2785 class CartesianProductEJA(FiniteDimensionalEJA):
2786 r"""
2787 The external (orthogonal) direct sum of two or more Euclidean
2788 Jordan algebras. Every Euclidean Jordan algebra decomposes into an
2789 orthogonal direct sum of simple Euclidean Jordan algebras which is
2790 then isometric to a Cartesian product, so no generality is lost by
2791 providing only this construction.
2792
2793 SETUP::
2794
2795 sage: from mjo.eja.eja_algebra import (random_eja,
2796 ....: CartesianProductEJA,
2797 ....: HadamardEJA,
2798 ....: JordanSpinEJA,
2799 ....: RealSymmetricEJA)
2800
2801 EXAMPLES:
2802
2803 The Jordan product is inherited from our factors and implemented by
2804 our CombinatorialFreeModule Cartesian product superclass::
2805
2806 sage: set_random_seed()
2807 sage: J1 = HadamardEJA(2)
2808 sage: J2 = RealSymmetricEJA(2)
2809 sage: J = cartesian_product([J1,J2])
2810 sage: x,y = J.random_elements(2)
2811 sage: x*y in J
2812 True
2813
2814 The ability to retrieve the original factors is implemented by our
2815 CombinatorialFreeModule Cartesian product superclass::
2816
2817 sage: J1 = HadamardEJA(2, field=QQ)
2818 sage: J2 = JordanSpinEJA(3, field=QQ)
2819 sage: J = cartesian_product([J1,J2])
2820 sage: J.cartesian_factors()
2821 (Euclidean Jordan algebra of dimension 2 over Rational Field,
2822 Euclidean Jordan algebra of dimension 3 over Rational Field)
2823
2824 You can provide more than two factors::
2825
2826 sage: J1 = HadamardEJA(2)
2827 sage: J2 = JordanSpinEJA(3)
2828 sage: J3 = RealSymmetricEJA(3)
2829 sage: cartesian_product([J1,J2,J3])
2830 Euclidean Jordan algebra of dimension 2 over Algebraic Real
2831 Field (+) Euclidean Jordan algebra of dimension 3 over Algebraic
2832 Real Field (+) Euclidean Jordan algebra of dimension 6 over
2833 Algebraic Real Field
2834
2835 Rank is additive on a Cartesian product::
2836
2837 sage: J1 = HadamardEJA(1)
2838 sage: J2 = RealSymmetricEJA(2)
2839 sage: J = cartesian_product([J1,J2])
2840 sage: J1.rank.clear_cache()
2841 sage: J2.rank.clear_cache()
2842 sage: J.rank.clear_cache()
2843 sage: J.rank()
2844 3
2845 sage: J.rank() == J1.rank() + J2.rank()
2846 True
2847
2848 The same rank computation works over the rationals, with whatever
2849 basis you like::
2850
2851 sage: J1 = HadamardEJA(1, field=QQ, orthonormalize=False)
2852 sage: J2 = RealSymmetricEJA(2, field=QQ, orthonormalize=False)
2853 sage: J = cartesian_product([J1,J2])
2854 sage: J1.rank.clear_cache()
2855 sage: J2.rank.clear_cache()
2856 sage: J.rank.clear_cache()
2857 sage: J.rank()
2858 3
2859 sage: J.rank() == J1.rank() + J2.rank()
2860 True
2861
2862 The product algebra will be associative if and only if all of its
2863 components are associative::
2864
2865 sage: J1 = HadamardEJA(2)
2866 sage: J1.is_associative()
2867 True
2868 sage: J2 = HadamardEJA(3)
2869 sage: J2.is_associative()
2870 True
2871 sage: J3 = RealSymmetricEJA(3)
2872 sage: J3.is_associative()
2873 False
2874 sage: CP1 = cartesian_product([J1,J2])
2875 sage: CP1.is_associative()
2876 True
2877 sage: CP2 = cartesian_product([J1,J3])
2878 sage: CP2.is_associative()
2879 False
2880
2881 Cartesian products of Cartesian products work::
2882
2883 sage: J1 = JordanSpinEJA(1)
2884 sage: J2 = JordanSpinEJA(1)
2885 sage: J3 = JordanSpinEJA(1)
2886 sage: J = cartesian_product([J1,cartesian_product([J2,J3])])
2887 sage: J.multiplication_table()
2888 +----++----+----+----+
2889 | * || b0 | b1 | b2 |
2890 +====++====+====+====+
2891 | b0 || b0 | 0 | 0 |
2892 +----++----+----+----+
2893 | b1 || 0 | b1 | 0 |
2894 +----++----+----+----+
2895 | b2 || 0 | 0 | b2 |
2896 +----++----+----+----+
2897 sage: HadamardEJA(3).multiplication_table()
2898 +----++----+----+----+
2899 | * || b0 | b1 | b2 |
2900 +====++====+====+====+
2901 | b0 || b0 | 0 | 0 |
2902 +----++----+----+----+
2903 | b1 || 0 | b1 | 0 |
2904 +----++----+----+----+
2905 | b2 || 0 | 0 | b2 |
2906 +----++----+----+----+
2907
2908 TESTS:
2909
2910 All factors must share the same base field::
2911
2912 sage: J1 = HadamardEJA(2, field=QQ)
2913 sage: J2 = RealSymmetricEJA(2)
2914 sage: CartesianProductEJA((J1,J2))
2915 Traceback (most recent call last):
2916 ...
2917 ValueError: all factors must share the same base field
2918
2919 The cached unit element is the same one that would be computed::
2920
2921 sage: set_random_seed() # long time
2922 sage: J1 = random_eja() # long time
2923 sage: J2 = random_eja() # long time
2924 sage: J = cartesian_product([J1,J2]) # long time
2925 sage: actual = J.one() # long time
2926 sage: J.one.clear_cache() # long time
2927 sage: expected = J.one() # long time
2928 sage: actual == expected # long time
2929 True
2930
2931 """
2932 Element = FiniteDimensionalEJAElement
2933
2934
2935 def __init__(self, factors, **kwargs):
2936 m = len(factors)
2937 if m == 0:
2938 return TrivialEJA()
2939
2940 self._sets = factors
2941
2942 field = factors[0].base_ring()
2943 if not all( J.base_ring() == field for J in factors ):
2944 raise ValueError("all factors must share the same base field")
2945
2946 associative = all( f.is_associative() for f in factors )
2947
2948 # Compute my matrix space. This category isn't perfect, but
2949 # is good enough for what we need to do.
2950 MS_cat = MagmaticAlgebras(field).FiniteDimensional().WithBasis()
2951 MS_cat = MS_cat.Unital().CartesianProducts()
2952 MS_factors = tuple( J.matrix_space() for J in factors )
2953 from sage.sets.cartesian_product import CartesianProduct
2954 MS = CartesianProduct(MS_factors, MS_cat)
2955
2956 basis = []
2957 zero = MS.zero()
2958 for i in range(m):
2959 for b in factors[i].matrix_basis():
2960 z = list(zero)
2961 z[i] = b
2962 basis.append(z)
2963
2964 basis = tuple( MS(b) for b in basis )
2965
2966 # Define jordan/inner products that operate on that matrix_basis.
2967 def jordan_product(x,y):
2968 return MS(tuple(
2969 (factors[i](x[i])*factors[i](y[i])).to_matrix()
2970 for i in range(m)
2971 ))
2972
2973 def inner_product(x, y):
2974 return sum(
2975 factors[i](x[i]).inner_product(factors[i](y[i]))
2976 for i in range(m)
2977 )
2978
2979 # There's no need to check the field since it already came
2980 # from an EJA. Likewise the axioms are guaranteed to be
2981 # satisfied, unless the guy writing this class sucks.
2982 #
2983 # If you want the basis to be orthonormalized, orthonormalize
2984 # the factors.
2985 FiniteDimensionalEJA.__init__(self,
2986 basis,
2987 jordan_product,
2988 inner_product,
2989 field=field,
2990 matrix_space=MS,
2991 orthonormalize=False,
2992 associative=associative,
2993 cartesian_product=True,
2994 check_field=False,
2995 check_axioms=False)
2996
2997 self.rank.set_cache(sum(J.rank() for J in factors))
2998 ones = tuple(J.one().to_matrix() for J in factors)
2999 self.one.set_cache(self(ones))
3000
3001 def cartesian_factors(self):
3002 # Copy/pasted from CombinatorialFreeModule_CartesianProduct.
3003 return self._sets
3004
3005 def cartesian_factor(self, i):
3006 r"""
3007 Return the ``i``th factor of this algebra.
3008 """
3009 return self._sets[i]
3010
3011 def _repr_(self):
3012 # Copy/pasted from CombinatorialFreeModule_CartesianProduct.
3013 from sage.categories.cartesian_product import cartesian_product
3014 return cartesian_product.symbol.join("%s" % factor
3015 for factor in self._sets)
3016
3017 def matrix_space(self):
3018 r"""
3019 Return the space that our matrix basis lives in as a Cartesian
3020 product.
3021
3022 We don't simply use the ``cartesian_product()`` functor here
3023 because it acts differently on SageMath MatrixSpaces and our
3024 custom MatrixAlgebras, which are CombinatorialFreeModules. We
3025 always want the result to be represented (and indexed) as
3026 an ordered tuple.
3027
3028 SETUP::
3029
3030 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
3031 ....: HadamardEJA,
3032 ....: OctonionHermitianEJA,
3033 ....: RealSymmetricEJA)
3034
3035 EXAMPLES::
3036
3037 sage: J1 = HadamardEJA(1)
3038 sage: J2 = RealSymmetricEJA(2)
3039 sage: J = cartesian_product([J1,J2])
3040 sage: J.matrix_space()
3041 The Cartesian product of (Full MatrixSpace of 1 by 1 dense
3042 matrices over Algebraic Real Field, Full MatrixSpace of 2
3043 by 2 dense matrices over Algebraic Real Field)
3044
3045 ::
3046
3047 sage: J1 = ComplexHermitianEJA(1)
3048 sage: J2 = ComplexHermitianEJA(1)
3049 sage: J = cartesian_product([J1,J2])
3050 sage: J.one().to_matrix()[0]
3051 +---+
3052 | 1 |
3053 +---+
3054 sage: J.one().to_matrix()[1]
3055 +---+
3056 | 1 |
3057 +---+
3058
3059 ::
3060
3061 sage: J1 = OctonionHermitianEJA(1)
3062 sage: J2 = OctonionHermitianEJA(1)
3063 sage: J = cartesian_product([J1,J2])
3064 sage: J.one().to_matrix()[0]
3065 +----+
3066 | e0 |
3067 +----+
3068 sage: J.one().to_matrix()[1]
3069 +----+
3070 | e0 |
3071 +----+
3072
3073 """
3074 return super().matrix_space()
3075
3076
3077 @cached_method
3078 def cartesian_projection(self, i):
3079 r"""
3080 SETUP::
3081
3082 sage: from mjo.eja.eja_algebra import (random_eja,
3083 ....: JordanSpinEJA,
3084 ....: HadamardEJA,
3085 ....: RealSymmetricEJA,
3086 ....: ComplexHermitianEJA)
3087
3088 EXAMPLES:
3089
3090 The projection morphisms are Euclidean Jordan algebra
3091 operators::
3092
3093 sage: J1 = HadamardEJA(2)
3094 sage: J2 = RealSymmetricEJA(2)
3095 sage: J = cartesian_product([J1,J2])
3096 sage: J.cartesian_projection(0)
3097 Linear operator between finite-dimensional Euclidean Jordan
3098 algebras represented by the matrix:
3099 [1 0 0 0 0]
3100 [0 1 0 0 0]
3101 Domain: Euclidean Jordan algebra of dimension 2 over Algebraic
3102 Real Field (+) Euclidean Jordan algebra of dimension 3 over
3103 Algebraic Real Field
3104 Codomain: Euclidean Jordan algebra of dimension 2 over Algebraic
3105 Real Field
3106 sage: J.cartesian_projection(1)
3107 Linear operator between finite-dimensional Euclidean Jordan
3108 algebras represented by the matrix:
3109 [0 0 1 0 0]
3110 [0 0 0 1 0]
3111 [0 0 0 0 1]
3112 Domain: Euclidean Jordan algebra of dimension 2 over Algebraic
3113 Real Field (+) Euclidean Jordan algebra of dimension 3 over
3114 Algebraic Real Field
3115 Codomain: Euclidean Jordan algebra of dimension 3 over Algebraic
3116 Real Field
3117
3118 The projections work the way you'd expect on the vector
3119 representation of an element::
3120
3121 sage: J1 = JordanSpinEJA(2)
3122 sage: J2 = ComplexHermitianEJA(2)
3123 sage: J = cartesian_product([J1,J2])
3124 sage: pi_left = J.cartesian_projection(0)
3125 sage: pi_right = J.cartesian_projection(1)
3126 sage: pi_left(J.one()).to_vector()
3127 (1, 0)
3128 sage: pi_right(J.one()).to_vector()
3129 (1, 0, 0, 1)
3130 sage: J.one().to_vector()
3131 (1, 0, 1, 0, 0, 1)
3132
3133 TESTS:
3134
3135 The answer never changes::
3136
3137 sage: set_random_seed()
3138 sage: J1 = random_eja()
3139 sage: J2 = random_eja()
3140 sage: J = cartesian_product([J1,J2])
3141 sage: P0 = J.cartesian_projection(0)
3142 sage: P1 = J.cartesian_projection(0)
3143 sage: P0 == P1
3144 True
3145
3146 """
3147 offset = sum( self.cartesian_factor(k).dimension()
3148 for k in range(i) )
3149 Ji = self.cartesian_factor(i)
3150 Pi = self._module_morphism(lambda j: Ji.monomial(j - offset),
3151 codomain=Ji)
3152
3153 return FiniteDimensionalEJAOperator(self,Ji,Pi.matrix())
3154
3155 @cached_method
3156 def cartesian_embedding(self, i):
3157 r"""
3158 SETUP::
3159
3160 sage: from mjo.eja.eja_algebra import (random_eja,
3161 ....: JordanSpinEJA,
3162 ....: HadamardEJA,
3163 ....: RealSymmetricEJA)
3164
3165 EXAMPLES:
3166
3167 The embedding morphisms are Euclidean Jordan algebra
3168 operators::
3169
3170 sage: J1 = HadamardEJA(2)
3171 sage: J2 = RealSymmetricEJA(2)
3172 sage: J = cartesian_product([J1,J2])
3173 sage: J.cartesian_embedding(0)
3174 Linear operator between finite-dimensional Euclidean Jordan
3175 algebras represented by the matrix:
3176 [1 0]
3177 [0 1]
3178 [0 0]
3179 [0 0]
3180 [0 0]
3181 Domain: Euclidean Jordan algebra of dimension 2 over
3182 Algebraic Real Field
3183 Codomain: Euclidean Jordan algebra of dimension 2 over
3184 Algebraic Real Field (+) Euclidean Jordan algebra of
3185 dimension 3 over Algebraic Real Field
3186 sage: J.cartesian_embedding(1)
3187 Linear operator between finite-dimensional Euclidean Jordan
3188 algebras represented by the matrix:
3189 [0 0 0]
3190 [0 0 0]
3191 [1 0 0]
3192 [0 1 0]
3193 [0 0 1]
3194 Domain: Euclidean Jordan algebra of dimension 3 over
3195 Algebraic Real Field
3196 Codomain: Euclidean Jordan algebra of dimension 2 over
3197 Algebraic Real Field (+) Euclidean Jordan algebra of
3198 dimension 3 over Algebraic Real Field
3199
3200 The embeddings work the way you'd expect on the vector
3201 representation of an element::
3202
3203 sage: J1 = JordanSpinEJA(3)
3204 sage: J2 = RealSymmetricEJA(2)
3205 sage: J = cartesian_product([J1,J2])
3206 sage: iota_left = J.cartesian_embedding(0)
3207 sage: iota_right = J.cartesian_embedding(1)
3208 sage: iota_left(J1.zero()) == J.zero()
3209 True
3210 sage: iota_right(J2.zero()) == J.zero()
3211 True
3212 sage: J1.one().to_vector()
3213 (1, 0, 0)
3214 sage: iota_left(J1.one()).to_vector()
3215 (1, 0, 0, 0, 0, 0)
3216 sage: J2.one().to_vector()
3217 (1, 0, 1)
3218 sage: iota_right(J2.one()).to_vector()
3219 (0, 0, 0, 1, 0, 1)
3220 sage: J.one().to_vector()
3221 (1, 0, 0, 1, 0, 1)
3222
3223 TESTS:
3224
3225 The answer never changes::
3226
3227 sage: set_random_seed()
3228 sage: J1 = random_eja()
3229 sage: J2 = random_eja()
3230 sage: J = cartesian_product([J1,J2])
3231 sage: E0 = J.cartesian_embedding(0)
3232 sage: E1 = J.cartesian_embedding(0)
3233 sage: E0 == E1
3234 True
3235
3236 Composing a projection with the corresponding inclusion should
3237 produce the identity map, and mismatching them should produce
3238 the zero map::
3239
3240 sage: set_random_seed()
3241 sage: J1 = random_eja()
3242 sage: J2 = random_eja()
3243 sage: J = cartesian_product([J1,J2])
3244 sage: iota_left = J.cartesian_embedding(0)
3245 sage: iota_right = J.cartesian_embedding(1)
3246 sage: pi_left = J.cartesian_projection(0)
3247 sage: pi_right = J.cartesian_projection(1)
3248 sage: pi_left*iota_left == J1.one().operator()
3249 True
3250 sage: pi_right*iota_right == J2.one().operator()
3251 True
3252 sage: (pi_left*iota_right).is_zero()
3253 True
3254 sage: (pi_right*iota_left).is_zero()
3255 True
3256
3257 """
3258 offset = sum( self.cartesian_factor(k).dimension()
3259 for k in range(i) )
3260 Ji = self.cartesian_factor(i)
3261 Ei = Ji._module_morphism(lambda j: self.monomial(j + offset),
3262 codomain=self)
3263 return FiniteDimensionalEJAOperator(Ji,self,Ei.matrix())
3264
3265
3266
3267 FiniteDimensionalEJA.CartesianProduct = CartesianProductEJA
3268
3269 class RationalBasisCartesianProductEJA(CartesianProductEJA,
3270 RationalBasisEJA):
3271 r"""
3272 A separate class for products of algebras for which we know a
3273 rational basis.
3274
3275 SETUP::
3276
3277 sage: from mjo.eja.eja_algebra import (HadamardEJA,
3278 ....: JordanSpinEJA,
3279 ....: OctonionHermitianEJA,
3280 ....: RealSymmetricEJA)
3281
3282 EXAMPLES:
3283
3284 This gives us fast characteristic polynomial computations in
3285 product algebras, too::
3286
3287
3288 sage: J1 = JordanSpinEJA(2)
3289 sage: J2 = RealSymmetricEJA(3)
3290 sage: J = cartesian_product([J1,J2])
3291 sage: J.characteristic_polynomial_of().degree()
3292 5
3293 sage: J.rank()
3294 5
3295
3296 TESTS:
3297
3298 The ``cartesian_product()`` function only uses the first factor to
3299 decide where the result will live; thus we have to be careful to
3300 check that all factors do indeed have a `_rational_algebra` member
3301 before we try to access it::
3302
3303 sage: J1 = OctonionHermitianEJA(1) # no rational basis
3304 sage: J2 = HadamardEJA(2)
3305 sage: cartesian_product([J1,J2])
3306 Euclidean Jordan algebra of dimension 1 over Algebraic Real Field
3307 (+) Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
3308 sage: cartesian_product([J2,J1])
3309 Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
3310 (+) Euclidean Jordan algebra of dimension 1 over Algebraic Real Field
3311
3312 """
3313 def __init__(self, algebras, **kwargs):
3314 CartesianProductEJA.__init__(self, algebras, **kwargs)
3315
3316 self._rational_algebra = None
3317 if self.vector_space().base_field() is not QQ:
3318 if all( hasattr(r, "_rational_algebra") for r in algebras ):
3319 self._rational_algebra = cartesian_product([
3320 r._rational_algebra for r in algebras
3321 ])
3322
3323
3324 RationalBasisEJA.CartesianProduct = RationalBasisCartesianProductEJA
3325
3326 def random_eja(max_dimension=None, *args, **kwargs):
3327 # Use the ConcreteEJA default as the total upper bound (regardless
3328 # of any whether or not any individual factors set a lower limit).
3329 if max_dimension is None:
3330 max_dimension = ConcreteEJA._max_random_instance_dimension()
3331 J1 = ConcreteEJA.random_instance(max_dimension, *args, **kwargs)
3332
3333
3334 # Roll the dice to see if we attempt a Cartesian product.
3335 dice_roll = ZZ.random_element(len(ConcreteEJA.__subclasses__()) + 1)
3336 new_max_dimension = max_dimension - J1.dimension()
3337 if new_max_dimension == 0 or dice_roll != 0:
3338 # If it's already as big as we're willing to tolerate, just
3339 # return it and don't worry about Cartesian products.
3340 return J1
3341 else:
3342 # Use random_eja() again so we can get more than two factors
3343 # if the sub-call also Decides on a cartesian product.
3344 J2 = random_eja(new_max_dimension, *args, **kwargs)
3345 return cartesian_product([J1,J2])