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`
36 Missing from this list is the algebra of three-by-three octononion
37 Hermitian matrices, as there is (as of yet) no implementation of the
38 octonions in SageMath. In addition to these, we provide two other
39 example constructions,
41 * :class:`HadamardEJA`
44 The Jordan spin algebra is a bilinear form algebra where the bilinear
45 form is the identity. The Hadamard EJA is simply a Cartesian product
46 of one-dimensional spin algebras. And last but not least, the trivial
47 EJA is exactly what you think. Cartesian products of these are also
48 supported using the usual ``cartesian_product()`` function; as a
49 result, we support (up to isomorphism) all Euclidean Jordan algebras
50 that don't involve octonions.
54 sage: from mjo.eja.eja_algebra import random_eja
59 Euclidean Jordan algebra of dimension...
62 from itertools
import repeat
64 from sage
.algebras
.quatalg
.quaternion_algebra
import QuaternionAlgebra
65 from sage
.categories
.magmatic_algebras
import MagmaticAlgebras
66 from sage
.categories
.sets_cat
import cartesian_product
67 from sage
.combinat
.free_module
import CombinatorialFreeModule
68 from sage
.matrix
.constructor
import matrix
69 from sage
.matrix
.matrix_space
import MatrixSpace
70 from sage
.misc
.cachefunc
import cached_method
71 from sage
.misc
.table
import table
72 from sage
.modules
.free_module
import FreeModule
, VectorSpace
73 from sage
.rings
.all
import (ZZ
, QQ
, AA
, QQbar
, RR
, RLF
, CLF
,
76 from mjo
.eja
.eja_element
import FiniteDimensionalEJAElement
77 from mjo
.eja
.eja_operator
import FiniteDimensionalEJAOperator
78 from mjo
.eja
.eja_utils
import _all2list
, _mat2vec
80 class FiniteDimensionalEJA(CombinatorialFreeModule
):
82 A finite-dimensional Euclidean Jordan algebra.
86 - ``basis`` -- a tuple; a tuple of basis elements in "matrix
87 form," which must be the same form as the arguments to
88 ``jordan_product`` and ``inner_product``. In reality, "matrix
89 form" can be either vectors, matrices, or a Cartesian product
90 (ordered tuple) of vectors or matrices. All of these would
91 ideally be vector spaces in sage with no special-casing
92 needed; but in reality we turn vectors into column-matrices
93 and Cartesian products `(a,b)` into column matrices
94 `(a,b)^{T}` after converting `a` and `b` themselves.
96 - ``jordan_product`` -- a function; afunction of two ``basis``
97 elements (in matrix form) that returns their jordan product,
98 also in matrix form; this will be applied to ``basis`` to
99 compute a multiplication table for the algebra.
101 - ``inner_product`` -- a function; a function of two ``basis``
102 elements (in matrix form) that returns their inner
103 product. This will be applied to ``basis`` to compute an
104 inner-product table (basically a matrix) for this algebra.
106 - ``field`` -- a subfield of the reals (default: ``AA``); the scalar
107 field for the algebra.
109 - ``orthonormalize`` -- boolean (default: ``True``); whether or
110 not to orthonormalize the basis. Doing so is expensive and
111 generally rules out using the rationals as your ``field``, but
112 is required for spectral decompositions.
116 sage: from mjo.eja.eja_algebra import random_eja
120 We should compute that an element subalgebra is associative even
121 if we circumvent the element method::
123 sage: set_random_seed()
124 sage: J = random_eja(field=QQ,orthonormalize=False)
125 sage: x = J.random_element()
126 sage: A = x.subalgebra_generated_by(orthonormalize=False)
127 sage: basis = tuple(b.superalgebra_element() for b in A.basis())
128 sage: J.subalgebra(basis, orthonormalize=False).is_associative()
132 Element
= FiniteDimensionalEJAElement
141 cartesian_product
=False,
149 if not field
.is_subring(RR
):
150 # Note: this does return true for the real algebraic
151 # field, the rationals, and any quadratic field where
152 # we've specified a real embedding.
153 raise ValueError("scalar field is not real")
155 from mjo
.eja
.eja_utils
import _change_ring
156 # If the basis given to us wasn't over the field that it's
157 # supposed to be over, fix that. Or, you know, crash.
158 basis
= tuple( _change_ring(b
, field
) for b
in basis
)
161 # Check commutativity of the Jordan and inner-products.
162 # This has to be done before we build the multiplication
163 # and inner-product tables/matrices, because we take
164 # advantage of symmetry in the process.
165 if not all( jordan_product(bi
,bj
) == jordan_product(bj
,bi
)
168 raise ValueError("Jordan product is not commutative")
170 if not all( inner_product(bi
,bj
) == inner_product(bj
,bi
)
173 raise ValueError("inner-product is not commutative")
176 category
= MagmaticAlgebras(field
).FiniteDimensional()
177 category
= category
.WithBasis().Unital().Commutative()
179 if associative
is None:
180 # We should figure it out. As with check_axioms, we have to do
181 # this without the help of the _jordan_product_is_associative()
182 # method because we need to know the category before we
183 # initialize the algebra.
184 associative
= all( jordan_product(jordan_product(bi
,bj
),bk
)
186 jordan_product(bi
,jordan_product(bj
,bk
))
192 # Element subalgebras can take advantage of this.
193 category
= category
.Associative()
194 if cartesian_product
:
195 # Use join() here because otherwise we only get the
196 # "Cartesian product of..." and not the things themselves.
197 category
= category
.join([category
,
198 category
.CartesianProducts()])
200 # Call the superclass constructor so that we can use its from_vector()
201 # method to build our multiplication table.
202 CombinatorialFreeModule
.__init
__(self
,
209 # Now comes all of the hard work. We'll be constructing an
210 # ambient vector space V that our (vectorized) basis lives in,
211 # as well as a subspace W of V spanned by those (vectorized)
212 # basis elements. The W-coordinates are the coefficients that
213 # we see in things like x = 1*e1 + 2*e2.
218 degree
= len(_all2list(basis
[0]))
220 # Build an ambient space that fits our matrix basis when
221 # written out as "long vectors."
222 V
= VectorSpace(field
, degree
)
224 # The matrix that will hole the orthonormal -> unorthonormal
225 # coordinate transformation.
226 self
._deortho
_matrix
= None
229 # Save a copy of the un-orthonormalized basis for later.
230 # Convert it to ambient V (vector) coordinates while we're
231 # at it, because we'd have to do it later anyway.
232 deortho_vector_basis
= tuple( V(_all2list(b
)) for b
in basis
)
234 from mjo
.eja
.eja_utils
import gram_schmidt
235 basis
= tuple(gram_schmidt(basis
, inner_product
))
237 # Save the (possibly orthonormalized) matrix basis for
239 self
._matrix
_basis
= basis
241 # Now create the vector space for the algebra, which will have
242 # its own set of non-ambient coordinates (in terms of the
244 vector_basis
= tuple( V(_all2list(b
)) for b
in basis
)
245 W
= V
.span_of_basis( vector_basis
, check
=check_axioms
)
248 # Now "W" is the vector space of our algebra coordinates. The
249 # variables "X1", "X2",... refer to the entries of vectors in
250 # W. Thus to convert back and forth between the orthonormal
251 # coordinates and the given ones, we need to stick the original
253 U
= V
.span_of_basis( deortho_vector_basis
, check
=check_axioms
)
254 self
._deortho
_matrix
= matrix( U
.coordinate_vector(q
)
255 for q
in vector_basis
)
258 # Now we actually compute the multiplication and inner-product
259 # tables/matrices using the possibly-orthonormalized basis.
260 self
._inner
_product
_matrix
= matrix
.identity(field
, n
)
261 self
._multiplication
_table
= [ [0 for j
in range(i
+1)]
264 # Note: the Jordan and inner-products are defined in terms
265 # of the ambient basis. It's important that their arguments
266 # are in ambient coordinates as well.
269 # ortho basis w.r.t. ambient coords
273 # The jordan product returns a matrixy answer, so we
274 # have to convert it to the algebra coordinates.
275 elt
= jordan_product(q_i
, q_j
)
276 elt
= W
.coordinate_vector(V(_all2list(elt
)))
277 self
._multiplication
_table
[i
][j
] = self
.from_vector(elt
)
279 if not orthonormalize
:
280 # If we're orthonormalizing the basis with respect
281 # to an inner-product, then the inner-product
282 # matrix with respect to the resulting basis is
283 # just going to be the identity.
284 ip
= inner_product(q_i
, q_j
)
285 self
._inner
_product
_matrix
[i
,j
] = ip
286 self
._inner
_product
_matrix
[j
,i
] = ip
288 self
._inner
_product
_matrix
._cache
= {'hermitian': True}
289 self
._inner
_product
_matrix
.set_immutable()
292 if not self
._is
_jordanian
():
293 raise ValueError("Jordan identity does not hold")
294 if not self
._inner
_product
_is
_associative
():
295 raise ValueError("inner product is not associative")
298 def _coerce_map_from_base_ring(self
):
300 Disable the map from the base ring into the algebra.
302 Performing a nonsense conversion like this automatically
303 is counterpedagogical. The fallback is to try the usual
304 element constructor, which should also fail.
308 sage: from mjo.eja.eja_algebra import random_eja
312 sage: set_random_seed()
313 sage: J = random_eja()
315 Traceback (most recent call last):
317 ValueError: not an element of this algebra
323 def product_on_basis(self
, i
, j
):
325 Returns the Jordan product of the `i` and `j`th basis elements.
327 This completely defines the Jordan product on the algebra, and
328 is used direclty by our superclass machinery to implement
333 sage: from mjo.eja.eja_algebra import random_eja
337 sage: set_random_seed()
338 sage: J = random_eja()
339 sage: n = J.dimension()
342 sage: ei_ej = J.zero()*J.zero()
344 ....: i = ZZ.random_element(n)
345 ....: j = ZZ.random_element(n)
346 ....: ei = J.monomial(i)
347 ....: ej = J.monomial(j)
348 ....: ei_ej = J.product_on_basis(i,j)
353 # We only stored the lower-triangular portion of the
354 # multiplication table.
356 return self
._multiplication
_table
[i
][j
]
358 return self
._multiplication
_table
[j
][i
]
360 def inner_product(self
, x
, y
):
362 The inner product associated with this Euclidean Jordan algebra.
364 Defaults to the trace inner product, but can be overridden by
365 subclasses if they are sure that the necessary properties are
370 sage: from mjo.eja.eja_algebra import (random_eja,
372 ....: BilinearFormEJA)
376 Our inner product is "associative," which means the following for
377 a symmetric bilinear form::
379 sage: set_random_seed()
380 sage: J = random_eja()
381 sage: x,y,z = J.random_elements(3)
382 sage: (x*y).inner_product(z) == y.inner_product(x*z)
387 Ensure that this is the usual inner product for the algebras
390 sage: set_random_seed()
391 sage: J = HadamardEJA.random_instance()
392 sage: x,y = J.random_elements(2)
393 sage: actual = x.inner_product(y)
394 sage: expected = x.to_vector().inner_product(y.to_vector())
395 sage: actual == expected
398 Ensure that this is one-half of the trace inner-product in a
399 BilinearFormEJA that isn't just the reals (when ``n`` isn't
400 one). This is in Faraut and Koranyi, and also my "On the
403 sage: set_random_seed()
404 sage: J = BilinearFormEJA.random_instance()
405 sage: n = J.dimension()
406 sage: x = J.random_element()
407 sage: y = J.random_element()
408 sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
412 B
= self
._inner
_product
_matrix
413 return (B
*x
.to_vector()).inner_product(y
.to_vector())
416 def is_associative(self
):
418 Return whether or not this algebra's Jordan product is associative.
422 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
426 sage: J = ComplexHermitianEJA(3, field=QQ, orthonormalize=False)
427 sage: J.is_associative()
429 sage: x = sum(J.gens())
430 sage: A = x.subalgebra_generated_by(orthonormalize=False)
431 sage: A.is_associative()
435 return "Associative" in self
.category().axioms()
437 def _is_commutative(self
):
439 Whether or not this algebra's multiplication table is commutative.
441 This method should of course always return ``True``, unless
442 this algebra was constructed with ``check_axioms=False`` and
443 passed an invalid multiplication table.
445 return all( x
*y
== y
*x
for x
in self
.gens() for y
in self
.gens() )
447 def _is_jordanian(self
):
449 Whether or not this algebra's multiplication table respects the
450 Jordan identity `(x^{2})(xy) = x(x^{2}y)`.
452 We only check one arrangement of `x` and `y`, so for a
453 ``True`` result to be truly true, you should also check
454 :meth:`_is_commutative`. This method should of course always
455 return ``True``, unless this algebra was constructed with
456 ``check_axioms=False`` and passed an invalid multiplication table.
458 return all( (self
.monomial(i
)**2)*(self
.monomial(i
)*self
.monomial(j
))
460 (self
.monomial(i
))*((self
.monomial(i
)**2)*self
.monomial(j
))
461 for i
in range(self
.dimension())
462 for j
in range(self
.dimension()) )
464 def _jordan_product_is_associative(self
):
466 Return whether or not this algebra's Jordan product is
467 associative; that is, whether or not `x*(y*z) = (x*y)*z`
470 This method should agree with :meth:`is_associative` unless
471 you lied about the value of the ``associative`` parameter
472 when you constructed the algebra.
476 sage: from mjo.eja.eja_algebra import (random_eja,
477 ....: RealSymmetricEJA,
478 ....: ComplexHermitianEJA,
479 ....: QuaternionHermitianEJA)
483 sage: J = RealSymmetricEJA(4, orthonormalize=False)
484 sage: J._jordan_product_is_associative()
486 sage: x = sum(J.gens())
487 sage: A = x.subalgebra_generated_by()
488 sage: A._jordan_product_is_associative()
493 sage: J = ComplexHermitianEJA(2,field=QQ,orthonormalize=False)
494 sage: J._jordan_product_is_associative()
496 sage: x = sum(J.gens())
497 sage: A = x.subalgebra_generated_by(orthonormalize=False)
498 sage: A._jordan_product_is_associative()
503 sage: J = QuaternionHermitianEJA(2)
504 sage: J._jordan_product_is_associative()
506 sage: x = sum(J.gens())
507 sage: A = x.subalgebra_generated_by()
508 sage: A._jordan_product_is_associative()
513 The values we've presupplied to the constructors agree with
516 sage: set_random_seed()
517 sage: J = random_eja()
518 sage: J.is_associative() == J._jordan_product_is_associative()
524 # Used to check whether or not something is zero.
527 # I don't know of any examples that make this magnitude
528 # necessary because I don't know how to make an
529 # associative algebra when the element subalgebra
530 # construction is unreliable (as it is over RDF; we can't
531 # find the degree of an element because we can't compute
532 # the rank of a matrix). But even multiplication of floats
533 # is non-associative, so *some* epsilon is needed... let's
534 # just take the one from _inner_product_is_associative?
537 for i
in range(self
.dimension()):
538 for j
in range(self
.dimension()):
539 for k
in range(self
.dimension()):
543 diff
= (x
*y
)*z
- x
*(y
*z
)
545 if diff
.norm() > epsilon
:
550 def _inner_product_is_associative(self
):
552 Return whether or not this algebra's inner product `B` is
553 associative; that is, whether or not `B(xy,z) = B(x,yz)`.
555 This method should of course always return ``True``, unless
556 this algebra was constructed with ``check_axioms=False`` and
557 passed an invalid Jordan or inner-product.
561 # Used to check whether or not something is zero.
564 # This choice is sufficient to allow the construction of
565 # QuaternionHermitianEJA(2, field=RDF) with check_axioms=True.
568 for i
in range(self
.dimension()):
569 for j
in range(self
.dimension()):
570 for k
in range(self
.dimension()):
574 diff
= (x
*y
).inner_product(z
) - x
.inner_product(y
*z
)
576 if diff
.abs() > epsilon
:
581 def _element_constructor_(self
, elt
):
583 Construct an element of this algebra from its vector or matrix
586 This gets called only after the parent element _call_ method
587 fails to find a coercion for the argument.
591 sage: from mjo.eja.eja_algebra import (random_eja,
594 ....: RealSymmetricEJA)
598 The identity in `S^n` is converted to the identity in the EJA::
600 sage: J = RealSymmetricEJA(3)
601 sage: I = matrix.identity(QQ,3)
602 sage: J(I) == J.one()
605 This skew-symmetric matrix can't be represented in the EJA::
607 sage: J = RealSymmetricEJA(3)
608 sage: A = matrix(QQ,3, lambda i,j: i-j)
610 Traceback (most recent call last):
612 ValueError: not an element of this algebra
614 Tuples work as well, provided that the matrix basis for the
615 algebra consists of them::
617 sage: J1 = HadamardEJA(3)
618 sage: J2 = RealSymmetricEJA(2)
619 sage: J = cartesian_product([J1,J2])
620 sage: J( (J1.matrix_basis()[1], J2.matrix_basis()[2]) )
625 Ensure that we can convert any element back and forth
626 faithfully between its matrix and algebra representations::
628 sage: set_random_seed()
629 sage: J = random_eja()
630 sage: x = J.random_element()
631 sage: J(x.to_matrix()) == x
634 We cannot coerce elements between algebras just because their
635 matrix representations are compatible::
637 sage: J1 = HadamardEJA(3)
638 sage: J2 = JordanSpinEJA(3)
640 Traceback (most recent call last):
642 ValueError: not an element of this algebra
644 Traceback (most recent call last):
646 ValueError: not an element of this algebra
648 msg
= "not an element of this algebra"
649 if elt
in self
.base_ring():
650 # Ensure that no base ring -> algebra coercion is performed
651 # by this method. There's some stupidity in sage that would
652 # otherwise propagate to this method; for example, sage thinks
653 # that the integer 3 belongs to the space of 2-by-2 matrices.
654 raise ValueError(msg
)
657 # Try to convert a vector into a column-matrix...
659 except (AttributeError, TypeError):
660 # and ignore failure, because we weren't really expecting
661 # a vector as an argument anyway.
664 if elt
not in self
.matrix_space():
665 raise ValueError(msg
)
667 # Thanks for nothing! Matrix spaces aren't vector spaces in
668 # Sage, so we have to figure out its matrix-basis coordinates
669 # ourselves. We use the basis space's ring instead of the
670 # element's ring because the basis space might be an algebraic
671 # closure whereas the base ring of the 3-by-3 identity matrix
672 # could be QQ instead of QQbar.
674 # And, we also have to handle Cartesian product bases (when
675 # the matrix basis consists of tuples) here. The "good news"
676 # is that we're already converting everything to long vectors,
677 # and that strategy works for tuples as well.
679 # We pass check=False because the matrix basis is "guaranteed"
680 # to be linearly independent... right? Ha ha.
682 V
= VectorSpace(self
.base_ring(), len(elt
))
683 W
= V
.span_of_basis( (V(_all2list(s
)) for s
in self
.matrix_basis()),
687 coords
= W
.coordinate_vector(V(elt
))
688 except ArithmeticError: # vector is not in free module
689 raise ValueError(msg
)
691 return self
.from_vector(coords
)
695 Return a string representation of ``self``.
699 sage: from mjo.eja.eja_algebra import JordanSpinEJA
703 Ensure that it says what we think it says::
705 sage: JordanSpinEJA(2, field=AA)
706 Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
707 sage: JordanSpinEJA(3, field=RDF)
708 Euclidean Jordan algebra of dimension 3 over Real Double Field
711 fmt
= "Euclidean Jordan algebra of dimension {} over {}"
712 return fmt
.format(self
.dimension(), self
.base_ring())
716 def characteristic_polynomial_of(self
):
718 Return the algebra's "characteristic polynomial of" function,
719 which is itself a multivariate polynomial that, when evaluated
720 at the coordinates of some algebra element, returns that
721 element's characteristic polynomial.
723 The resulting polynomial has `n+1` variables, where `n` is the
724 dimension of this algebra. The first `n` variables correspond to
725 the coordinates of an algebra element: when evaluated at the
726 coordinates of an algebra element with respect to a certain
727 basis, the result is a univariate polynomial (in the one
728 remaining variable ``t``), namely the characteristic polynomial
733 sage: from mjo.eja.eja_algebra import JordanSpinEJA, TrivialEJA
737 The characteristic polynomial in the spin algebra is given in
738 Alizadeh, Example 11.11::
740 sage: J = JordanSpinEJA(3)
741 sage: p = J.characteristic_polynomial_of(); p
742 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
743 sage: xvec = J.one().to_vector()
747 By definition, the characteristic polynomial is a monic
748 degree-zero polynomial in a rank-zero algebra. Note that
749 Cayley-Hamilton is indeed satisfied since the polynomial
750 ``1`` evaluates to the identity element of the algebra on
753 sage: J = TrivialEJA()
754 sage: J.characteristic_polynomial_of()
761 # The list of coefficient polynomials a_0, a_1, a_2, ..., a_(r-1).
762 a
= self
._charpoly
_coefficients
()
764 # We go to a bit of trouble here to reorder the
765 # indeterminates, so that it's easier to evaluate the
766 # characteristic polynomial at x's coordinates and get back
767 # something in terms of t, which is what we want.
768 S
= PolynomialRing(self
.base_ring(),'t')
772 S
= PolynomialRing(S
, R
.variable_names())
775 return (t
**r
+ sum( a
[k
]*(t
**k
) for k
in range(r
) ))
777 def coordinate_polynomial_ring(self
):
779 The multivariate polynomial ring in which this algebra's
780 :meth:`characteristic_polynomial_of` lives.
784 sage: from mjo.eja.eja_algebra import (HadamardEJA,
785 ....: RealSymmetricEJA)
789 sage: J = HadamardEJA(2)
790 sage: J.coordinate_polynomial_ring()
791 Multivariate Polynomial Ring in X1, X2...
792 sage: J = RealSymmetricEJA(3,field=QQ,orthonormalize=False)
793 sage: J.coordinate_polynomial_ring()
794 Multivariate Polynomial Ring in X1, X2, X3, X4, X5, X6...
797 var_names
= tuple( "X%d" % z
for z
in range(1, self
.dimension()+1) )
798 return PolynomialRing(self
.base_ring(), var_names
)
800 def inner_product(self
, x
, y
):
802 The inner product associated with this Euclidean Jordan algebra.
804 Defaults to the trace inner product, but can be overridden by
805 subclasses if they are sure that the necessary properties are
810 sage: from mjo.eja.eja_algebra import (random_eja,
812 ....: BilinearFormEJA)
816 Our inner product is "associative," which means the following for
817 a symmetric bilinear form::
819 sage: set_random_seed()
820 sage: J = random_eja()
821 sage: x,y,z = J.random_elements(3)
822 sage: (x*y).inner_product(z) == y.inner_product(x*z)
827 Ensure that this is the usual inner product for the algebras
830 sage: set_random_seed()
831 sage: J = HadamardEJA.random_instance()
832 sage: x,y = J.random_elements(2)
833 sage: actual = x.inner_product(y)
834 sage: expected = x.to_vector().inner_product(y.to_vector())
835 sage: actual == expected
838 Ensure that this is one-half of the trace inner-product in a
839 BilinearFormEJA that isn't just the reals (when ``n`` isn't
840 one). This is in Faraut and Koranyi, and also my "On the
843 sage: set_random_seed()
844 sage: J = BilinearFormEJA.random_instance()
845 sage: n = J.dimension()
846 sage: x = J.random_element()
847 sage: y = J.random_element()
848 sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
851 B
= self
._inner
_product
_matrix
852 return (B
*x
.to_vector()).inner_product(y
.to_vector())
855 def is_trivial(self
):
857 Return whether or not this algebra is trivial.
859 A trivial algebra contains only the zero element.
863 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
868 sage: J = ComplexHermitianEJA(3)
874 sage: J = TrivialEJA()
879 return self
.dimension() == 0
882 def multiplication_table(self
):
884 Return a visual representation of this algebra's multiplication
885 table (on basis elements).
889 sage: from mjo.eja.eja_algebra import JordanSpinEJA
893 sage: J = JordanSpinEJA(4)
894 sage: J.multiplication_table()
895 +----++----+----+----+----+
896 | * || e0 | e1 | e2 | e3 |
897 +====++====+====+====+====+
898 | e0 || e0 | e1 | e2 | e3 |
899 +----++----+----+----+----+
900 | e1 || e1 | e0 | 0 | 0 |
901 +----++----+----+----+----+
902 | e2 || e2 | 0 | e0 | 0 |
903 +----++----+----+----+----+
904 | e3 || e3 | 0 | 0 | e0 |
905 +----++----+----+----+----+
909 # Prepend the header row.
910 M
= [["*"] + list(self
.gens())]
912 # And to each subsequent row, prepend an entry that belongs to
913 # the left-side "header column."
914 M
+= [ [self
.monomial(i
)] + [ self
.monomial(i
)*self
.monomial(j
)
918 return table(M
, header_row
=True, header_column
=True, frame
=True)
921 def matrix_basis(self
):
923 Return an (often more natural) representation of this algebras
924 basis as an ordered tuple of matrices.
926 Every finite-dimensional Euclidean Jordan Algebra is a, up to
927 Jordan isomorphism, a direct sum of five simple
928 algebras---four of which comprise Hermitian matrices. And the
929 last type of algebra can of course be thought of as `n`-by-`1`
930 column matrices (ambiguusly called column vectors) to avoid
931 special cases. As a result, matrices (and column vectors) are
932 a natural representation format for Euclidean Jordan algebra
935 But, when we construct an algebra from a basis of matrices,
936 those matrix representations are lost in favor of coordinate
937 vectors *with respect to* that basis. We could eventually
938 convert back if we tried hard enough, but having the original
939 representations handy is valuable enough that we simply store
940 them and return them from this method.
942 Why implement this for non-matrix algebras? Avoiding special
943 cases for the :class:`BilinearFormEJA` pays with simplicity in
944 its own right. But mainly, we would like to be able to assume
945 that elements of a :class:`CartesianProductEJA` can be displayed
946 nicely, without having to have special classes for direct sums
947 one of whose components was a matrix algebra.
951 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
952 ....: RealSymmetricEJA)
956 sage: J = RealSymmetricEJA(2)
958 Finite family {0: e0, 1: e1, 2: e2}
959 sage: J.matrix_basis()
961 [1 0] [ 0 0.7071067811865475?] [0 0]
962 [0 0], [0.7071067811865475? 0], [0 1]
967 sage: J = JordanSpinEJA(2)
969 Finite family {0: e0, 1: e1}
970 sage: J.matrix_basis()
976 return self
._matrix
_basis
979 def matrix_space(self
):
981 Return the matrix space in which this algebra's elements live, if
982 we think of them as matrices (including column vectors of the
985 "By default" this will be an `n`-by-`1` column-matrix space,
986 except when the algebra is trivial. There it's `n`-by-`n`
987 (where `n` is zero), to ensure that two elements of the matrix
988 space (empty matrices) can be multiplied. For algebras of
989 matrices, this returns the space in which their
990 real embeddings live.
994 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
996 ....: QuaternionHermitianEJA,
1001 By default, the matrix representation is just a column-matrix
1002 equivalent to the vector representation::
1004 sage: J = JordanSpinEJA(3)
1005 sage: J.matrix_space()
1006 Full MatrixSpace of 3 by 1 dense matrices over Algebraic
1009 The matrix representation in the trivial algebra is
1010 zero-by-zero instead of the usual `n`-by-one::
1012 sage: J = TrivialEJA()
1013 sage: J.matrix_space()
1014 Full MatrixSpace of 0 by 0 dense matrices over Algebraic
1017 The matrix space for complex/quaternion Hermitian matrix EJA
1018 is the space in which their real-embeddings live, not the
1019 original complex/quaternion matrix space::
1021 sage: J = ComplexHermitianEJA(2,field=QQ,orthonormalize=False)
1022 sage: J.matrix_space()
1023 Full MatrixSpace of 4 by 4 dense matrices over Rational Field
1024 sage: J = QuaternionHermitianEJA(1,field=QQ,orthonormalize=False)
1025 sage: J.matrix_space()
1026 Full MatrixSpace of 4 by 4 dense matrices over Rational Field
1029 if self
.is_trivial():
1030 return MatrixSpace(self
.base_ring(), 0)
1032 return self
.matrix_basis()[0].parent()
1038 Return the unit element of this algebra.
1042 sage: from mjo.eja.eja_algebra import (HadamardEJA,
1047 We can compute unit element in the Hadamard EJA::
1049 sage: J = HadamardEJA(5)
1051 e0 + e1 + e2 + e3 + e4
1053 The unit element in the Hadamard EJA is inherited in the
1054 subalgebras generated by its elements::
1056 sage: J = HadamardEJA(5)
1058 e0 + e1 + e2 + e3 + e4
1059 sage: x = sum(J.gens())
1060 sage: A = x.subalgebra_generated_by(orthonormalize=False)
1063 sage: A.one().superalgebra_element()
1064 e0 + e1 + e2 + e3 + e4
1068 The identity element acts like the identity, regardless of
1069 whether or not we orthonormalize::
1071 sage: set_random_seed()
1072 sage: J = random_eja()
1073 sage: x = J.random_element()
1074 sage: J.one()*x == x and x*J.one() == x
1076 sage: A = x.subalgebra_generated_by()
1077 sage: y = A.random_element()
1078 sage: A.one()*y == y and y*A.one() == y
1083 sage: set_random_seed()
1084 sage: J = random_eja(field=QQ, orthonormalize=False)
1085 sage: x = J.random_element()
1086 sage: J.one()*x == x and x*J.one() == x
1088 sage: A = x.subalgebra_generated_by(orthonormalize=False)
1089 sage: y = A.random_element()
1090 sage: A.one()*y == y and y*A.one() == y
1093 The matrix of the unit element's operator is the identity,
1094 regardless of the base field and whether or not we
1097 sage: set_random_seed()
1098 sage: J = random_eja()
1099 sage: actual = J.one().operator().matrix()
1100 sage: expected = matrix.identity(J.base_ring(), J.dimension())
1101 sage: actual == expected
1103 sage: x = J.random_element()
1104 sage: A = x.subalgebra_generated_by()
1105 sage: actual = A.one().operator().matrix()
1106 sage: expected = matrix.identity(A.base_ring(), A.dimension())
1107 sage: actual == expected
1112 sage: set_random_seed()
1113 sage: J = random_eja(field=QQ, orthonormalize=False)
1114 sage: actual = J.one().operator().matrix()
1115 sage: expected = matrix.identity(J.base_ring(), J.dimension())
1116 sage: actual == expected
1118 sage: x = J.random_element()
1119 sage: A = x.subalgebra_generated_by(orthonormalize=False)
1120 sage: actual = A.one().operator().matrix()
1121 sage: expected = matrix.identity(A.base_ring(), A.dimension())
1122 sage: actual == expected
1125 Ensure that the cached unit element (often precomputed by
1126 hand) agrees with the computed one::
1128 sage: set_random_seed()
1129 sage: J = random_eja()
1130 sage: cached = J.one()
1131 sage: J.one.clear_cache()
1132 sage: J.one() == cached
1137 sage: set_random_seed()
1138 sage: J = random_eja(field=QQ, orthonormalize=False)
1139 sage: cached = J.one()
1140 sage: J.one.clear_cache()
1141 sage: J.one() == cached
1145 # We can brute-force compute the matrices of the operators
1146 # that correspond to the basis elements of this algebra.
1147 # If some linear combination of those basis elements is the
1148 # algebra identity, then the same linear combination of
1149 # their matrices has to be the identity matrix.
1151 # Of course, matrices aren't vectors in sage, so we have to
1152 # appeal to the "long vectors" isometry.
1153 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
1155 # Now we use basic linear algebra to find the coefficients,
1156 # of the matrices-as-vectors-linear-combination, which should
1157 # work for the original algebra basis too.
1158 A
= matrix(self
.base_ring(), oper_vecs
)
1160 # We used the isometry on the left-hand side already, but we
1161 # still need to do it for the right-hand side. Recall that we
1162 # wanted something that summed to the identity matrix.
1163 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
1165 # Now if there's an identity element in the algebra, this
1166 # should work. We solve on the left to avoid having to
1167 # transpose the matrix "A".
1168 return self
.from_vector(A
.solve_left(b
))
1171 def peirce_decomposition(self
, c
):
1173 The Peirce decomposition of this algebra relative to the
1176 In the future, this can be extended to a complete system of
1177 orthogonal idempotents.
1181 - ``c`` -- an idempotent of this algebra.
1185 A triple (J0, J5, J1) containing two subalgebras and one subspace
1188 - ``J0`` -- the algebra on the eigenspace of ``c.operator()``
1189 corresponding to the eigenvalue zero.
1191 - ``J5`` -- the eigenspace (NOT a subalgebra) of ``c.operator()``
1192 corresponding to the eigenvalue one-half.
1194 - ``J1`` -- the algebra on the eigenspace of ``c.operator()``
1195 corresponding to the eigenvalue one.
1197 These are the only possible eigenspaces for that operator, and this
1198 algebra is a direct sum of them. The spaces ``J0`` and ``J1`` are
1199 orthogonal, and are subalgebras of this algebra with the appropriate
1204 sage: from mjo.eja.eja_algebra import random_eja, RealSymmetricEJA
1208 The canonical example comes from the symmetric matrices, which
1209 decompose into diagonal and off-diagonal parts::
1211 sage: J = RealSymmetricEJA(3)
1212 sage: C = matrix(QQ, [ [1,0,0],
1216 sage: J0,J5,J1 = J.peirce_decomposition(c)
1218 Euclidean Jordan algebra of dimension 1...
1220 Vector space of degree 6 and dimension 2...
1222 Euclidean Jordan algebra of dimension 3...
1223 sage: J0.one().to_matrix()
1227 sage: orig_df = AA.options.display_format
1228 sage: AA.options.display_format = 'radical'
1229 sage: J.from_vector(J5.basis()[0]).to_matrix()
1233 sage: J.from_vector(J5.basis()[1]).to_matrix()
1237 sage: AA.options.display_format = orig_df
1238 sage: J1.one().to_matrix()
1245 Every algebra decomposes trivially with respect to its identity
1248 sage: set_random_seed()
1249 sage: J = random_eja()
1250 sage: J0,J5,J1 = J.peirce_decomposition(J.one())
1251 sage: J0.dimension() == 0 and J5.dimension() == 0
1253 sage: J1.superalgebra() == J and J1.dimension() == J.dimension()
1256 The decomposition is into eigenspaces, and its components are
1257 therefore necessarily orthogonal. Moreover, the identity
1258 elements in the two subalgebras are the projections onto their
1259 respective subspaces of the superalgebra's identity element::
1261 sage: set_random_seed()
1262 sage: J = random_eja()
1263 sage: x = J.random_element()
1264 sage: if not J.is_trivial():
1265 ....: while x.is_nilpotent():
1266 ....: x = J.random_element()
1267 sage: c = x.subalgebra_idempotent()
1268 sage: J0,J5,J1 = J.peirce_decomposition(c)
1270 sage: for (w,y,z) in zip(J0.basis(), J5.basis(), J1.basis()):
1271 ....: w = w.superalgebra_element()
1272 ....: y = J.from_vector(y)
1273 ....: z = z.superalgebra_element()
1274 ....: ipsum += w.inner_product(y).abs()
1275 ....: ipsum += w.inner_product(z).abs()
1276 ....: ipsum += y.inner_product(z).abs()
1279 sage: J1(c) == J1.one()
1281 sage: J0(J.one() - c) == J0.one()
1285 if not c
.is_idempotent():
1286 raise ValueError("element is not idempotent: %s" % c
)
1288 # Default these to what they should be if they turn out to be
1289 # trivial, because eigenspaces_left() won't return eigenvalues
1290 # corresponding to trivial spaces (e.g. it returns only the
1291 # eigenspace corresponding to lambda=1 if you take the
1292 # decomposition relative to the identity element).
1293 trivial
= self
.subalgebra(())
1294 J0
= trivial
# eigenvalue zero
1295 J5
= VectorSpace(self
.base_ring(), 0) # eigenvalue one-half
1296 J1
= trivial
# eigenvalue one
1298 for (eigval
, eigspace
) in c
.operator().matrix().right_eigenspaces():
1299 if eigval
== ~
(self
.base_ring()(2)):
1302 gens
= tuple( self
.from_vector(b
) for b
in eigspace
.basis() )
1303 subalg
= self
.subalgebra(gens
, check_axioms
=False)
1309 raise ValueError("unexpected eigenvalue: %s" % eigval
)
1314 def random_element(self
, thorough
=False):
1316 Return a random element of this algebra.
1318 Our algebra superclass method only returns a linear
1319 combination of at most two basis elements. We instead
1320 want the vector space "random element" method that
1321 returns a more diverse selection.
1325 - ``thorough`` -- (boolean; default False) whether or not we
1326 should generate irrational coefficients for the random
1327 element when our base ring is irrational; this slows the
1328 algebra operations to a crawl, but any truly random method
1332 # For a general base ring... maybe we can trust this to do the
1333 # right thing? Unlikely, but.
1334 V
= self
.vector_space()
1335 v
= V
.random_element()
1337 if self
.base_ring() is AA
:
1338 # The "random element" method of the algebraic reals is
1339 # stupid at the moment, and only returns integers between
1340 # -2 and 2, inclusive:
1342 # https://trac.sagemath.org/ticket/30875
1344 # Instead, we implement our own "random vector" method,
1345 # and then coerce that into the algebra. We use the vector
1346 # space degree here instead of the dimension because a
1347 # subalgebra could (for example) be spanned by only two
1348 # vectors, each with five coordinates. We need to
1349 # generate all five coordinates.
1351 v
*= QQbar
.random_element().real()
1353 v
*= QQ
.random_element()
1355 return self
.from_vector(V
.coordinate_vector(v
))
1357 def random_elements(self
, count
, thorough
=False):
1359 Return ``count`` random elements as a tuple.
1363 - ``thorough`` -- (boolean; default False) whether or not we
1364 should generate irrational coefficients for the random
1365 elements when our base ring is irrational; this slows the
1366 algebra operations to a crawl, but any truly random method
1371 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1375 sage: J = JordanSpinEJA(3)
1376 sage: x,y,z = J.random_elements(3)
1377 sage: all( [ x in J, y in J, z in J ])
1379 sage: len( J.random_elements(10) ) == 10
1383 return tuple( self
.random_element(thorough
)
1384 for idx
in range(count
) )
1388 def _charpoly_coefficients(self
):
1390 The `r` polynomial coefficients of the "characteristic polynomial
1395 sage: from mjo.eja.eja_algebra import random_eja
1399 The theory shows that these are all homogeneous polynomials of
1402 sage: set_random_seed()
1403 sage: J = random_eja()
1404 sage: all(p.is_homogeneous() for p in J._charpoly_coefficients())
1408 n
= self
.dimension()
1409 R
= self
.coordinate_polynomial_ring()
1411 F
= R
.fraction_field()
1414 # From a result in my book, these are the entries of the
1415 # basis representation of L_x.
1416 return sum( vars[k
]*self
.monomial(k
).operator().matrix()[i
,j
]
1419 L_x
= matrix(F
, n
, n
, L_x_i_j
)
1422 if self
.rank
.is_in_cache():
1424 # There's no need to pad the system with redundant
1425 # columns if we *know* they'll be redundant.
1428 # Compute an extra power in case the rank is equal to
1429 # the dimension (otherwise, we would stop at x^(r-1)).
1430 x_powers
= [ (L_x
**k
)*self
.one().to_vector()
1431 for k
in range(n
+1) ]
1432 A
= matrix
.column(F
, x_powers
[:n
])
1433 AE
= A
.extended_echelon_form()
1440 # The theory says that only the first "r" coefficients are
1441 # nonzero, and they actually live in the original polynomial
1442 # ring and not the fraction field. We negate them because in
1443 # the actual characteristic polynomial, they get moved to the
1444 # other side where x^r lives. We don't bother to trim A_rref
1445 # down to a square matrix and solve the resulting system,
1446 # because the upper-left r-by-r portion of A_rref is
1447 # guaranteed to be the identity matrix, so e.g.
1449 # A_rref.solve_right(Y)
1451 # would just be returning Y.
1452 return (-E
*b
)[:r
].change_ring(R
)
1457 Return the rank of this EJA.
1459 This is a cached method because we know the rank a priori for
1460 all of the algebras we can construct. Thus we can avoid the
1461 expensive ``_charpoly_coefficients()`` call unless we truly
1462 need to compute the whole characteristic polynomial.
1466 sage: from mjo.eja.eja_algebra import (HadamardEJA,
1467 ....: JordanSpinEJA,
1468 ....: RealSymmetricEJA,
1469 ....: ComplexHermitianEJA,
1470 ....: QuaternionHermitianEJA,
1475 The rank of the Jordan spin algebra is always two::
1477 sage: JordanSpinEJA(2).rank()
1479 sage: JordanSpinEJA(3).rank()
1481 sage: JordanSpinEJA(4).rank()
1484 The rank of the `n`-by-`n` Hermitian real, complex, or
1485 quaternion matrices is `n`::
1487 sage: RealSymmetricEJA(4).rank()
1489 sage: ComplexHermitianEJA(3).rank()
1491 sage: QuaternionHermitianEJA(2).rank()
1496 Ensure that every EJA that we know how to construct has a
1497 positive integer rank, unless the algebra is trivial in
1498 which case its rank will be zero::
1500 sage: set_random_seed()
1501 sage: J = random_eja()
1505 sage: r > 0 or (r == 0 and J.is_trivial())
1508 Ensure that computing the rank actually works, since the ranks
1509 of all simple algebras are known and will be cached by default::
1511 sage: set_random_seed() # long time
1512 sage: J = random_eja() # long time
1513 sage: cached = J.rank() # long time
1514 sage: J.rank.clear_cache() # long time
1515 sage: J.rank() == cached # long time
1519 return len(self
._charpoly
_coefficients
())
1522 def subalgebra(self
, basis
, **kwargs
):
1524 Create a subalgebra of this algebra from the given basis.
1526 from mjo
.eja
.eja_subalgebra
import FiniteDimensionalEJASubalgebra
1527 return FiniteDimensionalEJASubalgebra(self
, basis
, **kwargs
)
1530 def vector_space(self
):
1532 Return the vector space that underlies this algebra.
1536 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1540 sage: J = RealSymmetricEJA(2)
1541 sage: J.vector_space()
1542 Vector space of dimension 3 over...
1545 return self
.zero().to_vector().parent().ambient_vector_space()
1549 class RationalBasisEJA(FiniteDimensionalEJA
):
1551 New class for algebras whose supplied basis elements have all rational entries.
1555 sage: from mjo.eja.eja_algebra import BilinearFormEJA
1559 The supplied basis is orthonormalized by default::
1561 sage: B = matrix(QQ, [[1, 0, 0], [0, 25, -32], [0, -32, 41]])
1562 sage: J = BilinearFormEJA(B)
1563 sage: J.matrix_basis()
1580 # Abuse the check_field parameter to check that the entries of
1581 # out basis (in ambient coordinates) are in the field QQ.
1582 if not all( all(b_i
in QQ
for b_i
in b
.list()) for b
in basis
):
1583 raise TypeError("basis not rational")
1585 super().__init
__(basis
,
1589 check_field
=check_field
,
1592 self
._rational
_algebra
= None
1594 # There's no point in constructing the extra algebra if this
1595 # one is already rational.
1597 # Note: the same Jordan and inner-products work here,
1598 # because they are necessarily defined with respect to
1599 # ambient coordinates and not any particular basis.
1600 self
._rational
_algebra
= FiniteDimensionalEJA(
1605 associative
=self
.is_associative(),
1606 orthonormalize
=False,
1611 def _charpoly_coefficients(self
):
1615 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
1616 ....: JordanSpinEJA)
1620 The base ring of the resulting polynomial coefficients is what
1621 it should be, and not the rationals (unless the algebra was
1622 already over the rationals)::
1624 sage: J = JordanSpinEJA(3)
1625 sage: J._charpoly_coefficients()
1626 (X1^2 - X2^2 - X3^2, -2*X1)
1627 sage: a0 = J._charpoly_coefficients()[0]
1629 Algebraic Real Field
1630 sage: a0.base_ring()
1631 Algebraic Real Field
1634 if self
._rational
_algebra
is None:
1635 # There's no need to construct *another* algebra over the
1636 # rationals if this one is already over the
1637 # rationals. Likewise, if we never orthonormalized our
1638 # basis, we might as well just use the given one.
1639 return super()._charpoly
_coefficients
()
1641 # Do the computation over the rationals. The answer will be
1642 # the same, because all we've done is a change of basis.
1643 # Then, change back from QQ to our real base ring
1644 a
= ( a_i
.change_ring(self
.base_ring())
1645 for a_i
in self
._rational
_algebra
._charpoly
_coefficients
() )
1647 if self
._deortho
_matrix
is None:
1648 # This can happen if our base ring was, say, AA and we
1649 # chose not to (or didn't need to) orthonormalize. It's
1650 # still faster to do the computations over QQ even if
1651 # the numbers in the boxes stay the same.
1654 # Otherwise, convert the coordinate variables back to the
1655 # deorthonormalized ones.
1656 R
= self
.coordinate_polynomial_ring()
1657 from sage
.modules
.free_module_element
import vector
1658 X
= vector(R
, R
.gens())
1659 BX
= self
._deortho
_matrix
*X
1661 subs_dict
= { X[i]: BX[i] for i in range(len(X)) }
1662 return tuple( a_i
.subs(subs_dict
) for a_i
in a
)
1664 class ConcreteEJA(RationalBasisEJA
):
1666 A class for the Euclidean Jordan algebras that we know by name.
1668 These are the Jordan algebras whose basis, multiplication table,
1669 rank, and so on are known a priori. More to the point, they are
1670 the Euclidean Jordan algebras for which we are able to conjure up
1671 a "random instance."
1675 sage: from mjo.eja.eja_algebra import ConcreteEJA
1679 Our basis is normalized with respect to the algebra's inner
1680 product, unless we specify otherwise::
1682 sage: set_random_seed()
1683 sage: J = ConcreteEJA.random_instance()
1684 sage: all( b.norm() == 1 for b in J.gens() )
1687 Since our basis is orthonormal with respect to the algebra's inner
1688 product, and since we know that this algebra is an EJA, any
1689 left-multiplication operator's matrix will be symmetric because
1690 natural->EJA basis representation is an isometry and within the
1691 EJA the operator is self-adjoint by the Jordan axiom::
1693 sage: set_random_seed()
1694 sage: J = ConcreteEJA.random_instance()
1695 sage: x = J.random_element()
1696 sage: x.operator().is_self_adjoint()
1701 def _max_random_instance_size():
1703 Return an integer "size" that is an upper bound on the size of
1704 this algebra when it is used in a random test
1705 case. Unfortunately, the term "size" is ambiguous -- when
1706 dealing with `R^n` under either the Hadamard or Jordan spin
1707 product, the "size" refers to the dimension `n`. When dealing
1708 with a matrix algebra (real symmetric or complex/quaternion
1709 Hermitian), it refers to the size of the matrix, which is far
1710 less than the dimension of the underlying vector space.
1712 This method must be implemented in each subclass.
1714 raise NotImplementedError
1717 def random_instance(cls
, *args
, **kwargs
):
1719 Return a random instance of this type of algebra.
1721 This method should be implemented in each subclass.
1723 from sage
.misc
.prandom
import choice
1724 eja_class
= choice(cls
.__subclasses
__())
1726 # These all bubble up to the RationalBasisEJA superclass
1727 # constructor, so any (kw)args valid there are also valid
1729 return eja_class
.random_instance(*args
, **kwargs
)
1734 def dimension_over_reals():
1736 The dimension of this matrix's base ring over the reals.
1738 The reals are dimension one over themselves, obviously; that's
1739 just `\mathbb{R}^{1}`. Likewise, the complex numbers `a + bi`
1740 have dimension two. Finally, the quaternions have dimension
1741 four over the reals.
1743 This is used to determine the size of the matrix returned from
1744 :meth:`real_embed`, among other things.
1746 raise NotImplementedError
1749 def real_embed(cls
,M
):
1751 Embed the matrix ``M`` into a space of real matrices.
1753 The matrix ``M`` can have entries in any field at the moment:
1754 the real numbers, complex numbers, or quaternions. And although
1755 they are not a field, we can probably support octonions at some
1756 point, too. This function returns a real matrix that "acts like"
1757 the original with respect to matrix multiplication; i.e.
1759 real_embed(M*N) = real_embed(M)*real_embed(N)
1762 if M
.ncols() != M
.nrows():
1763 raise ValueError("the matrix 'M' must be square")
1768 def real_unembed(cls
,M
):
1770 The inverse of :meth:`real_embed`.
1772 if M
.ncols() != M
.nrows():
1773 raise ValueError("the matrix 'M' must be square")
1774 if not ZZ(M
.nrows()).mod(cls
.dimension_over_reals()).is_zero():
1775 raise ValueError("the matrix 'M' must be a real embedding")
1779 def jordan_product(X
,Y
):
1780 return (X
*Y
+ Y
*X
)/2
1783 def trace_inner_product(cls
,X
,Y
):
1785 Compute the trace inner-product of two real-embeddings.
1789 sage: from mjo.eja.eja_algebra import (RealSymmetricEJA,
1790 ....: ComplexHermitianEJA,
1791 ....: QuaternionHermitianEJA)
1795 This gives the same answer as it would if we computed the trace
1796 from the unembedded (original) matrices::
1798 sage: set_random_seed()
1799 sage: J = RealSymmetricEJA.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()
1806 sage: actual = J.trace_inner_product(Xe,Ye)
1807 sage: actual == expected
1812 sage: set_random_seed()
1813 sage: J = ComplexHermitianEJA.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().real()
1820 sage: actual = J.trace_inner_product(Xe,Ye)
1821 sage: actual == expected
1826 sage: set_random_seed()
1827 sage: J = QuaternionHermitianEJA.random_instance()
1828 sage: x,y = J.random_elements(2)
1829 sage: Xe = x.to_matrix()
1830 sage: Ye = y.to_matrix()
1831 sage: X = J.real_unembed(Xe)
1832 sage: Y = J.real_unembed(Ye)
1833 sage: expected = (X*Y).trace().coefficient_tuple()[0]
1834 sage: actual = J.trace_inner_product(Xe,Ye)
1835 sage: actual == expected
1839 # This does in fact compute the real part of the trace.
1840 # If we compute the trace of e.g. a complex matrix M,
1841 # then we do so by adding up its diagonal entries --
1842 # call them z_1 through z_n. The real embedding of z_1
1843 # will be a 2-by-2 REAL matrix [a, b; -b, a] whose trace
1844 # as a REAL matrix will be 2*a = 2*Re(z_1). And so forth.
1845 return (X
*Y
).trace()/cls
.dimension_over_reals()
1848 class RealMatrixEJA(MatrixEJA
):
1850 def dimension_over_reals():
1854 class RealSymmetricEJA(ConcreteEJA
, RealMatrixEJA
):
1856 The rank-n simple EJA consisting of real symmetric n-by-n
1857 matrices, the usual symmetric Jordan product, and the trace inner
1858 product. It has dimension `(n^2 + n)/2` over the reals.
1862 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1866 sage: J = RealSymmetricEJA(2)
1867 sage: e0, e1, e2 = J.gens()
1875 In theory, our "field" can be any subfield of the reals::
1877 sage: RealSymmetricEJA(2, field=RDF, check_axioms=True)
1878 Euclidean Jordan algebra of dimension 3 over Real Double Field
1879 sage: RealSymmetricEJA(2, field=RR, check_axioms=True)
1880 Euclidean Jordan algebra of dimension 3 over Real Field with
1881 53 bits of precision
1885 The dimension of this algebra is `(n^2 + n) / 2`::
1887 sage: set_random_seed()
1888 sage: n_max = RealSymmetricEJA._max_random_instance_size()
1889 sage: n = ZZ.random_element(1, n_max)
1890 sage: J = RealSymmetricEJA(n)
1891 sage: J.dimension() == (n^2 + n)/2
1894 The Jordan multiplication is what we think it is::
1896 sage: set_random_seed()
1897 sage: J = RealSymmetricEJA.random_instance()
1898 sage: x,y = J.random_elements(2)
1899 sage: actual = (x*y).to_matrix()
1900 sage: X = x.to_matrix()
1901 sage: Y = y.to_matrix()
1902 sage: expected = (X*Y + Y*X)/2
1903 sage: actual == expected
1905 sage: J(expected) == x*y
1908 We can change the generator prefix::
1910 sage: RealSymmetricEJA(3, prefix='q').gens()
1911 (q0, q1, q2, q3, q4, q5)
1913 We can construct the (trivial) algebra of rank zero::
1915 sage: RealSymmetricEJA(0)
1916 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1920 def _denormalized_basis(cls
, n
):
1922 Return a basis for the space of real symmetric n-by-n matrices.
1926 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1930 sage: set_random_seed()
1931 sage: n = ZZ.random_element(1,5)
1932 sage: B = RealSymmetricEJA._denormalized_basis(n)
1933 sage: all( M.is_symmetric() for M in B)
1937 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1941 for j
in range(i
+1):
1942 Eij
= matrix(ZZ
, n
, lambda k
,l
: k
==i
and l
==j
)
1946 Sij
= Eij
+ Eij
.transpose()
1952 def _max_random_instance_size():
1953 return 4 # Dimension 10
1956 def random_instance(cls
, **kwargs
):
1958 Return a random instance of this type of algebra.
1960 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1961 return cls(n
, **kwargs
)
1963 def __init__(self
, n
, **kwargs
):
1964 # We know this is a valid EJA, but will double-check
1965 # if the user passes check_axioms=True.
1966 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
1972 super().__init
__(self
._denormalized
_basis
(n
),
1973 self
.jordan_product
,
1974 self
.trace_inner_product
,
1975 associative
=associative
,
1978 # TODO: this could be factored out somehow, but is left here
1979 # because the MatrixEJA is not presently a subclass of the
1980 # FDEJA class that defines rank() and one().
1981 self
.rank
.set_cache(n
)
1982 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
1983 self
.one
.set_cache(self(idV
))
1987 class ComplexMatrixEJA(MatrixEJA
):
1988 # A manual dictionary-cache for the complex_extension() method,
1989 # since apparently @classmethods can't also be @cached_methods.
1990 _complex_extension
= {}
1993 def complex_extension(cls
,field
):
1995 The complex field that we embed/unembed, as an extension
1996 of the given ``field``.
1998 if field
in cls
._complex
_extension
:
1999 return cls
._complex
_extension
[field
]
2001 # Sage doesn't know how to adjoin the complex "i" (the root of
2002 # x^2 + 1) to a field in a general way. Here, we just enumerate
2003 # all of the cases that I have cared to support so far.
2005 # Sage doesn't know how to embed AA into QQbar, i.e. how
2006 # to adjoin sqrt(-1) to AA.
2008 elif not field
.is_exact():
2010 F
= field
.complex_field()
2012 # Works for QQ and... maybe some other fields.
2013 R
= PolynomialRing(field
, 'z')
2015 F
= field
.extension(z
**2 + 1, 'I', embedding
=CLF(-1).sqrt())
2017 cls
._complex
_extension
[field
] = F
2021 def dimension_over_reals():
2025 def real_embed(cls
,M
):
2027 Embed the n-by-n complex matrix ``M`` into the space of real
2028 matrices of size 2n-by-2n via the map the sends each entry `z = a +
2029 bi` to the block matrix ``[[a,b],[-b,a]]``.
2033 sage: from mjo.eja.eja_algebra import ComplexMatrixEJA
2037 sage: F = QuadraticField(-1, 'I')
2038 sage: x1 = F(4 - 2*i)
2039 sage: x2 = F(1 + 2*i)
2042 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
2043 sage: ComplexMatrixEJA.real_embed(M)
2052 Embedding is a homomorphism (isomorphism, in fact)::
2054 sage: set_random_seed()
2055 sage: n = ZZ.random_element(3)
2056 sage: F = QuadraticField(-1, 'I')
2057 sage: X = random_matrix(F, n)
2058 sage: Y = random_matrix(F, n)
2059 sage: Xe = ComplexMatrixEJA.real_embed(X)
2060 sage: Ye = ComplexMatrixEJA.real_embed(Y)
2061 sage: XYe = ComplexMatrixEJA.real_embed(X*Y)
2066 super().real_embed(M
)
2069 # We don't need any adjoined elements...
2070 field
= M
.base_ring().base_ring()
2076 blocks
.append(matrix(field
, 2, [ [ a
, b
],
2079 return matrix
.block(field
, n
, blocks
)
2083 def real_unembed(cls
,M
):
2085 The inverse of _embed_complex_matrix().
2089 sage: from mjo.eja.eja_algebra import ComplexMatrixEJA
2093 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
2094 ....: [-2, 1, -4, 3],
2095 ....: [ 9, 10, 11, 12],
2096 ....: [-10, 9, -12, 11] ])
2097 sage: ComplexMatrixEJA.real_unembed(A)
2099 [ 10*I + 9 12*I + 11]
2103 Unembedding is the inverse of embedding::
2105 sage: set_random_seed()
2106 sage: F = QuadraticField(-1, 'I')
2107 sage: M = random_matrix(F, 3)
2108 sage: Me = ComplexMatrixEJA.real_embed(M)
2109 sage: ComplexMatrixEJA.real_unembed(Me) == M
2113 super().real_unembed(M
)
2115 d
= cls
.dimension_over_reals()
2116 F
= cls
.complex_extension(M
.base_ring())
2119 # Go top-left to bottom-right (reading order), converting every
2120 # 2-by-2 block we see to a single complex element.
2122 for k
in range(n
/d
):
2123 for j
in range(n
/d
):
2124 submat
= M
[d
*k
:d
*k
+d
,d
*j
:d
*j
+d
]
2125 if submat
[0,0] != submat
[1,1]:
2126 raise ValueError('bad on-diagonal submatrix')
2127 if submat
[0,1] != -submat
[1,0]:
2128 raise ValueError('bad off-diagonal submatrix')
2129 z
= submat
[0,0] + submat
[0,1]*i
2132 return matrix(F
, n
/d
, elements
)
2135 class ComplexHermitianEJA(ConcreteEJA
, ComplexMatrixEJA
):
2137 The rank-n simple EJA consisting of complex Hermitian n-by-n
2138 matrices over the real numbers, the usual symmetric Jordan product,
2139 and the real-part-of-trace inner product. It has dimension `n^2` over
2144 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
2148 In theory, our "field" can be any subfield of the reals::
2150 sage: ComplexHermitianEJA(2, field=RDF, check_axioms=True)
2151 Euclidean Jordan algebra of dimension 4 over Real Double Field
2152 sage: ComplexHermitianEJA(2, field=RR, check_axioms=True)
2153 Euclidean Jordan algebra of dimension 4 over Real Field with
2154 53 bits of precision
2158 The dimension of this algebra is `n^2`::
2160 sage: set_random_seed()
2161 sage: n_max = ComplexHermitianEJA._max_random_instance_size()
2162 sage: n = ZZ.random_element(1, n_max)
2163 sage: J = ComplexHermitianEJA(n)
2164 sage: J.dimension() == n^2
2167 The Jordan multiplication is what we think it is::
2169 sage: set_random_seed()
2170 sage: J = ComplexHermitianEJA.random_instance()
2171 sage: x,y = J.random_elements(2)
2172 sage: actual = (x*y).to_matrix()
2173 sage: X = x.to_matrix()
2174 sage: Y = y.to_matrix()
2175 sage: expected = (X*Y + Y*X)/2
2176 sage: actual == expected
2178 sage: J(expected) == x*y
2181 We can change the generator prefix::
2183 sage: ComplexHermitianEJA(2, prefix='z').gens()
2186 We can construct the (trivial) algebra of rank zero::
2188 sage: ComplexHermitianEJA(0)
2189 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2194 def _denormalized_basis(cls
, n
):
2196 Returns a basis for the space of complex Hermitian n-by-n matrices.
2198 Why do we embed these? Basically, because all of numerical linear
2199 algebra assumes that you're working with vectors consisting of `n`
2200 entries from a field and scalars from the same field. There's no way
2201 to tell SageMath that (for example) the vectors contain complex
2202 numbers, while the scalar field is real.
2206 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
2210 sage: set_random_seed()
2211 sage: n = ZZ.random_element(1,5)
2212 sage: B = ComplexHermitianEJA._denormalized_basis(n)
2213 sage: all( M.is_symmetric() for M in B)
2218 R
= PolynomialRing(field
, 'z')
2220 F
= field
.extension(z
**2 + 1, 'I')
2223 # This is like the symmetric case, but we need to be careful:
2225 # * We want conjugate-symmetry, not just symmetry.
2226 # * The diagonal will (as a result) be real.
2229 Eij
= matrix
.zero(F
,n
)
2231 for j
in range(i
+1):
2235 Sij
= cls
.real_embed(Eij
)
2238 # The second one has a minus because it's conjugated.
2239 Eij
[j
,i
] = 1 # Eij = Eij + Eij.transpose()
2240 Sij_real
= cls
.real_embed(Eij
)
2242 # Eij = I*Eij - I*Eij.transpose()
2245 Sij_imag
= cls
.real_embed(Eij
)
2251 # Since we embedded these, we can drop back to the "field" that we
2252 # started with instead of the complex extension "F".
2253 return tuple( s
.change_ring(field
) for s
in S
)
2256 def __init__(self
, n
, **kwargs
):
2257 # We know this is a valid EJA, but will double-check
2258 # if the user passes check_axioms=True.
2259 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2265 super().__init
__(self
._denormalized
_basis
(n
),
2266 self
.jordan_product
,
2267 self
.trace_inner_product
,
2268 associative
=associative
,
2270 # TODO: this could be factored out somehow, but is left here
2271 # because the MatrixEJA is not presently a subclass of the
2272 # FDEJA class that defines rank() and one().
2273 self
.rank
.set_cache(n
)
2274 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
2275 self
.one
.set_cache(self(idV
))
2278 def _max_random_instance_size():
2279 return 3 # Dimension 9
2282 def random_instance(cls
, **kwargs
):
2284 Return a random instance of this type of algebra.
2286 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2287 return cls(n
, **kwargs
)
2289 class QuaternionMatrixEJA(MatrixEJA
):
2291 # A manual dictionary-cache for the quaternion_extension() method,
2292 # since apparently @classmethods can't also be @cached_methods.
2293 _quaternion_extension
= {}
2296 def quaternion_extension(cls
,field
):
2298 The quaternion field that we embed/unembed, as an extension
2299 of the given ``field``.
2301 if field
in cls
._quaternion
_extension
:
2302 return cls
._quaternion
_extension
[field
]
2304 Q
= QuaternionAlgebra(field
,-1,-1)
2306 cls
._quaternion
_extension
[field
] = Q
2310 def dimension_over_reals():
2314 def real_embed(cls
,M
):
2316 Embed the n-by-n quaternion matrix ``M`` into the space of real
2317 matrices of size 4n-by-4n by first sending each quaternion entry `z
2318 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
2319 c+di],[-c + di, a-bi]]`, and then embedding those into a real
2324 sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA
2328 sage: Q = QuaternionAlgebra(QQ,-1,-1)
2329 sage: i,j,k = Q.gens()
2330 sage: x = 1 + 2*i + 3*j + 4*k
2331 sage: M = matrix(Q, 1, [[x]])
2332 sage: QuaternionMatrixEJA.real_embed(M)
2338 Embedding is a homomorphism (isomorphism, in fact)::
2340 sage: set_random_seed()
2341 sage: n = ZZ.random_element(2)
2342 sage: Q = QuaternionAlgebra(QQ,-1,-1)
2343 sage: X = random_matrix(Q, n)
2344 sage: Y = random_matrix(Q, n)
2345 sage: Xe = QuaternionMatrixEJA.real_embed(X)
2346 sage: Ye = QuaternionMatrixEJA.real_embed(Y)
2347 sage: XYe = QuaternionMatrixEJA.real_embed(X*Y)
2352 super().real_embed(M
)
2353 quaternions
= M
.base_ring()
2356 F
= QuadraticField(-1, 'I')
2361 t
= z
.coefficient_tuple()
2366 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
2367 [-c
+ d
*i
, a
- b
*i
]])
2368 realM
= ComplexMatrixEJA
.real_embed(cplxM
)
2369 blocks
.append(realM
)
2371 # We should have real entries by now, so use the realest field
2372 # we've got for the return value.
2373 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
2378 def real_unembed(cls
,M
):
2380 The inverse of _embed_quaternion_matrix().
2384 sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA
2388 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
2389 ....: [-2, 1, -4, 3],
2390 ....: [-3, 4, 1, -2],
2391 ....: [-4, -3, 2, 1]])
2392 sage: QuaternionMatrixEJA.real_unembed(M)
2393 [1 + 2*i + 3*j + 4*k]
2397 Unembedding is the inverse of embedding::
2399 sage: set_random_seed()
2400 sage: Q = QuaternionAlgebra(QQ, -1, -1)
2401 sage: M = random_matrix(Q, 3)
2402 sage: Me = QuaternionMatrixEJA.real_embed(M)
2403 sage: QuaternionMatrixEJA.real_unembed(Me) == M
2407 super().real_unembed(M
)
2409 d
= cls
.dimension_over_reals()
2411 # Use the base ring of the matrix to ensure that its entries can be
2412 # multiplied by elements of the quaternion algebra.
2413 Q
= cls
.quaternion_extension(M
.base_ring())
2416 # Go top-left to bottom-right (reading order), converting every
2417 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
2420 for l
in range(n
/d
):
2421 for m
in range(n
/d
):
2422 submat
= ComplexMatrixEJA
.real_unembed(
2423 M
[d
*l
:d
*l
+d
,d
*m
:d
*m
+d
] )
2424 if submat
[0,0] != submat
[1,1].conjugate():
2425 raise ValueError('bad on-diagonal submatrix')
2426 if submat
[0,1] != -submat
[1,0].conjugate():
2427 raise ValueError('bad off-diagonal submatrix')
2428 z
= submat
[0,0].real()
2429 z
+= submat
[0,0].imag()*i
2430 z
+= submat
[0,1].real()*j
2431 z
+= submat
[0,1].imag()*k
2434 return matrix(Q
, n
/d
, elements
)
2437 class QuaternionHermitianEJA(ConcreteEJA
, QuaternionMatrixEJA
):
2439 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
2440 matrices, the usual symmetric Jordan product, and the
2441 real-part-of-trace inner product. It has dimension `2n^2 - n` over
2446 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2450 In theory, our "field" can be any subfield of the reals::
2452 sage: QuaternionHermitianEJA(2, field=RDF, check_axioms=True)
2453 Euclidean Jordan algebra of dimension 6 over Real Double Field
2454 sage: QuaternionHermitianEJA(2, field=RR, check_axioms=True)
2455 Euclidean Jordan algebra of dimension 6 over Real Field with
2456 53 bits of precision
2460 The dimension of this algebra is `2*n^2 - n`::
2462 sage: set_random_seed()
2463 sage: n_max = QuaternionHermitianEJA._max_random_instance_size()
2464 sage: n = ZZ.random_element(1, n_max)
2465 sage: J = QuaternionHermitianEJA(n)
2466 sage: J.dimension() == 2*(n^2) - n
2469 The Jordan multiplication is what we think it is::
2471 sage: set_random_seed()
2472 sage: J = QuaternionHermitianEJA.random_instance()
2473 sage: x,y = J.random_elements(2)
2474 sage: actual = (x*y).to_matrix()
2475 sage: X = x.to_matrix()
2476 sage: Y = y.to_matrix()
2477 sage: expected = (X*Y + Y*X)/2
2478 sage: actual == expected
2480 sage: J(expected) == x*y
2483 We can change the generator prefix::
2485 sage: QuaternionHermitianEJA(2, prefix='a').gens()
2486 (a0, a1, a2, a3, a4, a5)
2488 We can construct the (trivial) algebra of rank zero::
2490 sage: QuaternionHermitianEJA(0)
2491 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2495 def _denormalized_basis(cls
, n
):
2497 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
2499 Why do we embed these? Basically, because all of numerical
2500 linear algebra assumes that you're working with vectors consisting
2501 of `n` entries from a field and scalars from the same field. There's
2502 no way to tell SageMath that (for example) the vectors contain
2503 complex numbers, while the scalar field is real.
2507 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2511 sage: set_random_seed()
2512 sage: n = ZZ.random_element(1,5)
2513 sage: B = QuaternionHermitianEJA._denormalized_basis(n)
2514 sage: all( M.is_symmetric() for M in B )
2519 Q
= QuaternionAlgebra(QQ
,-1,-1)
2522 # This is like the symmetric case, but we need to be careful:
2524 # * We want conjugate-symmetry, not just symmetry.
2525 # * The diagonal will (as a result) be real.
2528 Eij
= matrix
.zero(Q
,n
)
2530 for j
in range(i
+1):
2534 Sij
= cls
.real_embed(Eij
)
2537 # The second, third, and fourth ones have a minus
2538 # because they're conjugated.
2539 # Eij = Eij + Eij.transpose()
2541 Sij_real
= cls
.real_embed(Eij
)
2543 # Eij = I*(Eij - Eij.transpose())
2546 Sij_I
= cls
.real_embed(Eij
)
2548 # Eij = J*(Eij - Eij.transpose())
2551 Sij_J
= cls
.real_embed(Eij
)
2553 # Eij = K*(Eij - Eij.transpose())
2556 Sij_K
= cls
.real_embed(Eij
)
2562 # Since we embedded these, we can drop back to the "field" that we
2563 # started with instead of the quaternion algebra "Q".
2564 return tuple( s
.change_ring(field
) for s
in S
)
2567 def __init__(self
, n
, **kwargs
):
2568 # We know this is a valid EJA, but will double-check
2569 # if the user passes check_axioms=True.
2570 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2576 super().__init
__(self
._denormalized
_basis
(n
),
2577 self
.jordan_product
,
2578 self
.trace_inner_product
,
2579 associative
=associative
,
2582 # TODO: this could be factored out somehow, but is left here
2583 # because the MatrixEJA is not presently a subclass of the
2584 # FDEJA class that defines rank() and one().
2585 self
.rank
.set_cache(n
)
2586 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
2587 self
.one
.set_cache(self(idV
))
2591 def _max_random_instance_size():
2593 The maximum rank of a random QuaternionHermitianEJA.
2595 return 2 # Dimension 6
2598 def random_instance(cls
, **kwargs
):
2600 Return a random instance of this type of algebra.
2602 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2603 return cls(n
, **kwargs
)
2606 class HadamardEJA(ConcreteEJA
):
2608 Return the Euclidean Jordan Algebra corresponding to the set
2609 `R^n` under the Hadamard product.
2611 Note: this is nothing more than the Cartesian product of ``n``
2612 copies of the spin algebra. Once Cartesian product algebras
2613 are implemented, this can go.
2617 sage: from mjo.eja.eja_algebra import HadamardEJA
2621 This multiplication table can be verified by hand::
2623 sage: J = HadamardEJA(3)
2624 sage: e0,e1,e2 = J.gens()
2640 We can change the generator prefix::
2642 sage: HadamardEJA(3, prefix='r').gens()
2646 def __init__(self
, n
, **kwargs
):
2648 jordan_product
= lambda x
,y
: x
2649 inner_product
= lambda x
,y
: x
2651 def jordan_product(x
,y
):
2653 return P( xi
*yi
for (xi
,yi
) in zip(x
,y
) )
2655 def inner_product(x
,y
):
2658 # New defaults for keyword arguments. Don't orthonormalize
2659 # because our basis is already orthonormal with respect to our
2660 # inner-product. Don't check the axioms, because we know this
2661 # is a valid EJA... but do double-check if the user passes
2662 # check_axioms=True. Note: we DON'T override the "check_field"
2663 # default here, because the user can pass in a field!
2664 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2665 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2667 column_basis
= tuple( b
.column() for b
in FreeModule(ZZ
, n
).basis() )
2668 super().__init
__(column_basis
,
2673 self
.rank
.set_cache(n
)
2676 self
.one
.set_cache( self
.zero() )
2678 self
.one
.set_cache( sum(self
.gens()) )
2681 def _max_random_instance_size():
2683 The maximum dimension of a random HadamardEJA.
2688 def random_instance(cls
, **kwargs
):
2690 Return a random instance of this type of algebra.
2692 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2693 return cls(n
, **kwargs
)
2696 class BilinearFormEJA(ConcreteEJA
):
2698 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2699 with the half-trace inner product and jordan product ``x*y =
2700 (<Bx,y>,y_bar>, x0*y_bar + y0*x_bar)`` where `B = 1 \times B22` is
2701 a symmetric positive-definite "bilinear form" matrix. Its
2702 dimension is the size of `B`, and it has rank two in dimensions
2703 larger than two. It reduces to the ``JordanSpinEJA`` when `B` is
2704 the identity matrix of order ``n``.
2706 We insist that the one-by-one upper-left identity block of `B` be
2707 passed in as well so that we can be passed a matrix of size zero
2708 to construct a trivial algebra.
2712 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
2713 ....: JordanSpinEJA)
2717 When no bilinear form is specified, the identity matrix is used,
2718 and the resulting algebra is the Jordan spin algebra::
2720 sage: B = matrix.identity(AA,3)
2721 sage: J0 = BilinearFormEJA(B)
2722 sage: J1 = JordanSpinEJA(3)
2723 sage: J0.multiplication_table() == J0.multiplication_table()
2726 An error is raised if the matrix `B` does not correspond to a
2727 positive-definite bilinear form::
2729 sage: B = matrix.random(QQ,2,3)
2730 sage: J = BilinearFormEJA(B)
2731 Traceback (most recent call last):
2733 ValueError: bilinear form is not positive-definite
2734 sage: B = matrix.zero(QQ,3)
2735 sage: J = BilinearFormEJA(B)
2736 Traceback (most recent call last):
2738 ValueError: bilinear form is not positive-definite
2742 We can create a zero-dimensional algebra::
2744 sage: B = matrix.identity(AA,0)
2745 sage: J = BilinearFormEJA(B)
2749 We can check the multiplication condition given in the Jordan, von
2750 Neumann, and Wigner paper (and also discussed on my "On the
2751 symmetry..." paper). Note that this relies heavily on the standard
2752 choice of basis, as does anything utilizing the bilinear form
2753 matrix. We opt not to orthonormalize the basis, because if we
2754 did, we would have to normalize the `s_{i}` in a similar manner::
2756 sage: set_random_seed()
2757 sage: n = ZZ.random_element(5)
2758 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2759 sage: B11 = matrix.identity(QQ,1)
2760 sage: B22 = M.transpose()*M
2761 sage: B = block_matrix(2,2,[ [B11,0 ],
2763 sage: J = BilinearFormEJA(B, orthonormalize=False)
2764 sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
2765 sage: V = J.vector_space()
2766 sage: sis = [ J( V([0] + (M.inverse()*ei).list()).column() )
2767 ....: for ei in eis ]
2768 sage: actual = [ sis[i]*sis[j]
2769 ....: for i in range(n-1)
2770 ....: for j in range(n-1) ]
2771 sage: expected = [ J.one() if i == j else J.zero()
2772 ....: for i in range(n-1)
2773 ....: for j in range(n-1) ]
2774 sage: actual == expected
2778 def __init__(self
, B
, **kwargs
):
2779 # The matrix "B" is supplied by the user in most cases,
2780 # so it makes sense to check whether or not its positive-
2781 # definite unless we are specifically asked not to...
2782 if ("check_axioms" not in kwargs
) or kwargs
["check_axioms"]:
2783 if not B
.is_positive_definite():
2784 raise ValueError("bilinear form is not positive-definite")
2786 # However, all of the other data for this EJA is computed
2787 # by us in manner that guarantees the axioms are
2788 # satisfied. So, again, unless we are specifically asked to
2789 # verify things, we'll skip the rest of the checks.
2790 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2792 def inner_product(x
,y
):
2793 return (y
.T
*B
*x
)[0,0]
2795 def jordan_product(x
,y
):
2801 z0
= inner_product(y
,x
)
2802 zbar
= y0
*xbar
+ x0
*ybar
2803 return P([z0
] + zbar
.list())
2806 column_basis
= tuple( b
.column() for b
in FreeModule(ZZ
, n
).basis() )
2808 # TODO: I haven't actually checked this, but it seems legit.
2813 super().__init
__(column_basis
,
2816 associative
=associative
,
2819 # The rank of this algebra is two, unless we're in a
2820 # one-dimensional ambient space (because the rank is bounded
2821 # by the ambient dimension).
2822 self
.rank
.set_cache(min(n
,2))
2825 self
.one
.set_cache( self
.zero() )
2827 self
.one
.set_cache( self
.monomial(0) )
2830 def _max_random_instance_size():
2832 The maximum dimension of a random BilinearFormEJA.
2837 def random_instance(cls
, **kwargs
):
2839 Return a random instance of this algebra.
2841 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2843 B
= matrix
.identity(ZZ
, n
)
2844 return cls(B
, **kwargs
)
2846 B11
= matrix
.identity(ZZ
, 1)
2847 M
= matrix
.random(ZZ
, n
-1)
2848 I
= matrix
.identity(ZZ
, n
-1)
2850 while alpha
.is_zero():
2851 alpha
= ZZ
.random_element().abs()
2852 B22
= M
.transpose()*M
+ alpha
*I
2854 from sage
.matrix
.special
import block_matrix
2855 B
= block_matrix(2,2, [ [B11
, ZZ(0) ],
2858 return cls(B
, **kwargs
)
2861 class JordanSpinEJA(BilinearFormEJA
):
2863 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2864 with the usual inner product and jordan product ``x*y =
2865 (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2870 sage: from mjo.eja.eja_algebra import JordanSpinEJA
2874 This multiplication table can be verified by hand::
2876 sage: J = JordanSpinEJA(4)
2877 sage: e0,e1,e2,e3 = J.gens()
2893 We can change the generator prefix::
2895 sage: JordanSpinEJA(2, prefix='B').gens()
2900 Ensure that we have the usual inner product on `R^n`::
2902 sage: set_random_seed()
2903 sage: J = JordanSpinEJA.random_instance()
2904 sage: x,y = J.random_elements(2)
2905 sage: actual = x.inner_product(y)
2906 sage: expected = x.to_vector().inner_product(y.to_vector())
2907 sage: actual == expected
2911 def __init__(self
, n
, **kwargs
):
2912 # This is a special case of the BilinearFormEJA with the
2913 # identity matrix as its bilinear form.
2914 B
= matrix
.identity(ZZ
, n
)
2916 # Don't orthonormalize because our basis is already
2917 # orthonormal with respect to our inner-product.
2918 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2920 # But also don't pass check_field=False here, because the user
2921 # can pass in a field!
2922 super().__init
__(B
, **kwargs
)
2925 def _max_random_instance_size():
2927 The maximum dimension of a random JordanSpinEJA.
2932 def random_instance(cls
, **kwargs
):
2934 Return a random instance of this type of algebra.
2936 Needed here to override the implementation for ``BilinearFormEJA``.
2938 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2939 return cls(n
, **kwargs
)
2942 class TrivialEJA(ConcreteEJA
):
2944 The trivial Euclidean Jordan algebra consisting of only a zero element.
2948 sage: from mjo.eja.eja_algebra import TrivialEJA
2952 sage: J = TrivialEJA()
2959 sage: 7*J.one()*12*J.one()
2961 sage: J.one().inner_product(J.one())
2963 sage: J.one().norm()
2965 sage: J.one().subalgebra_generated_by()
2966 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2971 def __init__(self
, **kwargs
):
2972 jordan_product
= lambda x
,y
: x
2973 inner_product
= lambda x
,y
: 0
2976 # New defaults for keyword arguments
2977 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2978 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2980 super().__init
__(basis
,
2986 # The rank is zero using my definition, namely the dimension of the
2987 # largest subalgebra generated by any element.
2988 self
.rank
.set_cache(0)
2989 self
.one
.set_cache( self
.zero() )
2992 def random_instance(cls
, **kwargs
):
2993 # We don't take a "size" argument so the superclass method is
2994 # inappropriate for us.
2995 return cls(**kwargs
)
2998 class CartesianProductEJA(FiniteDimensionalEJA
):
3000 The external (orthogonal) direct sum of two or more Euclidean
3001 Jordan algebras. Every Euclidean Jordan algebra decomposes into an
3002 orthogonal direct sum of simple Euclidean Jordan algebras which is
3003 then isometric to a Cartesian product, so no generality is lost by
3004 providing only this construction.
3008 sage: from mjo.eja.eja_algebra import (random_eja,
3009 ....: CartesianProductEJA,
3011 ....: JordanSpinEJA,
3012 ....: RealSymmetricEJA)
3016 The Jordan product is inherited from our factors and implemented by
3017 our CombinatorialFreeModule Cartesian product superclass::
3019 sage: set_random_seed()
3020 sage: J1 = HadamardEJA(2)
3021 sage: J2 = RealSymmetricEJA(2)
3022 sage: J = cartesian_product([J1,J2])
3023 sage: x,y = J.random_elements(2)
3027 The ability to retrieve the original factors is implemented by our
3028 CombinatorialFreeModule Cartesian product superclass::
3030 sage: J1 = HadamardEJA(2, field=QQ)
3031 sage: J2 = JordanSpinEJA(3, field=QQ)
3032 sage: J = cartesian_product([J1,J2])
3033 sage: J.cartesian_factors()
3034 (Euclidean Jordan algebra of dimension 2 over Rational Field,
3035 Euclidean Jordan algebra of dimension 3 over Rational Field)
3037 You can provide more than two factors::
3039 sage: J1 = HadamardEJA(2)
3040 sage: J2 = JordanSpinEJA(3)
3041 sage: J3 = RealSymmetricEJA(3)
3042 sage: cartesian_product([J1,J2,J3])
3043 Euclidean Jordan algebra of dimension 2 over Algebraic Real
3044 Field (+) Euclidean Jordan algebra of dimension 3 over Algebraic
3045 Real Field (+) Euclidean Jordan algebra of dimension 6 over
3046 Algebraic Real Field
3048 Rank is additive on a Cartesian product::
3050 sage: J1 = HadamardEJA(1)
3051 sage: J2 = RealSymmetricEJA(2)
3052 sage: J = cartesian_product([J1,J2])
3053 sage: J1.rank.clear_cache()
3054 sage: J2.rank.clear_cache()
3055 sage: J.rank.clear_cache()
3058 sage: J.rank() == J1.rank() + J2.rank()
3061 The same rank computation works over the rationals, with whatever
3064 sage: J1 = HadamardEJA(1, field=QQ, orthonormalize=False)
3065 sage: J2 = RealSymmetricEJA(2, field=QQ, orthonormalize=False)
3066 sage: J = cartesian_product([J1,J2])
3067 sage: J1.rank.clear_cache()
3068 sage: J2.rank.clear_cache()
3069 sage: J.rank.clear_cache()
3072 sage: J.rank() == J1.rank() + J2.rank()
3075 The product algebra will be associative if and only if all of its
3076 components are associative::
3078 sage: J1 = HadamardEJA(2)
3079 sage: J1.is_associative()
3081 sage: J2 = HadamardEJA(3)
3082 sage: J2.is_associative()
3084 sage: J3 = RealSymmetricEJA(3)
3085 sage: J3.is_associative()
3087 sage: CP1 = cartesian_product([J1,J2])
3088 sage: CP1.is_associative()
3090 sage: CP2 = cartesian_product([J1,J3])
3091 sage: CP2.is_associative()
3094 Cartesian products of Cartesian products work::
3096 sage: J1 = JordanSpinEJA(1)
3097 sage: J2 = JordanSpinEJA(1)
3098 sage: J3 = JordanSpinEJA(1)
3099 sage: J = cartesian_product([J1,cartesian_product([J2,J3])])
3100 sage: J.multiplication_table()
3101 +----++----+----+----+
3102 | * || e0 | e1 | e2 |
3103 +====++====+====+====+
3104 | e0 || e0 | 0 | 0 |
3105 +----++----+----+----+
3106 | e1 || 0 | e1 | 0 |
3107 +----++----+----+----+
3108 | e2 || 0 | 0 | e2 |
3109 +----++----+----+----+
3110 sage: HadamardEJA(3).multiplication_table()
3111 +----++----+----+----+
3112 | * || e0 | e1 | e2 |
3113 +====++====+====+====+
3114 | e0 || e0 | 0 | 0 |
3115 +----++----+----+----+
3116 | e1 || 0 | e1 | 0 |
3117 +----++----+----+----+
3118 | e2 || 0 | 0 | e2 |
3119 +----++----+----+----+
3123 All factors must share the same base field::
3125 sage: J1 = HadamardEJA(2, field=QQ)
3126 sage: J2 = RealSymmetricEJA(2)
3127 sage: CartesianProductEJA((J1,J2))
3128 Traceback (most recent call last):
3130 ValueError: all factors must share the same base field
3132 The cached unit element is the same one that would be computed::
3134 sage: set_random_seed() # long time
3135 sage: J1 = random_eja() # long time
3136 sage: J2 = random_eja() # long time
3137 sage: J = cartesian_product([J1,J2]) # long time
3138 sage: actual = J.one() # long time
3139 sage: J.one.clear_cache() # long time
3140 sage: expected = J.one() # long time
3141 sage: actual == expected # long time
3145 Element
= FiniteDimensionalEJAElement
3148 def __init__(self
, factors
, **kwargs
):
3153 self
._sets
= factors
3155 field
= factors
[0].base_ring()
3156 if not all( J
.base_ring() == field
for J
in factors
):
3157 raise ValueError("all factors must share the same base field")
3159 associative
= all( f
.is_associative() for f
in factors
)
3161 MS
= self
.matrix_space()
3165 for b
in factors
[i
].matrix_basis():
3170 basis
= tuple( MS(b
) for b
in basis
)
3172 # Define jordan/inner products that operate on that matrix_basis.
3173 def jordan_product(x
,y
):
3175 (factors
[i
](x
[i
])*factors
[i
](y
[i
])).to_matrix()
3179 def inner_product(x
, y
):
3181 factors
[i
](x
[i
]).inner_product(factors
[i
](y
[i
]))
3185 # There's no need to check the field since it already came
3186 # from an EJA. Likewise the axioms are guaranteed to be
3187 # satisfied, unless the guy writing this class sucks.
3189 # If you want the basis to be orthonormalized, orthonormalize
3191 FiniteDimensionalEJA
.__init
__(self
,
3196 orthonormalize
=False,
3197 associative
=associative
,
3198 cartesian_product
=True,
3202 ones
= tuple(J
.one().to_matrix() for J
in factors
)
3203 self
.one
.set_cache(self(ones
))
3204 self
.rank
.set_cache(sum(J
.rank() for J
in factors
))
3206 def cartesian_factors(self
):
3207 # Copy/pasted from CombinatorialFreeModule_CartesianProduct.
3210 def cartesian_factor(self
, i
):
3212 Return the ``i``th factor of this algebra.
3214 return self
._sets
[i
]
3217 # Copy/pasted from CombinatorialFreeModule_CartesianProduct.
3218 from sage
.categories
.cartesian_product
import cartesian_product
3219 return cartesian_product
.symbol
.join("%s" % factor
3220 for factor
in self
._sets
)
3222 def matrix_space(self
):
3224 Return the space that our matrix basis lives in as a Cartesian
3229 sage: from mjo.eja.eja_algebra import (HadamardEJA,
3230 ....: RealSymmetricEJA)
3234 sage: J1 = HadamardEJA(1)
3235 sage: J2 = RealSymmetricEJA(2)
3236 sage: J = cartesian_product([J1,J2])
3237 sage: J.matrix_space()
3238 The Cartesian product of (Full MatrixSpace of 1 by 1 dense
3239 matrices over Algebraic Real Field, Full MatrixSpace of 2
3240 by 2 dense matrices over Algebraic Real Field)
3243 from sage
.categories
.cartesian_product
import cartesian_product
3244 return cartesian_product( [J
.matrix_space()
3245 for J
in self
.cartesian_factors()] )
3248 def cartesian_projection(self
, i
):
3252 sage: from mjo.eja.eja_algebra import (random_eja,
3253 ....: JordanSpinEJA,
3255 ....: RealSymmetricEJA,
3256 ....: ComplexHermitianEJA)
3260 The projection morphisms are Euclidean Jordan algebra
3263 sage: J1 = HadamardEJA(2)
3264 sage: J2 = RealSymmetricEJA(2)
3265 sage: J = cartesian_product([J1,J2])
3266 sage: J.cartesian_projection(0)
3267 Linear operator between finite-dimensional Euclidean Jordan
3268 algebras represented by the matrix:
3271 Domain: Euclidean Jordan algebra of dimension 2 over Algebraic
3272 Real Field (+) Euclidean Jordan algebra of dimension 3 over
3273 Algebraic Real Field
3274 Codomain: Euclidean Jordan algebra of dimension 2 over Algebraic
3276 sage: J.cartesian_projection(1)
3277 Linear operator between finite-dimensional Euclidean Jordan
3278 algebras represented by the matrix:
3282 Domain: Euclidean Jordan algebra of dimension 2 over Algebraic
3283 Real Field (+) Euclidean Jordan algebra of dimension 3 over
3284 Algebraic Real Field
3285 Codomain: Euclidean Jordan algebra of dimension 3 over Algebraic
3288 The projections work the way you'd expect on the vector
3289 representation of an element::
3291 sage: J1 = JordanSpinEJA(2)
3292 sage: J2 = ComplexHermitianEJA(2)
3293 sage: J = cartesian_product([J1,J2])
3294 sage: pi_left = J.cartesian_projection(0)
3295 sage: pi_right = J.cartesian_projection(1)
3296 sage: pi_left(J.one()).to_vector()
3298 sage: pi_right(J.one()).to_vector()
3300 sage: J.one().to_vector()
3305 The answer never changes::
3307 sage: set_random_seed()
3308 sage: J1 = random_eja()
3309 sage: J2 = random_eja()
3310 sage: J = cartesian_product([J1,J2])
3311 sage: P0 = J.cartesian_projection(0)
3312 sage: P1 = J.cartesian_projection(0)
3317 offset
= sum( self
.cartesian_factor(k
).dimension()
3319 Ji
= self
.cartesian_factor(i
)
3320 Pi
= self
._module
_morphism
(lambda j
: Ji
.monomial(j
- offset
),
3323 return FiniteDimensionalEJAOperator(self
,Ji
,Pi
.matrix())
3326 def cartesian_embedding(self
, i
):
3330 sage: from mjo.eja.eja_algebra import (random_eja,
3331 ....: JordanSpinEJA,
3333 ....: RealSymmetricEJA)
3337 The embedding morphisms are Euclidean Jordan algebra
3340 sage: J1 = HadamardEJA(2)
3341 sage: J2 = RealSymmetricEJA(2)
3342 sage: J = cartesian_product([J1,J2])
3343 sage: J.cartesian_embedding(0)
3344 Linear operator between finite-dimensional Euclidean Jordan
3345 algebras represented by the matrix:
3351 Domain: Euclidean Jordan algebra of dimension 2 over
3352 Algebraic Real Field
3353 Codomain: Euclidean Jordan algebra of dimension 2 over
3354 Algebraic Real Field (+) Euclidean Jordan algebra of
3355 dimension 3 over Algebraic Real Field
3356 sage: J.cartesian_embedding(1)
3357 Linear operator between finite-dimensional Euclidean Jordan
3358 algebras represented by the matrix:
3364 Domain: Euclidean Jordan algebra of dimension 3 over
3365 Algebraic Real Field
3366 Codomain: Euclidean Jordan algebra of dimension 2 over
3367 Algebraic Real Field (+) Euclidean Jordan algebra of
3368 dimension 3 over Algebraic Real Field
3370 The embeddings work the way you'd expect on the vector
3371 representation of an element::
3373 sage: J1 = JordanSpinEJA(3)
3374 sage: J2 = RealSymmetricEJA(2)
3375 sage: J = cartesian_product([J1,J2])
3376 sage: iota_left = J.cartesian_embedding(0)
3377 sage: iota_right = J.cartesian_embedding(1)
3378 sage: iota_left(J1.zero()) == J.zero()
3380 sage: iota_right(J2.zero()) == J.zero()
3382 sage: J1.one().to_vector()
3384 sage: iota_left(J1.one()).to_vector()
3386 sage: J2.one().to_vector()
3388 sage: iota_right(J2.one()).to_vector()
3390 sage: J.one().to_vector()
3395 The answer never changes::
3397 sage: set_random_seed()
3398 sage: J1 = random_eja()
3399 sage: J2 = random_eja()
3400 sage: J = cartesian_product([J1,J2])
3401 sage: E0 = J.cartesian_embedding(0)
3402 sage: E1 = J.cartesian_embedding(0)
3406 Composing a projection with the corresponding inclusion should
3407 produce the identity map, and mismatching them should produce
3410 sage: set_random_seed()
3411 sage: J1 = random_eja()
3412 sage: J2 = random_eja()
3413 sage: J = cartesian_product([J1,J2])
3414 sage: iota_left = J.cartesian_embedding(0)
3415 sage: iota_right = J.cartesian_embedding(1)
3416 sage: pi_left = J.cartesian_projection(0)
3417 sage: pi_right = J.cartesian_projection(1)
3418 sage: pi_left*iota_left == J1.one().operator()
3420 sage: pi_right*iota_right == J2.one().operator()
3422 sage: (pi_left*iota_right).is_zero()
3424 sage: (pi_right*iota_left).is_zero()
3428 offset
= sum( self
.cartesian_factor(k
).dimension()
3430 Ji
= self
.cartesian_factor(i
)
3431 Ei
= Ji
._module
_morphism
(lambda j
: self
.monomial(j
+ offset
),
3433 return FiniteDimensionalEJAOperator(Ji
,self
,Ei
.matrix())
3437 FiniteDimensionalEJA
.CartesianProduct
= CartesianProductEJA
3439 class RationalBasisCartesianProductEJA(CartesianProductEJA
,
3442 A separate class for products of algebras for which we know a
3447 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
3448 ....: RealSymmetricEJA)
3452 This gives us fast characteristic polynomial computations in
3453 product algebras, too::
3456 sage: J1 = JordanSpinEJA(2)
3457 sage: J2 = RealSymmetricEJA(3)
3458 sage: J = cartesian_product([J1,J2])
3459 sage: J.characteristic_polynomial_of().degree()
3465 def __init__(self
, algebras
, **kwargs
):
3466 CartesianProductEJA
.__init
__(self
, algebras
, **kwargs
)
3468 self
._rational
_algebra
= None
3469 if self
.vector_space().base_field() is not QQ
:
3470 self
._rational
_algebra
= cartesian_product([
3471 r
._rational
_algebra
for r
in algebras
3475 RationalBasisEJA
.CartesianProduct
= RationalBasisCartesianProductEJA
3477 def random_eja(*args
, **kwargs
):
3478 J1
= ConcreteEJA
.random_instance(*args
, **kwargs
)
3480 # This might make Cartesian products appear roughly as often as
3481 # any other ConcreteEJA.
3482 if ZZ
.random_element(len(ConcreteEJA
.__subclasses
__()) + 1) == 0:
3483 # Use random_eja() again so we can get more than two factors.
3484 J2
= random_eja(*args
, **kwargs
)
3485 J
= cartesian_product([J1
,J2
])