2 Representations and constructions for Euclidean Jordan algebras.
4 A Euclidean Jordan algebra is a Jordan algebra that has some
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.
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
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.
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:
31 * :class:`BilinearFormEJA`
32 * :class:`RealSymmetricEJA`
33 * :class:`ComplexHermitianEJA`
34 * :class:`QuaternionHermitianEJA`
35 * :class:`OctonionHermitianEJA`
37 In addition to these, we provide two other example constructions,
39 * :class:`HadamardEJA`
42 The Jordan spin algebra is a bilinear form algebra where the bilinear
43 form is the identity. The Hadamard EJA is simply a Cartesian product
44 of one-dimensional spin algebras. And last but least, the trivial EJA
45 is exactly what you think it is; it could also be obtained by
46 constructing a dimension-zero instance of any of the other
47 algebras. Cartesian products of these are also supported using the
48 usual ``cartesian_product()`` function; as a result, we support (up to
49 isomorphism) all Euclidean Jordan algebras.
53 sage: from mjo.eja.eja_algebra import random_eja
58 Euclidean Jordan algebra of dimension...
61 from itertools
import repeat
63 from sage
.algebras
.quatalg
.quaternion_algebra
import QuaternionAlgebra
64 from sage
.categories
.magmatic_algebras
import MagmaticAlgebras
65 from sage
.categories
.sets_cat
import cartesian_product
66 from sage
.combinat
.free_module
import CombinatorialFreeModule
67 from sage
.matrix
.constructor
import matrix
68 from sage
.matrix
.matrix_space
import MatrixSpace
69 from sage
.misc
.cachefunc
import cached_method
70 from sage
.misc
.table
import table
71 from sage
.modules
.free_module
import FreeModule
, VectorSpace
72 from sage
.rings
.all
import (ZZ
, QQ
, AA
, QQbar
, RR
, RLF
, CLF
,
75 from mjo
.eja
.eja_element
import FiniteDimensionalEJAElement
76 from mjo
.eja
.eja_operator
import FiniteDimensionalEJAOperator
77 from mjo
.eja
.eja_utils
import _all2list
, _mat2vec
79 class FiniteDimensionalEJA(CombinatorialFreeModule
):
81 A finite-dimensional Euclidean Jordan algebra.
85 - ``basis`` -- a tuple; a tuple of basis elements in "matrix
86 form," which must be the same form as the arguments to
87 ``jordan_product`` and ``inner_product``. In reality, "matrix
88 form" can be either vectors, matrices, or a Cartesian product
89 (ordered tuple) of vectors or matrices. All of these would
90 ideally be vector spaces in sage with no special-casing
91 needed; but in reality we turn vectors into column-matrices
92 and Cartesian products `(a,b)` into column matrices
93 `(a,b)^{T}` after converting `a` and `b` themselves.
95 - ``jordan_product`` -- a function; afunction of two ``basis``
96 elements (in matrix form) that returns their jordan product,
97 also in matrix form; this will be applied to ``basis`` to
98 compute a multiplication table for the algebra.
100 - ``inner_product`` -- a function; a function of two ``basis``
101 elements (in matrix form) that returns their inner
102 product. This will be applied to ``basis`` to compute an
103 inner-product table (basically a matrix) for this algebra.
105 - ``field`` -- a subfield of the reals (default: ``AA``); the scalar
106 field for the algebra.
108 - ``orthonormalize`` -- boolean (default: ``True``); whether or
109 not to orthonormalize the basis. Doing so is expensive and
110 generally rules out using the rationals as your ``field``, but
111 is required for spectral decompositions.
115 sage: from mjo.eja.eja_algebra import random_eja
119 We should compute that an element subalgebra is associative even
120 if we circumvent the element method::
122 sage: set_random_seed()
123 sage: J = random_eja(field=QQ,orthonormalize=False)
124 sage: x = J.random_element()
125 sage: A = x.subalgebra_generated_by(orthonormalize=False)
126 sage: basis = tuple(b.superalgebra_element() for b in A.basis())
127 sage: J.subalgebra(basis, orthonormalize=False).is_associative()
131 Element
= FiniteDimensionalEJAElement
140 cartesian_product
=False,
148 if not field
.is_subring(RR
):
149 # Note: this does return true for the real algebraic
150 # field, the rationals, and any quadratic field where
151 # we've specified a real embedding.
152 raise ValueError("scalar field is not real")
155 # Check commutativity of the Jordan and inner-products.
156 # This has to be done before we build the multiplication
157 # and inner-product tables/matrices, because we take
158 # advantage of symmetry in the process.
159 if not all( jordan_product(bi
,bj
) == jordan_product(bj
,bi
)
162 raise ValueError("Jordan product is not commutative")
164 if not all( inner_product(bi
,bj
) == inner_product(bj
,bi
)
167 raise ValueError("inner-product is not commutative")
170 category
= MagmaticAlgebras(field
).FiniteDimensional()
171 category
= category
.WithBasis().Unital().Commutative()
173 if associative
is None:
174 # We should figure it out. As with check_axioms, we have to do
175 # this without the help of the _jordan_product_is_associative()
176 # method because we need to know the category before we
177 # initialize the algebra.
178 associative
= all( jordan_product(jordan_product(bi
,bj
),bk
)
180 jordan_product(bi
,jordan_product(bj
,bk
))
186 # Element subalgebras can take advantage of this.
187 category
= category
.Associative()
188 if cartesian_product
:
189 # Use join() here because otherwise we only get the
190 # "Cartesian product of..." and not the things themselves.
191 category
= category
.join([category
,
192 category
.CartesianProducts()])
194 # Call the superclass constructor so that we can use its from_vector()
195 # method to build our multiplication table.
196 CombinatorialFreeModule
.__init
__(self
,
203 # Now comes all of the hard work. We'll be constructing an
204 # ambient vector space V that our (vectorized) basis lives in,
205 # as well as a subspace W of V spanned by those (vectorized)
206 # basis elements. The W-coordinates are the coefficients that
207 # we see in things like x = 1*b1 + 2*b2.
212 degree
= len(_all2list(basis
[0]))
214 # Build an ambient space that fits our matrix basis when
215 # written out as "long vectors."
216 V
= VectorSpace(field
, degree
)
218 # The matrix that will hole the orthonormal -> unorthonormal
219 # coordinate transformation.
220 self
._deortho
_matrix
= None
223 # Save a copy of the un-orthonormalized basis for later.
224 # Convert it to ambient V (vector) coordinates while we're
225 # at it, because we'd have to do it later anyway.
226 deortho_vector_basis
= tuple( V(_all2list(b
)) for b
in basis
)
228 from mjo
.eja
.eja_utils
import gram_schmidt
229 basis
= tuple(gram_schmidt(basis
, inner_product
))
231 # Save the (possibly orthonormalized) matrix basis for
233 self
._matrix
_basis
= basis
235 # Now create the vector space for the algebra, which will have
236 # its own set of non-ambient coordinates (in terms of the
238 vector_basis
= tuple( V(_all2list(b
)) for b
in basis
)
239 W
= V
.span_of_basis( vector_basis
, check
=check_axioms
)
242 # Now "W" is the vector space of our algebra coordinates. The
243 # variables "X1", "X2",... refer to the entries of vectors in
244 # W. Thus to convert back and forth between the orthonormal
245 # coordinates and the given ones, we need to stick the original
247 U
= V
.span_of_basis( deortho_vector_basis
, check
=check_axioms
)
248 self
._deortho
_matrix
= matrix( U
.coordinate_vector(q
)
249 for q
in vector_basis
)
252 # Now we actually compute the multiplication and inner-product
253 # tables/matrices using the possibly-orthonormalized basis.
254 self
._inner
_product
_matrix
= matrix
.identity(field
, n
)
255 self
._multiplication
_table
= [ [0 for j
in range(i
+1)]
258 # Note: the Jordan and inner-products are defined in terms
259 # of the ambient basis. It's important that their arguments
260 # are in ambient coordinates as well.
263 # ortho basis w.r.t. ambient coords
267 # The jordan product returns a matrixy answer, so we
268 # have to convert it to the algebra coordinates.
269 elt
= jordan_product(q_i
, q_j
)
270 elt
= W
.coordinate_vector(V(_all2list(elt
)))
271 self
._multiplication
_table
[i
][j
] = self
.from_vector(elt
)
273 if not orthonormalize
:
274 # If we're orthonormalizing the basis with respect
275 # to an inner-product, then the inner-product
276 # matrix with respect to the resulting basis is
277 # just going to be the identity.
278 ip
= inner_product(q_i
, q_j
)
279 self
._inner
_product
_matrix
[i
,j
] = ip
280 self
._inner
_product
_matrix
[j
,i
] = ip
282 self
._inner
_product
_matrix
._cache
= {'hermitian': True}
283 self
._inner
_product
_matrix
.set_immutable()
286 if not self
._is
_jordanian
():
287 raise ValueError("Jordan identity does not hold")
288 if not self
._inner
_product
_is
_associative
():
289 raise ValueError("inner product is not associative")
292 def _coerce_map_from_base_ring(self
):
294 Disable the map from the base ring into the algebra.
296 Performing a nonsense conversion like this automatically
297 is counterpedagogical. The fallback is to try the usual
298 element constructor, which should also fail.
302 sage: from mjo.eja.eja_algebra import random_eja
306 sage: set_random_seed()
307 sage: J = random_eja()
309 Traceback (most recent call last):
311 ValueError: not an element of this algebra
317 def product_on_basis(self
, i
, j
):
319 Returns the Jordan product of the `i` and `j`th basis elements.
321 This completely defines the Jordan product on the algebra, and
322 is used direclty by our superclass machinery to implement
327 sage: from mjo.eja.eja_algebra import random_eja
331 sage: set_random_seed()
332 sage: J = random_eja()
333 sage: n = J.dimension()
336 sage: bi_bj = J.zero()*J.zero()
338 ....: i = ZZ.random_element(n)
339 ....: j = ZZ.random_element(n)
340 ....: bi = J.monomial(i)
341 ....: bj = J.monomial(j)
342 ....: bi_bj = J.product_on_basis(i,j)
347 # We only stored the lower-triangular portion of the
348 # multiplication table.
350 return self
._multiplication
_table
[i
][j
]
352 return self
._multiplication
_table
[j
][i
]
354 def inner_product(self
, x
, y
):
356 The inner product associated with this Euclidean Jordan algebra.
358 Defaults to the trace inner product, but can be overridden by
359 subclasses if they are sure that the necessary properties are
364 sage: from mjo.eja.eja_algebra import (random_eja,
366 ....: BilinearFormEJA)
370 Our inner product is "associative," which means the following for
371 a symmetric bilinear form::
373 sage: set_random_seed()
374 sage: J = random_eja()
375 sage: x,y,z = J.random_elements(3)
376 sage: (x*y).inner_product(z) == y.inner_product(x*z)
381 Ensure that this is the usual inner product for the algebras
384 sage: set_random_seed()
385 sage: J = HadamardEJA.random_instance()
386 sage: x,y = J.random_elements(2)
387 sage: actual = x.inner_product(y)
388 sage: expected = x.to_vector().inner_product(y.to_vector())
389 sage: actual == expected
392 Ensure that this is one-half of the trace inner-product in a
393 BilinearFormEJA that isn't just the reals (when ``n`` isn't
394 one). This is in Faraut and Koranyi, and also my "On the
397 sage: set_random_seed()
398 sage: J = BilinearFormEJA.random_instance()
399 sage: n = J.dimension()
400 sage: x = J.random_element()
401 sage: y = J.random_element()
402 sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
406 B
= self
._inner
_product
_matrix
407 return (B
*x
.to_vector()).inner_product(y
.to_vector())
410 def is_associative(self
):
412 Return whether or not this algebra's Jordan product is associative.
416 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
420 sage: J = ComplexHermitianEJA(3, field=QQ, orthonormalize=False)
421 sage: J.is_associative()
423 sage: x = sum(J.gens())
424 sage: A = x.subalgebra_generated_by(orthonormalize=False)
425 sage: A.is_associative()
429 return "Associative" in self
.category().axioms()
431 def _is_commutative(self
):
433 Whether or not this algebra's multiplication table is commutative.
435 This method should of course always return ``True``, unless
436 this algebra was constructed with ``check_axioms=False`` and
437 passed an invalid multiplication table.
439 return all( x
*y
== y
*x
for x
in self
.gens() for y
in self
.gens() )
441 def _is_jordanian(self
):
443 Whether or not this algebra's multiplication table respects the
444 Jordan identity `(x^{2})(xy) = x(x^{2}y)`.
446 We only check one arrangement of `x` and `y`, so for a
447 ``True`` result to be truly true, you should also check
448 :meth:`_is_commutative`. This method should of course always
449 return ``True``, unless this algebra was constructed with
450 ``check_axioms=False`` and passed an invalid multiplication table.
452 return all( (self
.monomial(i
)**2)*(self
.monomial(i
)*self
.monomial(j
))
454 (self
.monomial(i
))*((self
.monomial(i
)**2)*self
.monomial(j
))
455 for i
in range(self
.dimension())
456 for j
in range(self
.dimension()) )
458 def _jordan_product_is_associative(self
):
460 Return whether or not this algebra's Jordan product is
461 associative; that is, whether or not `x*(y*z) = (x*y)*z`
464 This method should agree with :meth:`is_associative` unless
465 you lied about the value of the ``associative`` parameter
466 when you constructed the algebra.
470 sage: from mjo.eja.eja_algebra import (random_eja,
471 ....: RealSymmetricEJA,
472 ....: ComplexHermitianEJA,
473 ....: QuaternionHermitianEJA)
477 sage: J = RealSymmetricEJA(4, orthonormalize=False)
478 sage: J._jordan_product_is_associative()
480 sage: x = sum(J.gens())
481 sage: A = x.subalgebra_generated_by()
482 sage: A._jordan_product_is_associative()
487 sage: J = ComplexHermitianEJA(2,field=QQ,orthonormalize=False)
488 sage: J._jordan_product_is_associative()
490 sage: x = sum(J.gens())
491 sage: A = x.subalgebra_generated_by(orthonormalize=False)
492 sage: A._jordan_product_is_associative()
497 sage: J = QuaternionHermitianEJA(2)
498 sage: J._jordan_product_is_associative()
500 sage: x = sum(J.gens())
501 sage: A = x.subalgebra_generated_by()
502 sage: A._jordan_product_is_associative()
507 The values we've presupplied to the constructors agree with
510 sage: set_random_seed()
511 sage: J = random_eja()
512 sage: J.is_associative() == J._jordan_product_is_associative()
518 # Used to check whether or not something is zero.
521 # I don't know of any examples that make this magnitude
522 # necessary because I don't know how to make an
523 # associative algebra when the element subalgebra
524 # construction is unreliable (as it is over RDF; we can't
525 # find the degree of an element because we can't compute
526 # the rank of a matrix). But even multiplication of floats
527 # is non-associative, so *some* epsilon is needed... let's
528 # just take the one from _inner_product_is_associative?
531 for i
in range(self
.dimension()):
532 for j
in range(self
.dimension()):
533 for k
in range(self
.dimension()):
537 diff
= (x
*y
)*z
- x
*(y
*z
)
539 if diff
.norm() > epsilon
:
544 def _inner_product_is_associative(self
):
546 Return whether or not this algebra's inner product `B` is
547 associative; that is, whether or not `B(xy,z) = B(x,yz)`.
549 This method should of course always return ``True``, unless
550 this algebra was constructed with ``check_axioms=False`` and
551 passed an invalid Jordan or inner-product.
555 # Used to check whether or not something is zero.
558 # This choice is sufficient to allow the construction of
559 # QuaternionHermitianEJA(2, field=RDF) with check_axioms=True.
562 for i
in range(self
.dimension()):
563 for j
in range(self
.dimension()):
564 for k
in range(self
.dimension()):
568 diff
= (x
*y
).inner_product(z
) - x
.inner_product(y
*z
)
570 if diff
.abs() > epsilon
:
575 def _element_constructor_(self
, elt
):
577 Construct an element of this algebra from its vector or matrix
580 This gets called only after the parent element _call_ method
581 fails to find a coercion for the argument.
585 sage: from mjo.eja.eja_algebra import (random_eja,
588 ....: RealSymmetricEJA)
592 The identity in `S^n` is converted to the identity in the EJA::
594 sage: J = RealSymmetricEJA(3)
595 sage: I = matrix.identity(QQ,3)
596 sage: J(I) == J.one()
599 This skew-symmetric matrix can't be represented in the EJA::
601 sage: J = RealSymmetricEJA(3)
602 sage: A = matrix(QQ,3, lambda i,j: i-j)
604 Traceback (most recent call last):
606 ValueError: not an element of this algebra
608 Tuples work as well, provided that the matrix basis for the
609 algebra consists of them::
611 sage: J1 = HadamardEJA(3)
612 sage: J2 = RealSymmetricEJA(2)
613 sage: J = cartesian_product([J1,J2])
614 sage: J( (J1.matrix_basis()[1], J2.matrix_basis()[2]) )
619 Ensure that we can convert any element back and forth
620 faithfully between its matrix and algebra representations::
622 sage: set_random_seed()
623 sage: J = random_eja()
624 sage: x = J.random_element()
625 sage: J(x.to_matrix()) == x
628 We cannot coerce elements between algebras just because their
629 matrix representations are compatible::
631 sage: J1 = HadamardEJA(3)
632 sage: J2 = JordanSpinEJA(3)
634 Traceback (most recent call last):
636 ValueError: not an element of this algebra
638 Traceback (most recent call last):
640 ValueError: not an element of this algebra
642 msg
= "not an element of this algebra"
643 if elt
in self
.base_ring():
644 # Ensure that no base ring -> algebra coercion is performed
645 # by this method. There's some stupidity in sage that would
646 # otherwise propagate to this method; for example, sage thinks
647 # that the integer 3 belongs to the space of 2-by-2 matrices.
648 raise ValueError(msg
)
651 # Try to convert a vector into a column-matrix...
653 except (AttributeError, TypeError):
654 # and ignore failure, because we weren't really expecting
655 # a vector as an argument anyway.
658 if elt
not in self
.matrix_space():
659 raise ValueError(msg
)
661 # Thanks for nothing! Matrix spaces aren't vector spaces in
662 # Sage, so we have to figure out its matrix-basis coordinates
663 # ourselves. We use the basis space's ring instead of the
664 # element's ring because the basis space might be an algebraic
665 # closure whereas the base ring of the 3-by-3 identity matrix
666 # could be QQ instead of QQbar.
668 # And, we also have to handle Cartesian product bases (when
669 # the matrix basis consists of tuples) here. The "good news"
670 # is that we're already converting everything to long vectors,
671 # and that strategy works for tuples as well.
673 # We pass check=False because the matrix basis is "guaranteed"
674 # to be linearly independent... right? Ha ha.
676 V
= VectorSpace(self
.base_ring(), len(elt
))
677 W
= V
.span_of_basis( (V(_all2list(s
)) for s
in self
.matrix_basis()),
681 coords
= W
.coordinate_vector(V(elt
))
682 except ArithmeticError: # vector is not in free module
683 raise ValueError(msg
)
685 return self
.from_vector(coords
)
689 Return a string representation of ``self``.
693 sage: from mjo.eja.eja_algebra import JordanSpinEJA
697 Ensure that it says what we think it says::
699 sage: JordanSpinEJA(2, field=AA)
700 Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
701 sage: JordanSpinEJA(3, field=RDF)
702 Euclidean Jordan algebra of dimension 3 over Real Double Field
705 fmt
= "Euclidean Jordan algebra of dimension {} over {}"
706 return fmt
.format(self
.dimension(), self
.base_ring())
710 def characteristic_polynomial_of(self
):
712 Return the algebra's "characteristic polynomial of" function,
713 which is itself a multivariate polynomial that, when evaluated
714 at the coordinates of some algebra element, returns that
715 element's characteristic polynomial.
717 The resulting polynomial has `n+1` variables, where `n` is the
718 dimension of this algebra. The first `n` variables correspond to
719 the coordinates of an algebra element: when evaluated at the
720 coordinates of an algebra element with respect to a certain
721 basis, the result is a univariate polynomial (in the one
722 remaining variable ``t``), namely the characteristic polynomial
727 sage: from mjo.eja.eja_algebra import JordanSpinEJA, TrivialEJA
731 The characteristic polynomial in the spin algebra is given in
732 Alizadeh, Example 11.11::
734 sage: J = JordanSpinEJA(3)
735 sage: p = J.characteristic_polynomial_of(); p
736 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
737 sage: xvec = J.one().to_vector()
741 By definition, the characteristic polynomial is a monic
742 degree-zero polynomial in a rank-zero algebra. Note that
743 Cayley-Hamilton is indeed satisfied since the polynomial
744 ``1`` evaluates to the identity element of the algebra on
747 sage: J = TrivialEJA()
748 sage: J.characteristic_polynomial_of()
755 # The list of coefficient polynomials a_0, a_1, a_2, ..., a_(r-1).
756 a
= self
._charpoly
_coefficients
()
758 # We go to a bit of trouble here to reorder the
759 # indeterminates, so that it's easier to evaluate the
760 # characteristic polynomial at x's coordinates and get back
761 # something in terms of t, which is what we want.
762 S
= PolynomialRing(self
.base_ring(),'t')
766 S
= PolynomialRing(S
, R
.variable_names())
769 return (t
**r
+ sum( a
[k
]*(t
**k
) for k
in range(r
) ))
771 def coordinate_polynomial_ring(self
):
773 The multivariate polynomial ring in which this algebra's
774 :meth:`characteristic_polynomial_of` lives.
778 sage: from mjo.eja.eja_algebra import (HadamardEJA,
779 ....: RealSymmetricEJA)
783 sage: J = HadamardEJA(2)
784 sage: J.coordinate_polynomial_ring()
785 Multivariate Polynomial Ring in X1, X2...
786 sage: J = RealSymmetricEJA(3,field=QQ,orthonormalize=False)
787 sage: J.coordinate_polynomial_ring()
788 Multivariate Polynomial Ring in X1, X2, X3, X4, X5, X6...
791 var_names
= tuple( "X%d" % z
for z
in range(1, self
.dimension()+1) )
792 return PolynomialRing(self
.base_ring(), var_names
)
794 def inner_product(self
, x
, y
):
796 The inner product associated with this Euclidean Jordan algebra.
798 Defaults to the trace inner product, but can be overridden by
799 subclasses if they are sure that the necessary properties are
804 sage: from mjo.eja.eja_algebra import (random_eja,
806 ....: BilinearFormEJA)
810 Our inner product is "associative," which means the following for
811 a symmetric bilinear form::
813 sage: set_random_seed()
814 sage: J = random_eja()
815 sage: x,y,z = J.random_elements(3)
816 sage: (x*y).inner_product(z) == y.inner_product(x*z)
821 Ensure that this is the usual inner product for the algebras
824 sage: set_random_seed()
825 sage: J = HadamardEJA.random_instance()
826 sage: x,y = J.random_elements(2)
827 sage: actual = x.inner_product(y)
828 sage: expected = x.to_vector().inner_product(y.to_vector())
829 sage: actual == expected
832 Ensure that this is one-half of the trace inner-product in a
833 BilinearFormEJA that isn't just the reals (when ``n`` isn't
834 one). This is in Faraut and Koranyi, and also my "On the
837 sage: set_random_seed()
838 sage: J = BilinearFormEJA.random_instance()
839 sage: n = J.dimension()
840 sage: x = J.random_element()
841 sage: y = J.random_element()
842 sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
845 B
= self
._inner
_product
_matrix
846 return (B
*x
.to_vector()).inner_product(y
.to_vector())
849 def is_trivial(self
):
851 Return whether or not this algebra is trivial.
853 A trivial algebra contains only the zero element.
857 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
862 sage: J = ComplexHermitianEJA(3)
868 sage: J = TrivialEJA()
873 return self
.dimension() == 0
876 def multiplication_table(self
):
878 Return a visual representation of this algebra's multiplication
879 table (on basis elements).
883 sage: from mjo.eja.eja_algebra import JordanSpinEJA
887 sage: J = JordanSpinEJA(4)
888 sage: J.multiplication_table()
889 +----++----+----+----+----+
890 | * || b0 | b1 | b2 | b3 |
891 +====++====+====+====+====+
892 | b0 || b0 | b1 | b2 | b3 |
893 +----++----+----+----+----+
894 | b1 || b1 | b0 | 0 | 0 |
895 +----++----+----+----+----+
896 | b2 || b2 | 0 | b0 | 0 |
897 +----++----+----+----+----+
898 | b3 || b3 | 0 | 0 | b0 |
899 +----++----+----+----+----+
903 # Prepend the header row.
904 M
= [["*"] + list(self
.gens())]
906 # And to each subsequent row, prepend an entry that belongs to
907 # the left-side "header column."
908 M
+= [ [self
.monomial(i
)] + [ self
.monomial(i
)*self
.monomial(j
)
912 return table(M
, header_row
=True, header_column
=True, frame
=True)
915 def matrix_basis(self
):
917 Return an (often more natural) representation of this algebras
918 basis as an ordered tuple of matrices.
920 Every finite-dimensional Euclidean Jordan Algebra is a, up to
921 Jordan isomorphism, a direct sum of five simple
922 algebras---four of which comprise Hermitian matrices. And the
923 last type of algebra can of course be thought of as `n`-by-`1`
924 column matrices (ambiguusly called column vectors) to avoid
925 special cases. As a result, matrices (and column vectors) are
926 a natural representation format for Euclidean Jordan algebra
929 But, when we construct an algebra from a basis of matrices,
930 those matrix representations are lost in favor of coordinate
931 vectors *with respect to* that basis. We could eventually
932 convert back if we tried hard enough, but having the original
933 representations handy is valuable enough that we simply store
934 them and return them from this method.
936 Why implement this for non-matrix algebras? Avoiding special
937 cases for the :class:`BilinearFormEJA` pays with simplicity in
938 its own right. But mainly, we would like to be able to assume
939 that elements of a :class:`CartesianProductEJA` can be displayed
940 nicely, without having to have special classes for direct sums
941 one of whose components was a matrix algebra.
945 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
946 ....: RealSymmetricEJA)
950 sage: J = RealSymmetricEJA(2)
952 Finite family {0: b0, 1: b1, 2: b2}
953 sage: J.matrix_basis()
955 [1 0] [ 0 0.7071067811865475?] [0 0]
956 [0 0], [0.7071067811865475? 0], [0 1]
961 sage: J = JordanSpinEJA(2)
963 Finite family {0: b0, 1: b1}
964 sage: J.matrix_basis()
970 return self
._matrix
_basis
973 def matrix_space(self
):
975 Return the matrix space in which this algebra's elements live, if
976 we think of them as matrices (including column vectors of the
979 "By default" this will be an `n`-by-`1` column-matrix space,
980 except when the algebra is trivial. There it's `n`-by-`n`
981 (where `n` is zero), to ensure that two elements of the matrix
982 space (empty matrices) can be multiplied. For algebras of
983 matrices, this returns the space in which their
984 real embeddings live.
988 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
990 ....: QuaternionHermitianEJA,
995 By default, the matrix representation is just a column-matrix
996 equivalent to the vector representation::
998 sage: J = JordanSpinEJA(3)
999 sage: J.matrix_space()
1000 Full MatrixSpace of 3 by 1 dense matrices over Algebraic
1003 The matrix representation in the trivial algebra is
1004 zero-by-zero instead of the usual `n`-by-one::
1006 sage: J = TrivialEJA()
1007 sage: J.matrix_space()
1008 Full MatrixSpace of 0 by 0 dense matrices over Algebraic
1011 The matrix space for complex/quaternion Hermitian matrix EJA
1012 is the space in which their real-embeddings live, not the
1013 original complex/quaternion matrix space::
1015 sage: J = ComplexHermitianEJA(2,field=QQ,orthonormalize=False)
1016 sage: J.matrix_space()
1017 Full MatrixSpace of 4 by 4 dense matrices over Rational Field
1018 sage: J = QuaternionHermitianEJA(1,field=QQ,orthonormalize=False)
1019 sage: J.matrix_space()
1020 Full MatrixSpace of 4 by 4 dense matrices over Rational Field
1023 if self
.is_trivial():
1024 return MatrixSpace(self
.base_ring(), 0)
1026 return self
.matrix_basis()[0].parent()
1032 Return the unit element of this algebra.
1036 sage: from mjo.eja.eja_algebra import (HadamardEJA,
1041 We can compute unit element in the Hadamard EJA::
1043 sage: J = HadamardEJA(5)
1045 b0 + b1 + b2 + b3 + b4
1047 The unit element in the Hadamard EJA is inherited in the
1048 subalgebras generated by its elements::
1050 sage: J = HadamardEJA(5)
1052 b0 + b1 + b2 + b3 + b4
1053 sage: x = sum(J.gens())
1054 sage: A = x.subalgebra_generated_by(orthonormalize=False)
1057 sage: A.one().superalgebra_element()
1058 b0 + b1 + b2 + b3 + b4
1062 The identity element acts like the identity, regardless of
1063 whether or not we orthonormalize::
1065 sage: set_random_seed()
1066 sage: J = random_eja()
1067 sage: x = J.random_element()
1068 sage: J.one()*x == x and x*J.one() == x
1070 sage: A = x.subalgebra_generated_by()
1071 sage: y = A.random_element()
1072 sage: A.one()*y == y and y*A.one() == y
1077 sage: set_random_seed()
1078 sage: J = random_eja(field=QQ, orthonormalize=False)
1079 sage: x = J.random_element()
1080 sage: J.one()*x == x and x*J.one() == x
1082 sage: A = x.subalgebra_generated_by(orthonormalize=False)
1083 sage: y = A.random_element()
1084 sage: A.one()*y == y and y*A.one() == y
1087 The matrix of the unit element's operator is the identity,
1088 regardless of the base field and whether or not we
1091 sage: set_random_seed()
1092 sage: J = random_eja()
1093 sage: actual = J.one().operator().matrix()
1094 sage: expected = matrix.identity(J.base_ring(), J.dimension())
1095 sage: actual == expected
1097 sage: x = J.random_element()
1098 sage: A = x.subalgebra_generated_by()
1099 sage: actual = A.one().operator().matrix()
1100 sage: expected = matrix.identity(A.base_ring(), A.dimension())
1101 sage: actual == expected
1106 sage: set_random_seed()
1107 sage: J = random_eja(field=QQ, orthonormalize=False)
1108 sage: actual = J.one().operator().matrix()
1109 sage: expected = matrix.identity(J.base_ring(), J.dimension())
1110 sage: actual == expected
1112 sage: x = J.random_element()
1113 sage: A = x.subalgebra_generated_by(orthonormalize=False)
1114 sage: actual = A.one().operator().matrix()
1115 sage: expected = matrix.identity(A.base_ring(), A.dimension())
1116 sage: actual == expected
1119 Ensure that the cached unit element (often precomputed by
1120 hand) agrees with the computed one::
1122 sage: set_random_seed()
1123 sage: J = random_eja()
1124 sage: cached = J.one()
1125 sage: J.one.clear_cache()
1126 sage: J.one() == cached
1131 sage: set_random_seed()
1132 sage: J = random_eja(field=QQ, orthonormalize=False)
1133 sage: cached = J.one()
1134 sage: J.one.clear_cache()
1135 sage: J.one() == cached
1139 # We can brute-force compute the matrices of the operators
1140 # that correspond to the basis elements of this algebra.
1141 # If some linear combination of those basis elements is the
1142 # algebra identity, then the same linear combination of
1143 # their matrices has to be the identity matrix.
1145 # Of course, matrices aren't vectors in sage, so we have to
1146 # appeal to the "long vectors" isometry.
1147 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
1149 # Now we use basic linear algebra to find the coefficients,
1150 # of the matrices-as-vectors-linear-combination, which should
1151 # work for the original algebra basis too.
1152 A
= matrix(self
.base_ring(), oper_vecs
)
1154 # We used the isometry on the left-hand side already, but we
1155 # still need to do it for the right-hand side. Recall that we
1156 # wanted something that summed to the identity matrix.
1157 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
1159 # Now if there's an identity element in the algebra, this
1160 # should work. We solve on the left to avoid having to
1161 # transpose the matrix "A".
1162 return self
.from_vector(A
.solve_left(b
))
1165 def peirce_decomposition(self
, c
):
1167 The Peirce decomposition of this algebra relative to the
1170 In the future, this can be extended to a complete system of
1171 orthogonal idempotents.
1175 - ``c`` -- an idempotent of this algebra.
1179 A triple (J0, J5, J1) containing two subalgebras and one subspace
1182 - ``J0`` -- the algebra on the eigenspace of ``c.operator()``
1183 corresponding to the eigenvalue zero.
1185 - ``J5`` -- the eigenspace (NOT a subalgebra) of ``c.operator()``
1186 corresponding to the eigenvalue one-half.
1188 - ``J1`` -- the algebra on the eigenspace of ``c.operator()``
1189 corresponding to the eigenvalue one.
1191 These are the only possible eigenspaces for that operator, and this
1192 algebra is a direct sum of them. The spaces ``J0`` and ``J1`` are
1193 orthogonal, and are subalgebras of this algebra with the appropriate
1198 sage: from mjo.eja.eja_algebra import random_eja, RealSymmetricEJA
1202 The canonical example comes from the symmetric matrices, which
1203 decompose into diagonal and off-diagonal parts::
1205 sage: J = RealSymmetricEJA(3)
1206 sage: C = matrix(QQ, [ [1,0,0],
1210 sage: J0,J5,J1 = J.peirce_decomposition(c)
1212 Euclidean Jordan algebra of dimension 1...
1214 Vector space of degree 6 and dimension 2...
1216 Euclidean Jordan algebra of dimension 3...
1217 sage: J0.one().to_matrix()
1221 sage: orig_df = AA.options.display_format
1222 sage: AA.options.display_format = 'radical'
1223 sage: J.from_vector(J5.basis()[0]).to_matrix()
1227 sage: J.from_vector(J5.basis()[1]).to_matrix()
1231 sage: AA.options.display_format = orig_df
1232 sage: J1.one().to_matrix()
1239 Every algebra decomposes trivially with respect to its identity
1242 sage: set_random_seed()
1243 sage: J = random_eja()
1244 sage: J0,J5,J1 = J.peirce_decomposition(J.one())
1245 sage: J0.dimension() == 0 and J5.dimension() == 0
1247 sage: J1.superalgebra() == J and J1.dimension() == J.dimension()
1250 The decomposition is into eigenspaces, and its components are
1251 therefore necessarily orthogonal. Moreover, the identity
1252 elements in the two subalgebras are the projections onto their
1253 respective subspaces of the superalgebra's identity element::
1255 sage: set_random_seed()
1256 sage: J = random_eja()
1257 sage: x = J.random_element()
1258 sage: if not J.is_trivial():
1259 ....: while x.is_nilpotent():
1260 ....: x = J.random_element()
1261 sage: c = x.subalgebra_idempotent()
1262 sage: J0,J5,J1 = J.peirce_decomposition(c)
1264 sage: for (w,y,z) in zip(J0.basis(), J5.basis(), J1.basis()):
1265 ....: w = w.superalgebra_element()
1266 ....: y = J.from_vector(y)
1267 ....: z = z.superalgebra_element()
1268 ....: ipsum += w.inner_product(y).abs()
1269 ....: ipsum += w.inner_product(z).abs()
1270 ....: ipsum += y.inner_product(z).abs()
1273 sage: J1(c) == J1.one()
1275 sage: J0(J.one() - c) == J0.one()
1279 if not c
.is_idempotent():
1280 raise ValueError("element is not idempotent: %s" % c
)
1282 # Default these to what they should be if they turn out to be
1283 # trivial, because eigenspaces_left() won't return eigenvalues
1284 # corresponding to trivial spaces (e.g. it returns only the
1285 # eigenspace corresponding to lambda=1 if you take the
1286 # decomposition relative to the identity element).
1287 trivial
= self
.subalgebra(())
1288 J0
= trivial
# eigenvalue zero
1289 J5
= VectorSpace(self
.base_ring(), 0) # eigenvalue one-half
1290 J1
= trivial
# eigenvalue one
1292 for (eigval
, eigspace
) in c
.operator().matrix().right_eigenspaces():
1293 if eigval
== ~
(self
.base_ring()(2)):
1296 gens
= tuple( self
.from_vector(b
) for b
in eigspace
.basis() )
1297 subalg
= self
.subalgebra(gens
, check_axioms
=False)
1303 raise ValueError("unexpected eigenvalue: %s" % eigval
)
1308 def random_element(self
, thorough
=False):
1310 Return a random element of this algebra.
1312 Our algebra superclass method only returns a linear
1313 combination of at most two basis elements. We instead
1314 want the vector space "random element" method that
1315 returns a more diverse selection.
1319 - ``thorough`` -- (boolean; default False) whether or not we
1320 should generate irrational coefficients for the random
1321 element when our base ring is irrational; this slows the
1322 algebra operations to a crawl, but any truly random method
1326 # For a general base ring... maybe we can trust this to do the
1327 # right thing? Unlikely, but.
1328 V
= self
.vector_space()
1329 v
= V
.random_element()
1331 if self
.base_ring() is AA
:
1332 # The "random element" method of the algebraic reals is
1333 # stupid at the moment, and only returns integers between
1334 # -2 and 2, inclusive:
1336 # https://trac.sagemath.org/ticket/30875
1338 # Instead, we implement our own "random vector" method,
1339 # and then coerce that into the algebra. We use the vector
1340 # space degree here instead of the dimension because a
1341 # subalgebra could (for example) be spanned by only two
1342 # vectors, each with five coordinates. We need to
1343 # generate all five coordinates.
1345 v
*= QQbar
.random_element().real()
1347 v
*= QQ
.random_element()
1349 return self
.from_vector(V
.coordinate_vector(v
))
1351 def random_elements(self
, count
, thorough
=False):
1353 Return ``count`` random elements as a tuple.
1357 - ``thorough`` -- (boolean; default False) whether or not we
1358 should generate irrational coefficients for the random
1359 elements when our base ring is irrational; this slows the
1360 algebra operations to a crawl, but any truly random method
1365 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1369 sage: J = JordanSpinEJA(3)
1370 sage: x,y,z = J.random_elements(3)
1371 sage: all( [ x in J, y in J, z in J ])
1373 sage: len( J.random_elements(10) ) == 10
1377 return tuple( self
.random_element(thorough
)
1378 for idx
in range(count
) )
1382 def _charpoly_coefficients(self
):
1384 The `r` polynomial coefficients of the "characteristic polynomial
1389 sage: from mjo.eja.eja_algebra import random_eja
1393 The theory shows that these are all homogeneous polynomials of
1396 sage: set_random_seed()
1397 sage: J = random_eja()
1398 sage: all(p.is_homogeneous() for p in J._charpoly_coefficients())
1402 n
= self
.dimension()
1403 R
= self
.coordinate_polynomial_ring()
1405 F
= R
.fraction_field()
1408 # From a result in my book, these are the entries of the
1409 # basis representation of L_x.
1410 return sum( vars[k
]*self
.monomial(k
).operator().matrix()[i
,j
]
1413 L_x
= matrix(F
, n
, n
, L_x_i_j
)
1416 if self
.rank
.is_in_cache():
1418 # There's no need to pad the system with redundant
1419 # columns if we *know* they'll be redundant.
1422 # Compute an extra power in case the rank is equal to
1423 # the dimension (otherwise, we would stop at x^(r-1)).
1424 x_powers
= [ (L_x
**k
)*self
.one().to_vector()
1425 for k
in range(n
+1) ]
1426 A
= matrix
.column(F
, x_powers
[:n
])
1427 AE
= A
.extended_echelon_form()
1434 # The theory says that only the first "r" coefficients are
1435 # nonzero, and they actually live in the original polynomial
1436 # ring and not the fraction field. We negate them because in
1437 # the actual characteristic polynomial, they get moved to the
1438 # other side where x^r lives. We don't bother to trim A_rref
1439 # down to a square matrix and solve the resulting system,
1440 # because the upper-left r-by-r portion of A_rref is
1441 # guaranteed to be the identity matrix, so e.g.
1443 # A_rref.solve_right(Y)
1445 # would just be returning Y.
1446 return (-E
*b
)[:r
].change_ring(R
)
1451 Return the rank of this EJA.
1453 This is a cached method because we know the rank a priori for
1454 all of the algebras we can construct. Thus we can avoid the
1455 expensive ``_charpoly_coefficients()`` call unless we truly
1456 need to compute the whole characteristic polynomial.
1460 sage: from mjo.eja.eja_algebra import (HadamardEJA,
1461 ....: JordanSpinEJA,
1462 ....: RealSymmetricEJA,
1463 ....: ComplexHermitianEJA,
1464 ....: QuaternionHermitianEJA,
1469 The rank of the Jordan spin algebra is always two::
1471 sage: JordanSpinEJA(2).rank()
1473 sage: JordanSpinEJA(3).rank()
1475 sage: JordanSpinEJA(4).rank()
1478 The rank of the `n`-by-`n` Hermitian real, complex, or
1479 quaternion matrices is `n`::
1481 sage: RealSymmetricEJA(4).rank()
1483 sage: ComplexHermitianEJA(3).rank()
1485 sage: QuaternionHermitianEJA(2).rank()
1490 Ensure that every EJA that we know how to construct has a
1491 positive integer rank, unless the algebra is trivial in
1492 which case its rank will be zero::
1494 sage: set_random_seed()
1495 sage: J = random_eja()
1499 sage: r > 0 or (r == 0 and J.is_trivial())
1502 Ensure that computing the rank actually works, since the ranks
1503 of all simple algebras are known and will be cached by default::
1505 sage: set_random_seed() # long time
1506 sage: J = random_eja() # long time
1507 sage: cached = J.rank() # long time
1508 sage: J.rank.clear_cache() # long time
1509 sage: J.rank() == cached # long time
1513 return len(self
._charpoly
_coefficients
())
1516 def subalgebra(self
, basis
, **kwargs
):
1518 Create a subalgebra of this algebra from the given basis.
1520 from mjo
.eja
.eja_subalgebra
import FiniteDimensionalEJASubalgebra
1521 return FiniteDimensionalEJASubalgebra(self
, basis
, **kwargs
)
1524 def vector_space(self
):
1526 Return the vector space that underlies this algebra.
1530 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1534 sage: J = RealSymmetricEJA(2)
1535 sage: J.vector_space()
1536 Vector space of dimension 3 over...
1539 return self
.zero().to_vector().parent().ambient_vector_space()
1543 class RationalBasisEJA(FiniteDimensionalEJA
):
1545 Algebras whose supplied basis elements have all rational entries.
1549 sage: from mjo.eja.eja_algebra import BilinearFormEJA
1553 The supplied basis is orthonormalized by default::
1555 sage: B = matrix(QQ, [[1, 0, 0], [0, 25, -32], [0, -32, 41]])
1556 sage: J = BilinearFormEJA(B)
1557 sage: J.matrix_basis()
1574 # Abuse the check_field parameter to check that the entries of
1575 # out basis (in ambient coordinates) are in the field QQ.
1576 if not all( all(b_i
in QQ
for b_i
in b
.list()) for b
in basis
):
1577 raise TypeError("basis not rational")
1579 super().__init
__(basis
,
1583 check_field
=check_field
,
1586 self
._rational
_algebra
= None
1588 # There's no point in constructing the extra algebra if this
1589 # one is already rational.
1591 # Note: the same Jordan and inner-products work here,
1592 # because they are necessarily defined with respect to
1593 # ambient coordinates and not any particular basis.
1594 self
._rational
_algebra
= FiniteDimensionalEJA(
1599 associative
=self
.is_associative(),
1600 orthonormalize
=False,
1605 def _charpoly_coefficients(self
):
1609 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
1610 ....: JordanSpinEJA)
1614 The base ring of the resulting polynomial coefficients is what
1615 it should be, and not the rationals (unless the algebra was
1616 already over the rationals)::
1618 sage: J = JordanSpinEJA(3)
1619 sage: J._charpoly_coefficients()
1620 (X1^2 - X2^2 - X3^2, -2*X1)
1621 sage: a0 = J._charpoly_coefficients()[0]
1623 Algebraic Real Field
1624 sage: a0.base_ring()
1625 Algebraic Real Field
1628 if self
._rational
_algebra
is None:
1629 # There's no need to construct *another* algebra over the
1630 # rationals if this one is already over the
1631 # rationals. Likewise, if we never orthonormalized our
1632 # basis, we might as well just use the given one.
1633 return super()._charpoly
_coefficients
()
1635 # Do the computation over the rationals. The answer will be
1636 # the same, because all we've done is a change of basis.
1637 # Then, change back from QQ to our real base ring
1638 a
= ( a_i
.change_ring(self
.base_ring())
1639 for a_i
in self
._rational
_algebra
._charpoly
_coefficients
() )
1641 if self
._deortho
_matrix
is None:
1642 # This can happen if our base ring was, say, AA and we
1643 # chose not to (or didn't need to) orthonormalize. It's
1644 # still faster to do the computations over QQ even if
1645 # the numbers in the boxes stay the same.
1648 # Otherwise, convert the coordinate variables back to the
1649 # deorthonormalized ones.
1650 R
= self
.coordinate_polynomial_ring()
1651 from sage
.modules
.free_module_element
import vector
1652 X
= vector(R
, R
.gens())
1653 BX
= self
._deortho
_matrix
*X
1655 subs_dict
= { X[i]: BX[i] for i in range(len(X)) }
1656 return tuple( a_i
.subs(subs_dict
) for a_i
in a
)
1658 class ConcreteEJA(FiniteDimensionalEJA
):
1660 A class for the Euclidean Jordan algebras that we know by name.
1662 These are the Jordan algebras whose basis, multiplication table,
1663 rank, and so on are known a priori. More to the point, they are
1664 the Euclidean Jordan algebras for which we are able to conjure up
1665 a "random instance."
1669 sage: from mjo.eja.eja_algebra import ConcreteEJA
1673 Our basis is normalized with respect to the algebra's inner
1674 product, unless we specify otherwise::
1676 sage: set_random_seed()
1677 sage: J = ConcreteEJA.random_instance()
1678 sage: all( b.norm() == 1 for b in J.gens() )
1681 Since our basis is orthonormal with respect to the algebra's inner
1682 product, and since we know that this algebra is an EJA, any
1683 left-multiplication operator's matrix will be symmetric because
1684 natural->EJA basis representation is an isometry and within the
1685 EJA the operator is self-adjoint by the Jordan axiom::
1687 sage: set_random_seed()
1688 sage: J = ConcreteEJA.random_instance()
1689 sage: x = J.random_element()
1690 sage: x.operator().is_self_adjoint()
1695 def _max_random_instance_size():
1697 Return an integer "size" that is an upper bound on the size of
1698 this algebra when it is used in a random test
1699 case. Unfortunately, the term "size" is ambiguous -- when
1700 dealing with `R^n` under either the Hadamard or Jordan spin
1701 product, the "size" refers to the dimension `n`. When dealing
1702 with a matrix algebra (real symmetric or complex/quaternion
1703 Hermitian), it refers to the size of the matrix, which is far
1704 less than the dimension of the underlying vector space.
1706 This method must be implemented in each subclass.
1708 raise NotImplementedError
1711 def random_instance(cls
, *args
, **kwargs
):
1713 Return a random instance of this type of algebra.
1715 This method should be implemented in each subclass.
1717 from sage
.misc
.prandom
import choice
1718 eja_class
= choice(cls
.__subclasses
__())
1720 # These all bubble up to the RationalBasisEJA superclass
1721 # constructor, so any (kw)args valid there are also valid
1723 return eja_class
.random_instance(*args
, **kwargs
)
1728 def jordan_product(X
,Y
):
1729 return (X
*Y
+ Y
*X
)/2
1732 def trace_inner_product(X
,Y
):
1734 A trace inner-product for matrices that aren't embedded in the
1735 reals. It takes MATRICES as arguments, not EJA elements.
1737 return (X
*Y
).trace().real()
1739 class RealEmbeddedMatrixEJA(MatrixEJA
):
1741 def dimension_over_reals():
1743 The dimension of this matrix's base ring over the reals.
1745 The reals are dimension one over themselves, obviously; that's
1746 just `\mathbb{R}^{1}`. Likewise, the complex numbers `a + bi`
1747 have dimension two. Finally, the quaternions have dimension
1748 four over the reals.
1750 This is used to determine the size of the matrix returned from
1751 :meth:`real_embed`, among other things.
1753 raise NotImplementedError
1756 def real_embed(cls
,M
):
1758 Embed the matrix ``M`` into a space of real matrices.
1760 The matrix ``M`` can have entries in any field at the moment:
1761 the real numbers, complex numbers, or quaternions. And although
1762 they are not a field, we can probably support octonions at some
1763 point, too. This function returns a real matrix that "acts like"
1764 the original with respect to matrix multiplication; i.e.
1766 real_embed(M*N) = real_embed(M)*real_embed(N)
1769 if M
.ncols() != M
.nrows():
1770 raise ValueError("the matrix 'M' must be square")
1775 def real_unembed(cls
,M
):
1777 The inverse of :meth:`real_embed`.
1779 if M
.ncols() != M
.nrows():
1780 raise ValueError("the matrix 'M' must be square")
1781 if not ZZ(M
.nrows()).mod(cls
.dimension_over_reals()).is_zero():
1782 raise ValueError("the matrix 'M' must be a real embedding")
1787 def trace_inner_product(cls
,X
,Y
):
1789 Compute the trace inner-product of two real-embeddings.
1793 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
1794 ....: QuaternionHermitianEJA)
1798 sage: set_random_seed()
1799 sage: J = ComplexHermitianEJA.random_instance()
1800 sage: x,y = J.random_elements(2)
1801 sage: Xe = x.to_matrix()
1802 sage: Ye = y.to_matrix()
1803 sage: X = J.real_unembed(Xe)
1804 sage: Y = J.real_unembed(Ye)
1805 sage: expected = (X*Y).trace().real()
1806 sage: actual = J.trace_inner_product(Xe,Ye)
1807 sage: actual == expected
1812 sage: set_random_seed()
1813 sage: J = QuaternionHermitianEJA.random_instance()
1814 sage: x,y = J.random_elements(2)
1815 sage: Xe = x.to_matrix()
1816 sage: Ye = y.to_matrix()
1817 sage: X = J.real_unembed(Xe)
1818 sage: Y = J.real_unembed(Ye)
1819 sage: expected = (X*Y).trace().coefficient_tuple()[0]
1820 sage: actual = J.trace_inner_product(Xe,Ye)
1821 sage: actual == expected
1825 # This does in fact compute the real part of the trace.
1826 # If we compute the trace of e.g. a complex matrix M,
1827 # then we do so by adding up its diagonal entries --
1828 # call them z_1 through z_n. The real embedding of z_1
1829 # will be a 2-by-2 REAL matrix [a, b; -b, a] whose trace
1830 # as a REAL matrix will be 2*a = 2*Re(z_1). And so forth.
1831 return (X
*Y
).trace()/cls
.dimension_over_reals()
1833 class RealSymmetricEJA(ConcreteEJA
, RationalBasisEJA
, MatrixEJA
):
1835 The rank-n simple EJA consisting of real symmetric n-by-n
1836 matrices, the usual symmetric Jordan product, and the trace inner
1837 product. It has dimension `(n^2 + n)/2` over the reals.
1841 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1845 sage: J = RealSymmetricEJA(2)
1846 sage: b0, b1, b2 = J.gens()
1854 In theory, our "field" can be any subfield of the reals::
1856 sage: RealSymmetricEJA(2, field=RDF, check_axioms=True)
1857 Euclidean Jordan algebra of dimension 3 over Real Double Field
1858 sage: RealSymmetricEJA(2, field=RR, check_axioms=True)
1859 Euclidean Jordan algebra of dimension 3 over Real Field with
1860 53 bits of precision
1864 The dimension of this algebra is `(n^2 + n) / 2`::
1866 sage: set_random_seed()
1867 sage: n_max = RealSymmetricEJA._max_random_instance_size()
1868 sage: n = ZZ.random_element(1, n_max)
1869 sage: J = RealSymmetricEJA(n)
1870 sage: J.dimension() == (n^2 + n)/2
1873 The Jordan multiplication is what we think it is::
1875 sage: set_random_seed()
1876 sage: J = RealSymmetricEJA.random_instance()
1877 sage: x,y = J.random_elements(2)
1878 sage: actual = (x*y).to_matrix()
1879 sage: X = x.to_matrix()
1880 sage: Y = y.to_matrix()
1881 sage: expected = (X*Y + Y*X)/2
1882 sage: actual == expected
1884 sage: J(expected) == x*y
1887 We can change the generator prefix::
1889 sage: RealSymmetricEJA(3, prefix='q').gens()
1890 (q0, q1, q2, q3, q4, q5)
1892 We can construct the (trivial) algebra of rank zero::
1894 sage: RealSymmetricEJA(0)
1895 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1899 def _denormalized_basis(cls
, n
, field
):
1901 Return a basis for the space of real symmetric n-by-n matrices.
1905 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1909 sage: set_random_seed()
1910 sage: n = ZZ.random_element(1,5)
1911 sage: B = RealSymmetricEJA._denormalized_basis(n,ZZ)
1912 sage: all( M.is_symmetric() for M in B)
1916 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1920 for j
in range(i
+1):
1921 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
1925 Sij
= Eij
+ Eij
.transpose()
1931 def _max_random_instance_size():
1932 return 4 # Dimension 10
1935 def random_instance(cls
, **kwargs
):
1937 Return a random instance of this type of algebra.
1939 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1940 return cls(n
, **kwargs
)
1942 def __init__(self
, n
, field
=AA
, **kwargs
):
1943 # We know this is a valid EJA, but will double-check
1944 # if the user passes check_axioms=True.
1945 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
1951 super().__init
__(self
._denormalized
_basis
(n
,field
),
1952 self
.jordan_product
,
1953 self
.trace_inner_product
,
1955 associative
=associative
,
1958 # TODO: this could be factored out somehow, but is left here
1959 # because the MatrixEJA is not presently a subclass of the
1960 # FDEJA class that defines rank() and one().
1961 self
.rank
.set_cache(n
)
1962 idV
= self
.matrix_space().one()
1963 self
.one
.set_cache(self(idV
))
1967 class ComplexMatrixEJA(RealEmbeddedMatrixEJA
):
1968 # A manual dictionary-cache for the complex_extension() method,
1969 # since apparently @classmethods can't also be @cached_methods.
1970 _complex_extension
= {}
1973 def complex_extension(cls
,field
):
1975 The complex field that we embed/unembed, as an extension
1976 of the given ``field``.
1978 if field
in cls
._complex
_extension
:
1979 return cls
._complex
_extension
[field
]
1981 # Sage doesn't know how to adjoin the complex "i" (the root of
1982 # x^2 + 1) to a field in a general way. Here, we just enumerate
1983 # all of the cases that I have cared to support so far.
1985 # Sage doesn't know how to embed AA into QQbar, i.e. how
1986 # to adjoin sqrt(-1) to AA.
1988 elif not field
.is_exact():
1990 F
= field
.complex_field()
1992 # Works for QQ and... maybe some other fields.
1993 R
= PolynomialRing(field
, 'z')
1995 F
= field
.extension(z
**2 + 1, 'I', embedding
=CLF(-1).sqrt())
1997 cls
._complex
_extension
[field
] = F
2001 def dimension_over_reals():
2005 def real_embed(cls
,M
):
2007 Embed the n-by-n complex matrix ``M`` into the space of real
2008 matrices of size 2n-by-2n via the map the sends each entry `z = a +
2009 bi` to the block matrix ``[[a,b],[-b,a]]``.
2013 sage: from mjo.eja.eja_algebra import ComplexMatrixEJA
2017 sage: F = QuadraticField(-1, 'I')
2018 sage: x1 = F(4 - 2*i)
2019 sage: x2 = F(1 + 2*i)
2022 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
2023 sage: ComplexMatrixEJA.real_embed(M)
2032 Embedding is a homomorphism (isomorphism, in fact)::
2034 sage: set_random_seed()
2035 sage: n = ZZ.random_element(3)
2036 sage: F = QuadraticField(-1, 'I')
2037 sage: X = random_matrix(F, n)
2038 sage: Y = random_matrix(F, n)
2039 sage: Xe = ComplexMatrixEJA.real_embed(X)
2040 sage: Ye = ComplexMatrixEJA.real_embed(Y)
2041 sage: XYe = ComplexMatrixEJA.real_embed(X*Y)
2046 super().real_embed(M
)
2049 # We don't need any adjoined elements...
2050 field
= M
.base_ring().base_ring()
2056 blocks
.append(matrix(field
, 2, [ [ a
, b
],
2059 return matrix
.block(field
, n
, blocks
)
2063 def real_unembed(cls
,M
):
2065 The inverse of _embed_complex_matrix().
2069 sage: from mjo.eja.eja_algebra import ComplexMatrixEJA
2073 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
2074 ....: [-2, 1, -4, 3],
2075 ....: [ 9, 10, 11, 12],
2076 ....: [-10, 9, -12, 11] ])
2077 sage: ComplexMatrixEJA.real_unembed(A)
2079 [ 10*I + 9 12*I + 11]
2083 Unembedding is the inverse of embedding::
2085 sage: set_random_seed()
2086 sage: F = QuadraticField(-1, 'I')
2087 sage: M = random_matrix(F, 3)
2088 sage: Me = ComplexMatrixEJA.real_embed(M)
2089 sage: ComplexMatrixEJA.real_unembed(Me) == M
2093 super().real_unembed(M
)
2095 d
= cls
.dimension_over_reals()
2096 F
= cls
.complex_extension(M
.base_ring())
2099 # Go top-left to bottom-right (reading order), converting every
2100 # 2-by-2 block we see to a single complex element.
2102 for k
in range(n
/d
):
2103 for j
in range(n
/d
):
2104 submat
= M
[d
*k
:d
*k
+d
,d
*j
:d
*j
+d
]
2105 if submat
[0,0] != submat
[1,1]:
2106 raise ValueError('bad on-diagonal submatrix')
2107 if submat
[0,1] != -submat
[1,0]:
2108 raise ValueError('bad off-diagonal submatrix')
2109 z
= submat
[0,0] + submat
[0,1]*i
2112 return matrix(F
, n
/d
, elements
)
2115 class ComplexHermitianEJA(ConcreteEJA
, RationalBasisEJA
, ComplexMatrixEJA
):
2117 The rank-n simple EJA consisting of complex Hermitian n-by-n
2118 matrices over the real numbers, the usual symmetric Jordan product,
2119 and the real-part-of-trace inner product. It has dimension `n^2` over
2124 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
2128 In theory, our "field" can be any subfield of the reals::
2130 sage: ComplexHermitianEJA(2, field=RDF, check_axioms=True)
2131 Euclidean Jordan algebra of dimension 4 over Real Double Field
2132 sage: ComplexHermitianEJA(2, field=RR, check_axioms=True)
2133 Euclidean Jordan algebra of dimension 4 over Real Field with
2134 53 bits of precision
2138 The dimension of this algebra is `n^2`::
2140 sage: set_random_seed()
2141 sage: n_max = ComplexHermitianEJA._max_random_instance_size()
2142 sage: n = ZZ.random_element(1, n_max)
2143 sage: J = ComplexHermitianEJA(n)
2144 sage: J.dimension() == n^2
2147 The Jordan multiplication is what we think it is::
2149 sage: set_random_seed()
2150 sage: J = ComplexHermitianEJA.random_instance()
2151 sage: x,y = J.random_elements(2)
2152 sage: actual = (x*y).to_matrix()
2153 sage: X = x.to_matrix()
2154 sage: Y = y.to_matrix()
2155 sage: expected = (X*Y + Y*X)/2
2156 sage: actual == expected
2158 sage: J(expected) == x*y
2161 We can change the generator prefix::
2163 sage: ComplexHermitianEJA(2, prefix='z').gens()
2166 We can construct the (trivial) algebra of rank zero::
2168 sage: ComplexHermitianEJA(0)
2169 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2174 def _denormalized_basis(cls
, n
, field
):
2176 Returns a basis for the space of complex Hermitian n-by-n matrices.
2178 Why do we embed these? Basically, because all of numerical linear
2179 algebra assumes that you're working with vectors consisting of `n`
2180 entries from a field and scalars from the same field. There's no way
2181 to tell SageMath that (for example) the vectors contain complex
2182 numbers, while the scalar field is real.
2186 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
2190 sage: set_random_seed()
2191 sage: n = ZZ.random_element(1,5)
2192 sage: B = ComplexHermitianEJA._denormalized_basis(n,ZZ)
2193 sage: all( M.is_symmetric() for M in B)
2197 R
= PolynomialRing(ZZ
, 'z')
2199 F
= ZZ
.extension(z
**2 + 1, 'I')
2202 # This is like the symmetric case, but we need to be careful:
2204 # * We want conjugate-symmetry, not just symmetry.
2205 # * The diagonal will (as a result) be real.
2208 Eij
= matrix
.zero(F
,n
)
2210 for j
in range(i
+1):
2214 Sij
= cls
.real_embed(Eij
)
2217 # The second one has a minus because it's conjugated.
2218 Eij
[j
,i
] = 1 # Eij = Eij + Eij.transpose()
2219 Sij_real
= cls
.real_embed(Eij
)
2221 # Eij = I*Eij - I*Eij.transpose()
2224 Sij_imag
= cls
.real_embed(Eij
)
2230 # Since we embedded the entries, we can drop back to the
2231 # desired real "field" instead of the extension "F".
2232 return tuple( s
.change_ring(field
) for s
in S
)
2235 def __init__(self
, n
, field
=AA
, **kwargs
):
2236 # We know this is a valid EJA, but will double-check
2237 # if the user passes check_axioms=True.
2238 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2244 super().__init
__(self
._denormalized
_basis
(n
,field
),
2245 self
.jordan_product
,
2246 self
.trace_inner_product
,
2248 associative
=associative
,
2250 # TODO: this could be factored out somehow, but is left here
2251 # because the MatrixEJA is not presently a subclass of the
2252 # FDEJA class that defines rank() and one().
2253 self
.rank
.set_cache(n
)
2254 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
2255 self
.one
.set_cache(self(idV
))
2258 def _max_random_instance_size():
2259 return 3 # Dimension 9
2262 def random_instance(cls
, **kwargs
):
2264 Return a random instance of this type of algebra.
2266 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2267 return cls(n
, **kwargs
)
2269 class QuaternionMatrixEJA(RealEmbeddedMatrixEJA
):
2271 # A manual dictionary-cache for the quaternion_extension() method,
2272 # since apparently @classmethods can't also be @cached_methods.
2273 _quaternion_extension
= {}
2276 def quaternion_extension(cls
,field
):
2278 The quaternion field that we embed/unembed, as an extension
2279 of the given ``field``.
2281 if field
in cls
._quaternion
_extension
:
2282 return cls
._quaternion
_extension
[field
]
2284 Q
= QuaternionAlgebra(field
,-1,-1)
2286 cls
._quaternion
_extension
[field
] = Q
2290 def dimension_over_reals():
2294 def real_embed(cls
,M
):
2296 Embed the n-by-n quaternion matrix ``M`` into the space of real
2297 matrices of size 4n-by-4n by first sending each quaternion entry `z
2298 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
2299 c+di],[-c + di, a-bi]]`, and then embedding those into a real
2304 sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA
2308 sage: Q = QuaternionAlgebra(QQ,-1,-1)
2309 sage: i,j,k = Q.gens()
2310 sage: x = 1 + 2*i + 3*j + 4*k
2311 sage: M = matrix(Q, 1, [[x]])
2312 sage: QuaternionMatrixEJA.real_embed(M)
2318 Embedding is a homomorphism (isomorphism, in fact)::
2320 sage: set_random_seed()
2321 sage: n = ZZ.random_element(2)
2322 sage: Q = QuaternionAlgebra(QQ,-1,-1)
2323 sage: X = random_matrix(Q, n)
2324 sage: Y = random_matrix(Q, n)
2325 sage: Xe = QuaternionMatrixEJA.real_embed(X)
2326 sage: Ye = QuaternionMatrixEJA.real_embed(Y)
2327 sage: XYe = QuaternionMatrixEJA.real_embed(X*Y)
2332 super().real_embed(M
)
2333 quaternions
= M
.base_ring()
2336 F
= QuadraticField(-1, 'I')
2341 t
= z
.coefficient_tuple()
2346 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
2347 [-c
+ d
*i
, a
- b
*i
]])
2348 realM
= ComplexMatrixEJA
.real_embed(cplxM
)
2349 blocks
.append(realM
)
2351 # We should have real entries by now, so use the realest field
2352 # we've got for the return value.
2353 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
2358 def real_unembed(cls
,M
):
2360 The inverse of _embed_quaternion_matrix().
2364 sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA
2368 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
2369 ....: [-2, 1, -4, 3],
2370 ....: [-3, 4, 1, -2],
2371 ....: [-4, -3, 2, 1]])
2372 sage: QuaternionMatrixEJA.real_unembed(M)
2373 [1 + 2*i + 3*j + 4*k]
2377 Unembedding is the inverse of embedding::
2379 sage: set_random_seed()
2380 sage: Q = QuaternionAlgebra(QQ, -1, -1)
2381 sage: M = random_matrix(Q, 3)
2382 sage: Me = QuaternionMatrixEJA.real_embed(M)
2383 sage: QuaternionMatrixEJA.real_unembed(Me) == M
2387 super().real_unembed(M
)
2389 d
= cls
.dimension_over_reals()
2391 # Use the base ring of the matrix to ensure that its entries can be
2392 # multiplied by elements of the quaternion algebra.
2393 Q
= cls
.quaternion_extension(M
.base_ring())
2396 # Go top-left to bottom-right (reading order), converting every
2397 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
2400 for l
in range(n
/d
):
2401 for m
in range(n
/d
):
2402 submat
= ComplexMatrixEJA
.real_unembed(
2403 M
[d
*l
:d
*l
+d
,d
*m
:d
*m
+d
] )
2404 if submat
[0,0] != submat
[1,1].conjugate():
2405 raise ValueError('bad on-diagonal submatrix')
2406 if submat
[0,1] != -submat
[1,0].conjugate():
2407 raise ValueError('bad off-diagonal submatrix')
2408 z
= submat
[0,0].real()
2409 z
+= submat
[0,0].imag()*i
2410 z
+= submat
[0,1].real()*j
2411 z
+= submat
[0,1].imag()*k
2414 return matrix(Q
, n
/d
, elements
)
2417 class QuaternionHermitianEJA(ConcreteEJA
,
2419 QuaternionMatrixEJA
):
2421 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
2422 matrices, the usual symmetric Jordan product, and the
2423 real-part-of-trace inner product. It has dimension `2n^2 - n` over
2428 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2432 In theory, our "field" can be any subfield of the reals::
2434 sage: QuaternionHermitianEJA(2, field=RDF, check_axioms=True)
2435 Euclidean Jordan algebra of dimension 6 over Real Double Field
2436 sage: QuaternionHermitianEJA(2, field=RR, check_axioms=True)
2437 Euclidean Jordan algebra of dimension 6 over Real Field with
2438 53 bits of precision
2442 The dimension of this algebra is `2*n^2 - n`::
2444 sage: set_random_seed()
2445 sage: n_max = QuaternionHermitianEJA._max_random_instance_size()
2446 sage: n = ZZ.random_element(1, n_max)
2447 sage: J = QuaternionHermitianEJA(n)
2448 sage: J.dimension() == 2*(n^2) - n
2451 The Jordan multiplication is what we think it is::
2453 sage: set_random_seed()
2454 sage: J = QuaternionHermitianEJA.random_instance()
2455 sage: x,y = J.random_elements(2)
2456 sage: actual = (x*y).to_matrix()
2457 sage: X = x.to_matrix()
2458 sage: Y = y.to_matrix()
2459 sage: expected = (X*Y + Y*X)/2
2460 sage: actual == expected
2462 sage: J(expected) == x*y
2465 We can change the generator prefix::
2467 sage: QuaternionHermitianEJA(2, prefix='a').gens()
2468 (a0, a1, a2, a3, a4, a5)
2470 We can construct the (trivial) algebra of rank zero::
2472 sage: QuaternionHermitianEJA(0)
2473 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2477 def _denormalized_basis(cls
, n
, field
):
2479 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
2481 Why do we embed these? Basically, because all of numerical
2482 linear algebra assumes that you're working with vectors consisting
2483 of `n` entries from a field and scalars from the same field. There's
2484 no way to tell SageMath that (for example) the vectors contain
2485 complex numbers, while the scalar field is real.
2489 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2493 sage: set_random_seed()
2494 sage: n = ZZ.random_element(1,5)
2495 sage: B = QuaternionHermitianEJA._denormalized_basis(n,ZZ)
2496 sage: all( M.is_symmetric() for M in B )
2500 Q
= QuaternionAlgebra(QQ
,-1,-1)
2503 # This is like the symmetric case, but we need to be careful:
2505 # * We want conjugate-symmetry, not just symmetry.
2506 # * The diagonal will (as a result) be real.
2509 Eij
= matrix
.zero(Q
,n
)
2511 for j
in range(i
+1):
2515 Sij
= cls
.real_embed(Eij
)
2518 # The second, third, and fourth ones have a minus
2519 # because they're conjugated.
2520 # Eij = Eij + Eij.transpose()
2522 Sij_real
= cls
.real_embed(Eij
)
2524 # Eij = I*(Eij - Eij.transpose())
2527 Sij_I
= cls
.real_embed(Eij
)
2529 # Eij = J*(Eij - Eij.transpose())
2532 Sij_J
= cls
.real_embed(Eij
)
2534 # Eij = K*(Eij - Eij.transpose())
2537 Sij_K
= cls
.real_embed(Eij
)
2543 # Since we embedded the entries, we can drop back to the
2544 # desired real "field" instead of the quaternion algebra "Q".
2545 return tuple( s
.change_ring(field
) for s
in S
)
2548 def __init__(self
, n
, field
=AA
, **kwargs
):
2549 # We know this is a valid EJA, but will double-check
2550 # if the user passes check_axioms=True.
2551 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2557 super().__init
__(self
._denormalized
_basis
(n
,field
),
2558 self
.jordan_product
,
2559 self
.trace_inner_product
,
2561 associative
=associative
,
2564 # TODO: this could be factored out somehow, but is left here
2565 # because the MatrixEJA is not presently a subclass of the
2566 # FDEJA class that defines rank() and one().
2567 self
.rank
.set_cache(n
)
2568 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
2569 self
.one
.set_cache(self(idV
))
2573 def _max_random_instance_size():
2575 The maximum rank of a random QuaternionHermitianEJA.
2577 return 2 # Dimension 6
2580 def random_instance(cls
, **kwargs
):
2582 Return a random instance of this type of algebra.
2584 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2585 return cls(n
, **kwargs
)
2587 class OctonionHermitianEJA(ConcreteEJA
, MatrixEJA
):
2591 sage: from mjo.eja.eja_algebra import (FiniteDimensionalEJA,
2592 ....: OctonionHermitianEJA)
2596 The 3-by-3 algebra satisfies the axioms of an EJA::
2598 sage: OctonionHermitianEJA(3, # long time
2599 ....: field=QQ, # long time
2600 ....: orthonormalize=False, # long time
2601 ....: check_axioms=True) # long time
2602 Euclidean Jordan algebra of dimension 27 over Rational Field
2604 After a change-of-basis, the 2-by-2 algebra has the same
2605 multiplication table as the ten-dimensional Jordan spin algebra::
2607 sage: b = OctonionHermitianEJA._denormalized_basis(2,QQ)
2608 sage: basis = (b[0] + b[9],) + b[1:9] + (b[0] - b[9],)
2609 sage: jp = OctonionHermitianEJA.jordan_product
2610 sage: ip = OctonionHermitianEJA.trace_inner_product
2611 sage: J = FiniteDimensionalEJA(basis,
2615 ....: orthonormalize=False)
2616 sage: J.multiplication_table()
2617 +----++----+----+----+----+----+----+----+----+----+----+
2618 | * || b0 | b1 | b2 | b3 | b4 | b5 | b6 | b7 | b8 | b9 |
2619 +====++====+====+====+====+====+====+====+====+====+====+
2620 | b0 || b0 | b1 | b2 | b3 | b4 | b5 | b6 | b7 | b8 | b9 |
2621 +----++----+----+----+----+----+----+----+----+----+----+
2622 | b1 || b1 | b0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
2623 +----++----+----+----+----+----+----+----+----+----+----+
2624 | b2 || b2 | 0 | b0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
2625 +----++----+----+----+----+----+----+----+----+----+----+
2626 | b3 || b3 | 0 | 0 | b0 | 0 | 0 | 0 | 0 | 0 | 0 |
2627 +----++----+----+----+----+----+----+----+----+----+----+
2628 | b4 || b4 | 0 | 0 | 0 | b0 | 0 | 0 | 0 | 0 | 0 |
2629 +----++----+----+----+----+----+----+----+----+----+----+
2630 | b5 || b5 | 0 | 0 | 0 | 0 | b0 | 0 | 0 | 0 | 0 |
2631 +----++----+----+----+----+----+----+----+----+----+----+
2632 | b6 || b6 | 0 | 0 | 0 | 0 | 0 | b0 | 0 | 0 | 0 |
2633 +----++----+----+----+----+----+----+----+----+----+----+
2634 | b7 || b7 | 0 | 0 | 0 | 0 | 0 | 0 | b0 | 0 | 0 |
2635 +----++----+----+----+----+----+----+----+----+----+----+
2636 | b8 || b8 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | b0 | 0 |
2637 +----++----+----+----+----+----+----+----+----+----+----+
2638 | b9 || b9 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | b0 |
2639 +----++----+----+----+----+----+----+----+----+----+----+
2643 We can actually construct the 27-dimensional Albert algebra,
2644 and we get the right unit element if we recompute it::
2646 sage: J = OctonionHermitianEJA(3, # long time
2647 ....: field=QQ, # long time
2648 ....: orthonormalize=False) # long time
2649 sage: J.one.clear_cache() # long time
2650 sage: J.one() # long time
2652 sage: J.one().to_matrix() # long time
2661 The 2-by-2 algebra is isomorphic to the ten-dimensional Jordan
2662 spin algebra, but just to be sure, we recompute its rank::
2664 sage: J = OctonionHermitianEJA(2, # long time
2665 ....: field=QQ, # long time
2666 ....: orthonormalize=False) # long time
2667 sage: J.rank.clear_cache() # long time
2668 sage: J.rank() # long time
2673 def _max_random_instance_size():
2675 The maximum rank of a random QuaternionHermitianEJA.
2677 return 1 # Dimension 1
2680 def random_instance(cls
, **kwargs
):
2682 Return a random instance of this type of algebra.
2684 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2685 return cls(n
, **kwargs
)
2687 def __init__(self
, n
, field
=AA
, **kwargs
):
2689 # Otherwise we don't get an EJA.
2690 raise ValueError("n cannot exceed 3")
2692 # We know this is a valid EJA, but will double-check
2693 # if the user passes check_axioms=True.
2694 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2696 super().__init
__(self
._denormalized
_basis
(n
,field
),
2697 self
.jordan_product
,
2698 self
.trace_inner_product
,
2702 # TODO: this could be factored out somehow, but is left here
2703 # because the MatrixEJA is not presently a subclass of the
2704 # FDEJA class that defines rank() and one().
2705 self
.rank
.set_cache(n
)
2706 idV
= self
.matrix_space().one()
2707 self
.one
.set_cache(self(idV
))
2711 def _denormalized_basis(cls
, n
, field
):
2713 Returns a basis for the space of octonion Hermitian n-by-n
2718 sage: from mjo.eja.eja_algebra import OctonionHermitianEJA
2722 sage: B = OctonionHermitianEJA._denormalized_basis(3,QQ)
2723 sage: all( M.is_hermitian() for M in B )
2729 from mjo
.octonions
import OctonionMatrixAlgebra
2730 MS
= OctonionMatrixAlgebra(n
, scalars
=field
)
2731 es
= MS
.entry_algebra().gens()
2735 for j
in range(i
+1):
2737 E_ii
= MS
.monomial( (i
,j
,es
[0]) )
2741 E_ij
= MS
.monomial( (i
,j
,e
) )
2743 # If the conjugate has a negative sign in front
2744 # of it, (j,i,ec) won't be a monomial!
2745 if (j
,i
,ec
) in MS
.indices():
2746 E_ij
+= MS
.monomial( (j
,i
,ec
) )
2748 E_ij
-= MS
.monomial( (j
,i
,-ec
) )
2751 return tuple( basis
)
2754 def trace_inner_product(X
,Y
):
2756 The octonions don't know that the reals are embedded in them,
2757 so we have to take the e0 component ourselves.
2761 sage: from mjo.eja.eja_algebra import OctonionHermitianEJA
2765 sage: J = OctonionHermitianEJA(2,field=QQ,orthonormalize=False)
2766 sage: I = J.one().to_matrix()
2767 sage: J.trace_inner_product(I, -I)
2771 return (X
*Y
).trace().real().coefficient(0)
2773 class HadamardEJA(ConcreteEJA
, RationalBasisEJA
):
2775 Return the Euclidean Jordan Algebra corresponding to the set
2776 `R^n` under the Hadamard product.
2778 Note: this is nothing more than the Cartesian product of ``n``
2779 copies of the spin algebra. Once Cartesian product algebras
2780 are implemented, this can go.
2784 sage: from mjo.eja.eja_algebra import HadamardEJA
2788 This multiplication table can be verified by hand::
2790 sage: J = HadamardEJA(3)
2791 sage: b0,b1,b2 = J.gens()
2807 We can change the generator prefix::
2809 sage: HadamardEJA(3, prefix='r').gens()
2813 def __init__(self
, n
, field
=AA
, **kwargs
):
2815 jordan_product
= lambda x
,y
: x
2816 inner_product
= lambda x
,y
: x
2818 def jordan_product(x
,y
):
2820 return P( xi
*yi
for (xi
,yi
) in zip(x
,y
) )
2822 def inner_product(x
,y
):
2825 # New defaults for keyword arguments. Don't orthonormalize
2826 # because our basis is already orthonormal with respect to our
2827 # inner-product. Don't check the axioms, because we know this
2828 # is a valid EJA... but do double-check if the user passes
2829 # check_axioms=True. Note: we DON'T override the "check_field"
2830 # default here, because the user can pass in a field!
2831 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2832 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2834 column_basis
= tuple( b
.column()
2835 for b
in FreeModule(field
, n
).basis() )
2836 super().__init
__(column_basis
,
2842 self
.rank
.set_cache(n
)
2845 self
.one
.set_cache( self
.zero() )
2847 self
.one
.set_cache( sum(self
.gens()) )
2850 def _max_random_instance_size():
2852 The maximum dimension of a random HadamardEJA.
2857 def random_instance(cls
, **kwargs
):
2859 Return a random instance of this type of algebra.
2861 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2862 return cls(n
, **kwargs
)
2865 class BilinearFormEJA(ConcreteEJA
, RationalBasisEJA
):
2867 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2868 with the half-trace inner product and jordan product ``x*y =
2869 (<Bx,y>,y_bar>, x0*y_bar + y0*x_bar)`` where `B = 1 \times B22` is
2870 a symmetric positive-definite "bilinear form" matrix. Its
2871 dimension is the size of `B`, and it has rank two in dimensions
2872 larger than two. It reduces to the ``JordanSpinEJA`` when `B` is
2873 the identity matrix of order ``n``.
2875 We insist that the one-by-one upper-left identity block of `B` be
2876 passed in as well so that we can be passed a matrix of size zero
2877 to construct a trivial algebra.
2881 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
2882 ....: JordanSpinEJA)
2886 When no bilinear form is specified, the identity matrix is used,
2887 and the resulting algebra is the Jordan spin algebra::
2889 sage: B = matrix.identity(AA,3)
2890 sage: J0 = BilinearFormEJA(B)
2891 sage: J1 = JordanSpinEJA(3)
2892 sage: J0.multiplication_table() == J0.multiplication_table()
2895 An error is raised if the matrix `B` does not correspond to a
2896 positive-definite bilinear form::
2898 sage: B = matrix.random(QQ,2,3)
2899 sage: J = BilinearFormEJA(B)
2900 Traceback (most recent call last):
2902 ValueError: bilinear form is not positive-definite
2903 sage: B = matrix.zero(QQ,3)
2904 sage: J = BilinearFormEJA(B)
2905 Traceback (most recent call last):
2907 ValueError: bilinear form is not positive-definite
2911 We can create a zero-dimensional algebra::
2913 sage: B = matrix.identity(AA,0)
2914 sage: J = BilinearFormEJA(B)
2918 We can check the multiplication condition given in the Jordan, von
2919 Neumann, and Wigner paper (and also discussed on my "On the
2920 symmetry..." paper). Note that this relies heavily on the standard
2921 choice of basis, as does anything utilizing the bilinear form
2922 matrix. We opt not to orthonormalize the basis, because if we
2923 did, we would have to normalize the `s_{i}` in a similar manner::
2925 sage: set_random_seed()
2926 sage: n = ZZ.random_element(5)
2927 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2928 sage: B11 = matrix.identity(QQ,1)
2929 sage: B22 = M.transpose()*M
2930 sage: B = block_matrix(2,2,[ [B11,0 ],
2932 sage: J = BilinearFormEJA(B, orthonormalize=False)
2933 sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
2934 sage: V = J.vector_space()
2935 sage: sis = [ J( V([0] + (M.inverse()*ei).list()).column() )
2936 ....: for ei in eis ]
2937 sage: actual = [ sis[i]*sis[j]
2938 ....: for i in range(n-1)
2939 ....: for j in range(n-1) ]
2940 sage: expected = [ J.one() if i == j else J.zero()
2941 ....: for i in range(n-1)
2942 ....: for j in range(n-1) ]
2943 sage: actual == expected
2947 def __init__(self
, B
, field
=AA
, **kwargs
):
2948 # The matrix "B" is supplied by the user in most cases,
2949 # so it makes sense to check whether or not its positive-
2950 # definite unless we are specifically asked not to...
2951 if ("check_axioms" not in kwargs
) or kwargs
["check_axioms"]:
2952 if not B
.is_positive_definite():
2953 raise ValueError("bilinear form is not positive-definite")
2955 # However, all of the other data for this EJA is computed
2956 # by us in manner that guarantees the axioms are
2957 # satisfied. So, again, unless we are specifically asked to
2958 # verify things, we'll skip the rest of the checks.
2959 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2961 def inner_product(x
,y
):
2962 return (y
.T
*B
*x
)[0,0]
2964 def jordan_product(x
,y
):
2970 z0
= inner_product(y
,x
)
2971 zbar
= y0
*xbar
+ x0
*ybar
2972 return P([z0
] + zbar
.list())
2975 column_basis
= tuple( b
.column()
2976 for b
in FreeModule(field
, n
).basis() )
2978 # TODO: I haven't actually checked this, but it seems legit.
2983 super().__init
__(column_basis
,
2987 associative
=associative
,
2990 # The rank of this algebra is two, unless we're in a
2991 # one-dimensional ambient space (because the rank is bounded
2992 # by the ambient dimension).
2993 self
.rank
.set_cache(min(n
,2))
2996 self
.one
.set_cache( self
.zero() )
2998 self
.one
.set_cache( self
.monomial(0) )
3001 def _max_random_instance_size():
3003 The maximum dimension of a random BilinearFormEJA.
3008 def random_instance(cls
, **kwargs
):
3010 Return a random instance of this algebra.
3012 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
3014 B
= matrix
.identity(ZZ
, n
)
3015 return cls(B
, **kwargs
)
3017 B11
= matrix
.identity(ZZ
, 1)
3018 M
= matrix
.random(ZZ
, n
-1)
3019 I
= matrix
.identity(ZZ
, n
-1)
3021 while alpha
.is_zero():
3022 alpha
= ZZ
.random_element().abs()
3023 B22
= M
.transpose()*M
+ alpha
*I
3025 from sage
.matrix
.special
import block_matrix
3026 B
= block_matrix(2,2, [ [B11
, ZZ(0) ],
3029 return cls(B
, **kwargs
)
3032 class JordanSpinEJA(BilinearFormEJA
):
3034 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
3035 with the usual inner product and jordan product ``x*y =
3036 (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
3041 sage: from mjo.eja.eja_algebra import JordanSpinEJA
3045 This multiplication table can be verified by hand::
3047 sage: J = JordanSpinEJA(4)
3048 sage: b0,b1,b2,b3 = J.gens()
3064 We can change the generator prefix::
3066 sage: JordanSpinEJA(2, prefix='B').gens()
3071 Ensure that we have the usual inner product on `R^n`::
3073 sage: set_random_seed()
3074 sage: J = JordanSpinEJA.random_instance()
3075 sage: x,y = J.random_elements(2)
3076 sage: actual = x.inner_product(y)
3077 sage: expected = x.to_vector().inner_product(y.to_vector())
3078 sage: actual == expected
3082 def __init__(self
, n
, *args
, **kwargs
):
3083 # This is a special case of the BilinearFormEJA with the
3084 # identity matrix as its bilinear form.
3085 B
= matrix
.identity(ZZ
, n
)
3087 # Don't orthonormalize because our basis is already
3088 # orthonormal with respect to our inner-product.
3089 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
3091 # But also don't pass check_field=False here, because the user
3092 # can pass in a field!
3093 super().__init
__(B
, *args
, **kwargs
)
3096 def _max_random_instance_size():
3098 The maximum dimension of a random JordanSpinEJA.
3103 def random_instance(cls
, **kwargs
):
3105 Return a random instance of this type of algebra.
3107 Needed here to override the implementation for ``BilinearFormEJA``.
3109 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
3110 return cls(n
, **kwargs
)
3113 class TrivialEJA(ConcreteEJA
, RationalBasisEJA
):
3115 The trivial Euclidean Jordan algebra consisting of only a zero element.
3119 sage: from mjo.eja.eja_algebra import TrivialEJA
3123 sage: J = TrivialEJA()
3130 sage: 7*J.one()*12*J.one()
3132 sage: J.one().inner_product(J.one())
3134 sage: J.one().norm()
3136 sage: J.one().subalgebra_generated_by()
3137 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
3142 def __init__(self
, **kwargs
):
3143 jordan_product
= lambda x
,y
: x
3144 inner_product
= lambda x
,y
: 0
3147 # New defaults for keyword arguments
3148 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
3149 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
3151 super().__init
__(basis
,
3157 # The rank is zero using my definition, namely the dimension of the
3158 # largest subalgebra generated by any element.
3159 self
.rank
.set_cache(0)
3160 self
.one
.set_cache( self
.zero() )
3163 def random_instance(cls
, **kwargs
):
3164 # We don't take a "size" argument so the superclass method is
3165 # inappropriate for us.
3166 return cls(**kwargs
)
3169 class CartesianProductEJA(FiniteDimensionalEJA
):
3171 The external (orthogonal) direct sum of two or more Euclidean
3172 Jordan algebras. Every Euclidean Jordan algebra decomposes into an
3173 orthogonal direct sum of simple Euclidean Jordan algebras which is
3174 then isometric to a Cartesian product, so no generality is lost by
3175 providing only this construction.
3179 sage: from mjo.eja.eja_algebra import (random_eja,
3180 ....: CartesianProductEJA,
3182 ....: JordanSpinEJA,
3183 ....: RealSymmetricEJA)
3187 The Jordan product is inherited from our factors and implemented by
3188 our CombinatorialFreeModule Cartesian product superclass::
3190 sage: set_random_seed()
3191 sage: J1 = HadamardEJA(2)
3192 sage: J2 = RealSymmetricEJA(2)
3193 sage: J = cartesian_product([J1,J2])
3194 sage: x,y = J.random_elements(2)
3198 The ability to retrieve the original factors is implemented by our
3199 CombinatorialFreeModule Cartesian product superclass::
3201 sage: J1 = HadamardEJA(2, field=QQ)
3202 sage: J2 = JordanSpinEJA(3, field=QQ)
3203 sage: J = cartesian_product([J1,J2])
3204 sage: J.cartesian_factors()
3205 (Euclidean Jordan algebra of dimension 2 over Rational Field,
3206 Euclidean Jordan algebra of dimension 3 over Rational Field)
3208 You can provide more than two factors::
3210 sage: J1 = HadamardEJA(2)
3211 sage: J2 = JordanSpinEJA(3)
3212 sage: J3 = RealSymmetricEJA(3)
3213 sage: cartesian_product([J1,J2,J3])
3214 Euclidean Jordan algebra of dimension 2 over Algebraic Real
3215 Field (+) Euclidean Jordan algebra of dimension 3 over Algebraic
3216 Real Field (+) Euclidean Jordan algebra of dimension 6 over
3217 Algebraic Real Field
3219 Rank is additive on a Cartesian product::
3221 sage: J1 = HadamardEJA(1)
3222 sage: J2 = RealSymmetricEJA(2)
3223 sage: J = cartesian_product([J1,J2])
3224 sage: J1.rank.clear_cache()
3225 sage: J2.rank.clear_cache()
3226 sage: J.rank.clear_cache()
3229 sage: J.rank() == J1.rank() + J2.rank()
3232 The same rank computation works over the rationals, with whatever
3235 sage: J1 = HadamardEJA(1, field=QQ, orthonormalize=False)
3236 sage: J2 = RealSymmetricEJA(2, field=QQ, orthonormalize=False)
3237 sage: J = cartesian_product([J1,J2])
3238 sage: J1.rank.clear_cache()
3239 sage: J2.rank.clear_cache()
3240 sage: J.rank.clear_cache()
3243 sage: J.rank() == J1.rank() + J2.rank()
3246 The product algebra will be associative if and only if all of its
3247 components are associative::
3249 sage: J1 = HadamardEJA(2)
3250 sage: J1.is_associative()
3252 sage: J2 = HadamardEJA(3)
3253 sage: J2.is_associative()
3255 sage: J3 = RealSymmetricEJA(3)
3256 sage: J3.is_associative()
3258 sage: CP1 = cartesian_product([J1,J2])
3259 sage: CP1.is_associative()
3261 sage: CP2 = cartesian_product([J1,J3])
3262 sage: CP2.is_associative()
3265 Cartesian products of Cartesian products work::
3267 sage: J1 = JordanSpinEJA(1)
3268 sage: J2 = JordanSpinEJA(1)
3269 sage: J3 = JordanSpinEJA(1)
3270 sage: J = cartesian_product([J1,cartesian_product([J2,J3])])
3271 sage: J.multiplication_table()
3272 +----++----+----+----+
3273 | * || b0 | b1 | b2 |
3274 +====++====+====+====+
3275 | b0 || b0 | 0 | 0 |
3276 +----++----+----+----+
3277 | b1 || 0 | b1 | 0 |
3278 +----++----+----+----+
3279 | b2 || 0 | 0 | b2 |
3280 +----++----+----+----+
3281 sage: HadamardEJA(3).multiplication_table()
3282 +----++----+----+----+
3283 | * || b0 | b1 | b2 |
3284 +====++====+====+====+
3285 | b0 || b0 | 0 | 0 |
3286 +----++----+----+----+
3287 | b1 || 0 | b1 | 0 |
3288 +----++----+----+----+
3289 | b2 || 0 | 0 | b2 |
3290 +----++----+----+----+
3294 All factors must share the same base field::
3296 sage: J1 = HadamardEJA(2, field=QQ)
3297 sage: J2 = RealSymmetricEJA(2)
3298 sage: CartesianProductEJA((J1,J2))
3299 Traceback (most recent call last):
3301 ValueError: all factors must share the same base field
3303 The cached unit element is the same one that would be computed::
3305 sage: set_random_seed() # long time
3306 sage: J1 = random_eja() # long time
3307 sage: J2 = random_eja() # long time
3308 sage: J = cartesian_product([J1,J2]) # long time
3309 sage: actual = J.one() # long time
3310 sage: J.one.clear_cache() # long time
3311 sage: expected = J.one() # long time
3312 sage: actual == expected # long time
3316 Element
= FiniteDimensionalEJAElement
3319 def __init__(self
, factors
, **kwargs
):
3324 self
._sets
= factors
3326 field
= factors
[0].base_ring()
3327 if not all( J
.base_ring() == field
for J
in factors
):
3328 raise ValueError("all factors must share the same base field")
3330 associative
= all( f
.is_associative() for f
in factors
)
3332 MS
= self
.matrix_space()
3336 for b
in factors
[i
].matrix_basis():
3341 basis
= tuple( MS(b
) for b
in basis
)
3343 # Define jordan/inner products that operate on that matrix_basis.
3344 def jordan_product(x
,y
):
3346 (factors
[i
](x
[i
])*factors
[i
](y
[i
])).to_matrix()
3350 def inner_product(x
, y
):
3352 factors
[i
](x
[i
]).inner_product(factors
[i
](y
[i
]))
3356 # There's no need to check the field since it already came
3357 # from an EJA. Likewise the axioms are guaranteed to be
3358 # satisfied, unless the guy writing this class sucks.
3360 # If you want the basis to be orthonormalized, orthonormalize
3362 FiniteDimensionalEJA
.__init
__(self
,
3367 orthonormalize
=False,
3368 associative
=associative
,
3369 cartesian_product
=True,
3373 ones
= tuple(J
.one().to_matrix() for J
in factors
)
3374 self
.one
.set_cache(self(ones
))
3375 self
.rank
.set_cache(sum(J
.rank() for J
in factors
))
3377 def cartesian_factors(self
):
3378 # Copy/pasted from CombinatorialFreeModule_CartesianProduct.
3381 def cartesian_factor(self
, i
):
3383 Return the ``i``th factor of this algebra.
3385 return self
._sets
[i
]
3388 # Copy/pasted from CombinatorialFreeModule_CartesianProduct.
3389 from sage
.categories
.cartesian_product
import cartesian_product
3390 return cartesian_product
.symbol
.join("%s" % factor
3391 for factor
in self
._sets
)
3393 def matrix_space(self
):
3395 Return the space that our matrix basis lives in as a Cartesian
3398 We don't simply use the ``cartesian_product()`` functor here
3399 because it acts differently on SageMath MatrixSpaces and our
3400 custom MatrixAlgebras, which are CombinatorialFreeModules. We
3401 always want the result to be represented (and indexed) as
3406 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
3408 ....: OctonionHermitianEJA,
3409 ....: RealSymmetricEJA)
3413 sage: J1 = HadamardEJA(1)
3414 sage: J2 = RealSymmetricEJA(2)
3415 sage: J = cartesian_product([J1,J2])
3416 sage: J.matrix_space()
3417 The Cartesian product of (Full MatrixSpace of 1 by 1 dense
3418 matrices over Algebraic Real Field, Full MatrixSpace of 2
3419 by 2 dense matrices over Algebraic Real Field)
3423 sage: J1 = ComplexHermitianEJA(1)
3424 sage: J2 = ComplexHermitianEJA(1)
3425 sage: J = cartesian_product([J1,J2])
3426 sage: J.one().to_matrix()[0]
3429 sage: J.one().to_matrix()[1]
3435 sage: J1 = OctonionHermitianEJA(1)
3436 sage: J2 = OctonionHermitianEJA(1)
3437 sage: J = cartesian_product([J1,J2])
3438 sage: J.one().to_matrix()[0]
3442 sage: J.one().to_matrix()[1]
3448 scalars
= self
.cartesian_factor(0).base_ring()
3450 # This category isn't perfect, but is good enough for what we
3452 cat
= MagmaticAlgebras(scalars
).FiniteDimensional().WithBasis()
3453 cat
= cat
.Unital().CartesianProducts()
3454 factors
= tuple( J
.matrix_space() for J
in self
.cartesian_factors() )
3456 from sage
.sets
.cartesian_product
import CartesianProduct
3457 return CartesianProduct(factors
, cat
)
3461 def cartesian_projection(self
, i
):
3465 sage: from mjo.eja.eja_algebra import (random_eja,
3466 ....: JordanSpinEJA,
3468 ....: RealSymmetricEJA,
3469 ....: ComplexHermitianEJA)
3473 The projection morphisms are Euclidean Jordan algebra
3476 sage: J1 = HadamardEJA(2)
3477 sage: J2 = RealSymmetricEJA(2)
3478 sage: J = cartesian_product([J1,J2])
3479 sage: J.cartesian_projection(0)
3480 Linear operator between finite-dimensional Euclidean Jordan
3481 algebras represented by the matrix:
3484 Domain: Euclidean Jordan algebra of dimension 2 over Algebraic
3485 Real Field (+) Euclidean Jordan algebra of dimension 3 over
3486 Algebraic Real Field
3487 Codomain: Euclidean Jordan algebra of dimension 2 over Algebraic
3489 sage: J.cartesian_projection(1)
3490 Linear operator between finite-dimensional Euclidean Jordan
3491 algebras represented by the matrix:
3495 Domain: Euclidean Jordan algebra of dimension 2 over Algebraic
3496 Real Field (+) Euclidean Jordan algebra of dimension 3 over
3497 Algebraic Real Field
3498 Codomain: Euclidean Jordan algebra of dimension 3 over Algebraic
3501 The projections work the way you'd expect on the vector
3502 representation of an element::
3504 sage: J1 = JordanSpinEJA(2)
3505 sage: J2 = ComplexHermitianEJA(2)
3506 sage: J = cartesian_product([J1,J2])
3507 sage: pi_left = J.cartesian_projection(0)
3508 sage: pi_right = J.cartesian_projection(1)
3509 sage: pi_left(J.one()).to_vector()
3511 sage: pi_right(J.one()).to_vector()
3513 sage: J.one().to_vector()
3518 The answer never changes::
3520 sage: set_random_seed()
3521 sage: J1 = random_eja()
3522 sage: J2 = random_eja()
3523 sage: J = cartesian_product([J1,J2])
3524 sage: P0 = J.cartesian_projection(0)
3525 sage: P1 = J.cartesian_projection(0)
3530 offset
= sum( self
.cartesian_factor(k
).dimension()
3532 Ji
= self
.cartesian_factor(i
)
3533 Pi
= self
._module
_morphism
(lambda j
: Ji
.monomial(j
- offset
),
3536 return FiniteDimensionalEJAOperator(self
,Ji
,Pi
.matrix())
3539 def cartesian_embedding(self
, i
):
3543 sage: from mjo.eja.eja_algebra import (random_eja,
3544 ....: JordanSpinEJA,
3546 ....: RealSymmetricEJA)
3550 The embedding morphisms are Euclidean Jordan algebra
3553 sage: J1 = HadamardEJA(2)
3554 sage: J2 = RealSymmetricEJA(2)
3555 sage: J = cartesian_product([J1,J2])
3556 sage: J.cartesian_embedding(0)
3557 Linear operator between finite-dimensional Euclidean Jordan
3558 algebras represented by the matrix:
3564 Domain: Euclidean Jordan algebra of dimension 2 over
3565 Algebraic Real Field
3566 Codomain: Euclidean Jordan algebra of dimension 2 over
3567 Algebraic Real Field (+) Euclidean Jordan algebra of
3568 dimension 3 over Algebraic Real Field
3569 sage: J.cartesian_embedding(1)
3570 Linear operator between finite-dimensional Euclidean Jordan
3571 algebras represented by the matrix:
3577 Domain: Euclidean Jordan algebra of dimension 3 over
3578 Algebraic Real Field
3579 Codomain: Euclidean Jordan algebra of dimension 2 over
3580 Algebraic Real Field (+) Euclidean Jordan algebra of
3581 dimension 3 over Algebraic Real Field
3583 The embeddings work the way you'd expect on the vector
3584 representation of an element::
3586 sage: J1 = JordanSpinEJA(3)
3587 sage: J2 = RealSymmetricEJA(2)
3588 sage: J = cartesian_product([J1,J2])
3589 sage: iota_left = J.cartesian_embedding(0)
3590 sage: iota_right = J.cartesian_embedding(1)
3591 sage: iota_left(J1.zero()) == J.zero()
3593 sage: iota_right(J2.zero()) == J.zero()
3595 sage: J1.one().to_vector()
3597 sage: iota_left(J1.one()).to_vector()
3599 sage: J2.one().to_vector()
3601 sage: iota_right(J2.one()).to_vector()
3603 sage: J.one().to_vector()
3608 The answer never changes::
3610 sage: set_random_seed()
3611 sage: J1 = random_eja()
3612 sage: J2 = random_eja()
3613 sage: J = cartesian_product([J1,J2])
3614 sage: E0 = J.cartesian_embedding(0)
3615 sage: E1 = J.cartesian_embedding(0)
3619 Composing a projection with the corresponding inclusion should
3620 produce the identity map, and mismatching them should produce
3623 sage: set_random_seed()
3624 sage: J1 = random_eja()
3625 sage: J2 = random_eja()
3626 sage: J = cartesian_product([J1,J2])
3627 sage: iota_left = J.cartesian_embedding(0)
3628 sage: iota_right = J.cartesian_embedding(1)
3629 sage: pi_left = J.cartesian_projection(0)
3630 sage: pi_right = J.cartesian_projection(1)
3631 sage: pi_left*iota_left == J1.one().operator()
3633 sage: pi_right*iota_right == J2.one().operator()
3635 sage: (pi_left*iota_right).is_zero()
3637 sage: (pi_right*iota_left).is_zero()
3641 offset
= sum( self
.cartesian_factor(k
).dimension()
3643 Ji
= self
.cartesian_factor(i
)
3644 Ei
= Ji
._module
_morphism
(lambda j
: self
.monomial(j
+ offset
),
3646 return FiniteDimensionalEJAOperator(Ji
,self
,Ei
.matrix())
3650 FiniteDimensionalEJA
.CartesianProduct
= CartesianProductEJA
3652 class RationalBasisCartesianProductEJA(CartesianProductEJA
,
3655 A separate class for products of algebras for which we know a
3660 sage: from mjo.eja.eja_algebra import (HadamardEJA,
3661 ....: JordanSpinEJA,
3662 ....: OctonionHermitianEJA,
3663 ....: RealSymmetricEJA)
3667 This gives us fast characteristic polynomial computations in
3668 product algebras, too::
3671 sage: J1 = JordanSpinEJA(2)
3672 sage: J2 = RealSymmetricEJA(3)
3673 sage: J = cartesian_product([J1,J2])
3674 sage: J.characteristic_polynomial_of().degree()
3681 The ``cartesian_product()`` function only uses the first factor to
3682 decide where the result will live; thus we have to be careful to
3683 check that all factors do indeed have a `_rational_algebra` member
3684 before we try to access it::
3686 sage: J1 = OctonionHermitianEJA(1) # no rational basis
3687 sage: J2 = HadamardEJA(2)
3688 sage: cartesian_product([J1,J2])
3689 Euclidean Jordan algebra of dimension 1 over Algebraic Real Field
3690 (+) Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
3691 sage: cartesian_product([J2,J1])
3692 Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
3693 (+) Euclidean Jordan algebra of dimension 1 over Algebraic Real Field
3696 def __init__(self
, algebras
, **kwargs
):
3697 CartesianProductEJA
.__init
__(self
, algebras
, **kwargs
)
3699 self
._rational
_algebra
= None
3700 if self
.vector_space().base_field() is not QQ
:
3701 if all( hasattr(r
, "_rational_algebra") for r
in algebras
):
3702 self
._rational
_algebra
= cartesian_product([
3703 r
._rational
_algebra
for r
in algebras
3707 RationalBasisEJA
.CartesianProduct
= RationalBasisCartesianProductEJA
3709 def random_eja(*args
, **kwargs
):
3710 J1
= ConcreteEJA
.random_instance(*args
, **kwargs
)
3712 # This might make Cartesian products appear roughly as often as
3713 # any other ConcreteEJA.
3714 if ZZ
.random_element(len(ConcreteEJA
.__subclasses
__()) + 1) == 0:
3715 # Use random_eja() again so we can get more than two factors.
3716 J2
= random_eja(*args
, **kwargs
)
3717 J
= cartesian_product([J1
,J2
])