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 CombinatorialFreeModule_CartesianProduct
)
69 from sage
.matrix
.constructor
import matrix
70 from sage
.matrix
.matrix_space
import MatrixSpace
71 from sage
.misc
.cachefunc
import cached_method
72 from sage
.misc
.table
import table
73 from sage
.modules
.free_module
import FreeModule
, VectorSpace
74 from sage
.rings
.all
import (ZZ
, QQ
, AA
, QQbar
, RR
, RLF
, CLF
,
77 from mjo
.eja
.eja_element
import FiniteDimensionalEJAElement
78 from mjo
.eja
.eja_operator
import FiniteDimensionalEJAOperator
79 from mjo
.eja
.eja_utils
import _all2list
, _mat2vec
81 class FiniteDimensionalEJA(CombinatorialFreeModule
):
83 A finite-dimensional Euclidean Jordan algebra.
87 - ``basis`` -- a tuple; a tuple of basis elements in "matrix
88 form," which must be the same form as the arguments to
89 ``jordan_product`` and ``inner_product``. In reality, "matrix
90 form" can be either vectors, matrices, or a Cartesian product
91 (ordered tuple) of vectors or matrices. All of these would
92 ideally be vector spaces in sage with no special-casing
93 needed; but in reality we turn vectors into column-matrices
94 and Cartesian products `(a,b)` into column matrices
95 `(a,b)^{T}` after converting `a` and `b` themselves.
97 - ``jordan_product`` -- a function; afunction of two ``basis``
98 elements (in matrix form) that returns their jordan product,
99 also in matrix form; this will be applied to ``basis`` to
100 compute a multiplication table for the algebra.
102 - ``inner_product`` -- a function; a function of two ``basis``
103 elements (in matrix form) that returns their inner
104 product. This will be applied to ``basis`` to compute an
105 inner-product table (basically a matrix) for this algebra.
107 - ``field`` -- a subfield of the reals (default: ``AA``); the scalar
108 field for the algebra.
110 - ``orthonormalize`` -- boolean (default: ``True``); whether or
111 not to orthonormalize the basis. Doing so is expensive and
112 generally rules out using the rationals as your ``field``, but
113 is required for spectral decompositions.
117 sage: from mjo.eja.eja_algebra import random_eja
121 We should compute that an element subalgebra is associative even
122 if we circumvent the element method::
124 sage: set_random_seed()
125 sage: J = random_eja(field=QQ,orthonormalize=False)
126 sage: x = J.random_element()
127 sage: A = x.subalgebra_generated_by(orthonormalize=False)
128 sage: basis = tuple(b.superalgebra_element() for b in A.basis())
129 sage: J.subalgebra(basis, orthonormalize=False).is_associative()
133 Element
= FiniteDimensionalEJAElement
142 cartesian_product
=False,
147 # Keep track of whether or not the matrix basis consists of
148 # tuples, since we need special cases for them damned near
149 # everywhere. This is INDEPENDENT of whether or not the
150 # algebra is a cartesian product, since a subalgebra of a
151 # cartesian product will have a basis of tuples, but will not
152 # in general itself be a cartesian product algebra.
153 self
._matrix
_basis
_is
_cartesian
= False
156 if hasattr(basis
[0], 'cartesian_factors'):
157 self
._matrix
_basis
_is
_cartesian
= True
160 if not field
.is_subring(RR
):
161 # Note: this does return true for the real algebraic
162 # field, the rationals, and any quadratic field where
163 # we've specified a real embedding.
164 raise ValueError("scalar field is not real")
166 # If the basis given to us wasn't over the field that it's
167 # supposed to be over, fix that. Or, you know, crash.
168 if not cartesian_product
:
169 # The field for a cartesian product algebra comes from one
170 # of its factors and is the same for all factors, so
171 # there's no need to "reapply" it on product algebras.
172 if self
._matrix
_basis
_is
_cartesian
:
173 # OK since if n == 0, the basis does not consist of tuples.
174 P
= basis
[0].parent()
175 basis
= tuple( P(tuple(b_i
.change_ring(field
) for b_i
in b
))
178 basis
= tuple( b
.change_ring(field
) for b
in basis
)
182 # Check commutativity of the Jordan and inner-products.
183 # This has to be done before we build the multiplication
184 # and inner-product tables/matrices, because we take
185 # advantage of symmetry in the process.
186 if not all( jordan_product(bi
,bj
) == jordan_product(bj
,bi
)
189 raise ValueError("Jordan product is not commutative")
191 if not all( inner_product(bi
,bj
) == inner_product(bj
,bi
)
194 raise ValueError("inner-product is not commutative")
197 category
= MagmaticAlgebras(field
).FiniteDimensional()
198 category
= category
.WithBasis().Unital().Commutative()
200 if associative
is None:
201 # We should figure it out. As with check_axioms, we have to do
202 # this without the help of the _jordan_product_is_associative()
203 # method because we need to know the category before we
204 # initialize the algebra.
205 associative
= all( jordan_product(jordan_product(bi
,bj
),bk
)
207 jordan_product(bi
,jordan_product(bj
,bk
))
213 # Element subalgebras can take advantage of this.
214 category
= category
.Associative()
215 if cartesian_product
:
216 # Use join() here because otherwise we only get the
217 # "Cartesian product of..." and not the things themselves.
218 category
= category
.join([category
,
219 category
.CartesianProducts()])
221 # Call the superclass constructor so that we can use its from_vector()
222 # method to build our multiplication table.
223 CombinatorialFreeModule
.__init
__(self
,
230 # Now comes all of the hard work. We'll be constructing an
231 # ambient vector space V that our (vectorized) basis lives in,
232 # as well as a subspace W of V spanned by those (vectorized)
233 # basis elements. The W-coordinates are the coefficients that
234 # we see in things like x = 1*e1 + 2*e2.
239 degree
= len(_all2list(basis
[0]))
241 # Build an ambient space that fits our matrix basis when
242 # written out as "long vectors."
243 V
= VectorSpace(field
, degree
)
245 # The matrix that will hole the orthonormal -> unorthonormal
246 # coordinate transformation.
247 self
._deortho
_matrix
= None
250 # Save a copy of the un-orthonormalized basis for later.
251 # Convert it to ambient V (vector) coordinates while we're
252 # at it, because we'd have to do it later anyway.
253 deortho_vector_basis
= tuple( V(_all2list(b
)) for b
in basis
)
255 from mjo
.eja
.eja_utils
import gram_schmidt
256 basis
= tuple(gram_schmidt(basis
, inner_product
))
258 # Save the (possibly orthonormalized) matrix basis for
260 self
._matrix
_basis
= basis
262 # Now create the vector space for the algebra, which will have
263 # its own set of non-ambient coordinates (in terms of the
265 vector_basis
= tuple( V(_all2list(b
)) for b
in basis
)
266 W
= V
.span_of_basis( vector_basis
, check
=check_axioms
)
269 # Now "W" is the vector space of our algebra coordinates. The
270 # variables "X1", "X2",... refer to the entries of vectors in
271 # W. Thus to convert back and forth between the orthonormal
272 # coordinates and the given ones, we need to stick the original
274 U
= V
.span_of_basis( deortho_vector_basis
, check
=check_axioms
)
275 self
._deortho
_matrix
= matrix( U
.coordinate_vector(q
)
276 for q
in vector_basis
)
279 # Now we actually compute the multiplication and inner-product
280 # tables/matrices using the possibly-orthonormalized basis.
281 self
._inner
_product
_matrix
= matrix
.identity(field
, n
)
282 self
._multiplication
_table
= [ [0 for j
in range(i
+1)]
285 # Note: the Jordan and inner-products are defined in terms
286 # of the ambient basis. It's important that their arguments
287 # are in ambient coordinates as well.
290 # ortho basis w.r.t. ambient coords
294 # The jordan product returns a matrixy answer, so we
295 # have to convert it to the algebra coordinates.
296 elt
= jordan_product(q_i
, q_j
)
297 elt
= W
.coordinate_vector(V(_all2list(elt
)))
298 self
._multiplication
_table
[i
][j
] = self
.from_vector(elt
)
300 if not orthonormalize
:
301 # If we're orthonormalizing the basis with respect
302 # to an inner-product, then the inner-product
303 # matrix with respect to the resulting basis is
304 # just going to be the identity.
305 ip
= inner_product(q_i
, q_j
)
306 self
._inner
_product
_matrix
[i
,j
] = ip
307 self
._inner
_product
_matrix
[j
,i
] = ip
309 self
._inner
_product
_matrix
._cache
= {'hermitian': True}
310 self
._inner
_product
_matrix
.set_immutable()
313 if not self
._is
_jordanian
():
314 raise ValueError("Jordan identity does not hold")
315 if not self
._inner
_product
_is
_associative
():
316 raise ValueError("inner product is not associative")
319 def _coerce_map_from_base_ring(self
):
321 Disable the map from the base ring into the algebra.
323 Performing a nonsense conversion like this automatically
324 is counterpedagogical. The fallback is to try the usual
325 element constructor, which should also fail.
329 sage: from mjo.eja.eja_algebra import random_eja
333 sage: set_random_seed()
334 sage: J = random_eja()
336 Traceback (most recent call last):
338 ValueError: not an element of this algebra
344 def product_on_basis(self
, i
, j
):
346 Returns the Jordan product of the `i` and `j`th basis elements.
348 This completely defines the Jordan product on the algebra, and
349 is used direclty by our superclass machinery to implement
354 sage: from mjo.eja.eja_algebra import random_eja
358 sage: set_random_seed()
359 sage: J = random_eja()
360 sage: n = J.dimension()
363 sage: ei_ej = J.zero()*J.zero()
365 ....: i = ZZ.random_element(n)
366 ....: j = ZZ.random_element(n)
367 ....: ei = J.monomial(i)
368 ....: ej = J.monomial(j)
369 ....: ei_ej = J.product_on_basis(i,j)
374 # We only stored the lower-triangular portion of the
375 # multiplication table.
377 return self
._multiplication
_table
[i
][j
]
379 return self
._multiplication
_table
[j
][i
]
381 def inner_product(self
, x
, y
):
383 The inner product associated with this Euclidean Jordan algebra.
385 Defaults to the trace inner product, but can be overridden by
386 subclasses if they are sure that the necessary properties are
391 sage: from mjo.eja.eja_algebra import (random_eja,
393 ....: BilinearFormEJA)
397 Our inner product is "associative," which means the following for
398 a symmetric bilinear form::
400 sage: set_random_seed()
401 sage: J = random_eja()
402 sage: x,y,z = J.random_elements(3)
403 sage: (x*y).inner_product(z) == y.inner_product(x*z)
408 Ensure that this is the usual inner product for the algebras
411 sage: set_random_seed()
412 sage: J = HadamardEJA.random_instance()
413 sage: x,y = J.random_elements(2)
414 sage: actual = x.inner_product(y)
415 sage: expected = x.to_vector().inner_product(y.to_vector())
416 sage: actual == expected
419 Ensure that this is one-half of the trace inner-product in a
420 BilinearFormEJA that isn't just the reals (when ``n`` isn't
421 one). This is in Faraut and Koranyi, and also my "On the
424 sage: set_random_seed()
425 sage: J = BilinearFormEJA.random_instance()
426 sage: n = J.dimension()
427 sage: x = J.random_element()
428 sage: y = J.random_element()
429 sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
433 B
= self
._inner
_product
_matrix
434 return (B
*x
.to_vector()).inner_product(y
.to_vector())
437 def is_associative(self
):
439 Return whether or not this algebra's Jordan product is associative.
443 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
447 sage: J = ComplexHermitianEJA(3, field=QQ, orthonormalize=False)
448 sage: J.is_associative()
450 sage: x = sum(J.gens())
451 sage: A = x.subalgebra_generated_by(orthonormalize=False)
452 sage: A.is_associative()
456 return "Associative" in self
.category().axioms()
458 def _is_commutative(self
):
460 Whether or not this algebra's multiplication table is commutative.
462 This method should of course always return ``True``, unless
463 this algebra was constructed with ``check_axioms=False`` and
464 passed an invalid multiplication table.
466 return all( x
*y
== y
*x
for x
in self
.gens() for y
in self
.gens() )
468 def _is_jordanian(self
):
470 Whether or not this algebra's multiplication table respects the
471 Jordan identity `(x^{2})(xy) = x(x^{2}y)`.
473 We only check one arrangement of `x` and `y`, so for a
474 ``True`` result to be truly true, you should also check
475 :meth:`_is_commutative`. This method should of course always
476 return ``True``, unless this algebra was constructed with
477 ``check_axioms=False`` and passed an invalid multiplication table.
479 return all( (self
.monomial(i
)**2)*(self
.monomial(i
)*self
.monomial(j
))
481 (self
.monomial(i
))*((self
.monomial(i
)**2)*self
.monomial(j
))
482 for i
in range(self
.dimension())
483 for j
in range(self
.dimension()) )
485 def _jordan_product_is_associative(self
):
487 Return whether or not this algebra's Jordan product is
488 associative; that is, whether or not `x*(y*z) = (x*y)*z`
491 This method should agree with :meth:`is_associative` unless
492 you lied about the value of the ``associative`` parameter
493 when you constructed the algebra.
497 sage: from mjo.eja.eja_algebra import (random_eja,
498 ....: RealSymmetricEJA,
499 ....: ComplexHermitianEJA,
500 ....: QuaternionHermitianEJA)
504 sage: J = RealSymmetricEJA(4, orthonormalize=False)
505 sage: J._jordan_product_is_associative()
507 sage: x = sum(J.gens())
508 sage: A = x.subalgebra_generated_by()
509 sage: A._jordan_product_is_associative()
514 sage: J = ComplexHermitianEJA(2,field=QQ,orthonormalize=False)
515 sage: J._jordan_product_is_associative()
517 sage: x = sum(J.gens())
518 sage: A = x.subalgebra_generated_by(orthonormalize=False)
519 sage: A._jordan_product_is_associative()
524 sage: J = QuaternionHermitianEJA(2)
525 sage: J._jordan_product_is_associative()
527 sage: x = sum(J.gens())
528 sage: A = x.subalgebra_generated_by()
529 sage: A._jordan_product_is_associative()
534 The values we've presupplied to the constructors agree with
537 sage: set_random_seed()
538 sage: J = random_eja()
539 sage: J.is_associative() == J._jordan_product_is_associative()
545 # Used to check whether or not something is zero.
548 # I don't know of any examples that make this magnitude
549 # necessary because I don't know how to make an
550 # associative algebra when the element subalgebra
551 # construction is unreliable (as it is over RDF; we can't
552 # find the degree of an element because we can't compute
553 # the rank of a matrix). But even multiplication of floats
554 # is non-associative, so *some* epsilon is needed... let's
555 # just take the one from _inner_product_is_associative?
558 for i
in range(self
.dimension()):
559 for j
in range(self
.dimension()):
560 for k
in range(self
.dimension()):
564 diff
= (x
*y
)*z
- x
*(y
*z
)
566 if diff
.norm() > epsilon
:
571 def _inner_product_is_associative(self
):
573 Return whether or not this algebra's inner product `B` is
574 associative; that is, whether or not `B(xy,z) = B(x,yz)`.
576 This method should of course always return ``True``, unless
577 this algebra was constructed with ``check_axioms=False`` and
578 passed an invalid Jordan or inner-product.
582 # Used to check whether or not something is zero.
585 # This choice is sufficient to allow the construction of
586 # QuaternionHermitianEJA(2, field=RDF) with check_axioms=True.
589 for i
in range(self
.dimension()):
590 for j
in range(self
.dimension()):
591 for k
in range(self
.dimension()):
595 diff
= (x
*y
).inner_product(z
) - x
.inner_product(y
*z
)
597 if diff
.abs() > epsilon
:
602 def _element_constructor_(self
, elt
):
604 Construct an element of this algebra from its vector or matrix
607 This gets called only after the parent element _call_ method
608 fails to find a coercion for the argument.
612 sage: from mjo.eja.eja_algebra import (random_eja,
615 ....: RealSymmetricEJA)
619 The identity in `S^n` is converted to the identity in the EJA::
621 sage: J = RealSymmetricEJA(3)
622 sage: I = matrix.identity(QQ,3)
623 sage: J(I) == J.one()
626 This skew-symmetric matrix can't be represented in the EJA::
628 sage: J = RealSymmetricEJA(3)
629 sage: A = matrix(QQ,3, lambda i,j: i-j)
631 Traceback (most recent call last):
633 ValueError: not an element of this algebra
635 Tuples work as well, provided that the matrix basis for the
636 algebra consists of them::
638 sage: J1 = HadamardEJA(3)
639 sage: J2 = RealSymmetricEJA(2)
640 sage: J = cartesian_product([J1,J2])
641 sage: J( (J1.matrix_basis()[1], J2.matrix_basis()[2]) )
646 Ensure that we can convert any element back and forth
647 faithfully between its matrix and algebra representations::
649 sage: set_random_seed()
650 sage: J = random_eja()
651 sage: x = J.random_element()
652 sage: J(x.to_matrix()) == x
655 We cannot coerce elements between algebras just because their
656 matrix representations are compatible::
658 sage: J1 = HadamardEJA(3)
659 sage: J2 = JordanSpinEJA(3)
661 Traceback (most recent call last):
663 ValueError: not an element of this algebra
665 Traceback (most recent call last):
667 ValueError: not an element of this algebra
669 msg
= "not an element of this algebra"
670 if elt
in self
.base_ring():
671 # Ensure that no base ring -> algebra coercion is performed
672 # by this method. There's some stupidity in sage that would
673 # otherwise propagate to this method; for example, sage thinks
674 # that the integer 3 belongs to the space of 2-by-2 matrices.
675 raise ValueError(msg
)
678 # Try to convert a vector into a column-matrix...
680 except (AttributeError, TypeError):
681 # and ignore failure, because we weren't really expecting
682 # a vector as an argument anyway.
685 if elt
not in self
.matrix_space():
686 raise ValueError(msg
)
688 # Thanks for nothing! Matrix spaces aren't vector spaces in
689 # Sage, so we have to figure out its matrix-basis coordinates
690 # ourselves. We use the basis space's ring instead of the
691 # element's ring because the basis space might be an algebraic
692 # closure whereas the base ring of the 3-by-3 identity matrix
693 # could be QQ instead of QQbar.
695 # And, we also have to handle Cartesian product bases (when
696 # the matrix basis consists of tuples) here. The "good news"
697 # is that we're already converting everything to long vectors,
698 # and that strategy works for tuples as well.
700 # We pass check=False because the matrix basis is "guaranteed"
701 # to be linearly independent... right? Ha ha.
703 V
= VectorSpace(self
.base_ring(), len(elt
))
704 W
= V
.span_of_basis( (V(_all2list(s
)) for s
in self
.matrix_basis()),
708 coords
= W
.coordinate_vector(V(elt
))
709 except ArithmeticError: # vector is not in free module
710 raise ValueError(msg
)
712 return self
.from_vector(coords
)
716 Return a string representation of ``self``.
720 sage: from mjo.eja.eja_algebra import JordanSpinEJA
724 Ensure that it says what we think it says::
726 sage: JordanSpinEJA(2, field=AA)
727 Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
728 sage: JordanSpinEJA(3, field=RDF)
729 Euclidean Jordan algebra of dimension 3 over Real Double Field
732 fmt
= "Euclidean Jordan algebra of dimension {} over {}"
733 return fmt
.format(self
.dimension(), self
.base_ring())
737 def characteristic_polynomial_of(self
):
739 Return the algebra's "characteristic polynomial of" function,
740 which is itself a multivariate polynomial that, when evaluated
741 at the coordinates of some algebra element, returns that
742 element's characteristic polynomial.
744 The resulting polynomial has `n+1` variables, where `n` is the
745 dimension of this algebra. The first `n` variables correspond to
746 the coordinates of an algebra element: when evaluated at the
747 coordinates of an algebra element with respect to a certain
748 basis, the result is a univariate polynomial (in the one
749 remaining variable ``t``), namely the characteristic polynomial
754 sage: from mjo.eja.eja_algebra import JordanSpinEJA, TrivialEJA
758 The characteristic polynomial in the spin algebra is given in
759 Alizadeh, Example 11.11::
761 sage: J = JordanSpinEJA(3)
762 sage: p = J.characteristic_polynomial_of(); p
763 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
764 sage: xvec = J.one().to_vector()
768 By definition, the characteristic polynomial is a monic
769 degree-zero polynomial in a rank-zero algebra. Note that
770 Cayley-Hamilton is indeed satisfied since the polynomial
771 ``1`` evaluates to the identity element of the algebra on
774 sage: J = TrivialEJA()
775 sage: J.characteristic_polynomial_of()
782 # The list of coefficient polynomials a_0, a_1, a_2, ..., a_(r-1).
783 a
= self
._charpoly
_coefficients
()
785 # We go to a bit of trouble here to reorder the
786 # indeterminates, so that it's easier to evaluate the
787 # characteristic polynomial at x's coordinates and get back
788 # something in terms of t, which is what we want.
789 S
= PolynomialRing(self
.base_ring(),'t')
793 S
= PolynomialRing(S
, R
.variable_names())
796 return (t
**r
+ sum( a
[k
]*(t
**k
) for k
in range(r
) ))
798 def coordinate_polynomial_ring(self
):
800 The multivariate polynomial ring in which this algebra's
801 :meth:`characteristic_polynomial_of` lives.
805 sage: from mjo.eja.eja_algebra import (HadamardEJA,
806 ....: RealSymmetricEJA)
810 sage: J = HadamardEJA(2)
811 sage: J.coordinate_polynomial_ring()
812 Multivariate Polynomial Ring in X1, X2...
813 sage: J = RealSymmetricEJA(3,field=QQ,orthonormalize=False)
814 sage: J.coordinate_polynomial_ring()
815 Multivariate Polynomial Ring in X1, X2, X3, X4, X5, X6...
818 var_names
= tuple( "X%d" % z
for z
in range(1, self
.dimension()+1) )
819 return PolynomialRing(self
.base_ring(), var_names
)
821 def inner_product(self
, x
, y
):
823 The inner product associated with this Euclidean Jordan algebra.
825 Defaults to the trace inner product, but can be overridden by
826 subclasses if they are sure that the necessary properties are
831 sage: from mjo.eja.eja_algebra import (random_eja,
833 ....: BilinearFormEJA)
837 Our inner product is "associative," which means the following for
838 a symmetric bilinear form::
840 sage: set_random_seed()
841 sage: J = random_eja()
842 sage: x,y,z = J.random_elements(3)
843 sage: (x*y).inner_product(z) == y.inner_product(x*z)
848 Ensure that this is the usual inner product for the algebras
851 sage: set_random_seed()
852 sage: J = HadamardEJA.random_instance()
853 sage: x,y = J.random_elements(2)
854 sage: actual = x.inner_product(y)
855 sage: expected = x.to_vector().inner_product(y.to_vector())
856 sage: actual == expected
859 Ensure that this is one-half of the trace inner-product in a
860 BilinearFormEJA that isn't just the reals (when ``n`` isn't
861 one). This is in Faraut and Koranyi, and also my "On the
864 sage: set_random_seed()
865 sage: J = BilinearFormEJA.random_instance()
866 sage: n = J.dimension()
867 sage: x = J.random_element()
868 sage: y = J.random_element()
869 sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
872 B
= self
._inner
_product
_matrix
873 return (B
*x
.to_vector()).inner_product(y
.to_vector())
876 def is_trivial(self
):
878 Return whether or not this algebra is trivial.
880 A trivial algebra contains only the zero element.
884 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
889 sage: J = ComplexHermitianEJA(3)
895 sage: J = TrivialEJA()
900 return self
.dimension() == 0
903 def multiplication_table(self
):
905 Return a visual representation of this algebra's multiplication
906 table (on basis elements).
910 sage: from mjo.eja.eja_algebra import JordanSpinEJA
914 sage: J = JordanSpinEJA(4)
915 sage: J.multiplication_table()
916 +----++----+----+----+----+
917 | * || e0 | e1 | e2 | e3 |
918 +====++====+====+====+====+
919 | e0 || e0 | e1 | e2 | e3 |
920 +----++----+----+----+----+
921 | e1 || e1 | e0 | 0 | 0 |
922 +----++----+----+----+----+
923 | e2 || e2 | 0 | e0 | 0 |
924 +----++----+----+----+----+
925 | e3 || e3 | 0 | 0 | e0 |
926 +----++----+----+----+----+
930 # Prepend the header row.
931 M
= [["*"] + list(self
.gens())]
933 # And to each subsequent row, prepend an entry that belongs to
934 # the left-side "header column."
935 M
+= [ [self
.monomial(i
)] + [ self
.monomial(i
)*self
.monomial(j
)
939 return table(M
, header_row
=True, header_column
=True, frame
=True)
942 def matrix_basis(self
):
944 Return an (often more natural) representation of this algebras
945 basis as an ordered tuple of matrices.
947 Every finite-dimensional Euclidean Jordan Algebra is a, up to
948 Jordan isomorphism, a direct sum of five simple
949 algebras---four of which comprise Hermitian matrices. And the
950 last type of algebra can of course be thought of as `n`-by-`1`
951 column matrices (ambiguusly called column vectors) to avoid
952 special cases. As a result, matrices (and column vectors) are
953 a natural representation format for Euclidean Jordan algebra
956 But, when we construct an algebra from a basis of matrices,
957 those matrix representations are lost in favor of coordinate
958 vectors *with respect to* that basis. We could eventually
959 convert back if we tried hard enough, but having the original
960 representations handy is valuable enough that we simply store
961 them and return them from this method.
963 Why implement this for non-matrix algebras? Avoiding special
964 cases for the :class:`BilinearFormEJA` pays with simplicity in
965 its own right. But mainly, we would like to be able to assume
966 that elements of a :class:`CartesianProductEJA` can be displayed
967 nicely, without having to have special classes for direct sums
968 one of whose components was a matrix algebra.
972 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
973 ....: RealSymmetricEJA)
977 sage: J = RealSymmetricEJA(2)
979 Finite family {0: e0, 1: e1, 2: e2}
980 sage: J.matrix_basis()
982 [1 0] [ 0 0.7071067811865475?] [0 0]
983 [0 0], [0.7071067811865475? 0], [0 1]
988 sage: J = JordanSpinEJA(2)
990 Finite family {0: e0, 1: e1}
991 sage: J.matrix_basis()
997 return self
._matrix
_basis
1000 def matrix_space(self
):
1002 Return the matrix space in which this algebra's elements live, if
1003 we think of them as matrices (including column vectors of the
1006 "By default" this will be an `n`-by-`1` column-matrix space,
1007 except when the algebra is trivial. There it's `n`-by-`n`
1008 (where `n` is zero), to ensure that two elements of the matrix
1009 space (empty matrices) can be multiplied. For algebras of
1010 matrices, this returns the space in which their
1011 real embeddings live.
1015 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
1016 ....: JordanSpinEJA,
1017 ....: QuaternionHermitianEJA,
1022 By default, the matrix representation is just a column-matrix
1023 equivalent to the vector representation::
1025 sage: J = JordanSpinEJA(3)
1026 sage: J.matrix_space()
1027 Full MatrixSpace of 3 by 1 dense matrices over Algebraic
1030 The matrix representation in the trivial algebra is
1031 zero-by-zero instead of the usual `n`-by-one::
1033 sage: J = TrivialEJA()
1034 sage: J.matrix_space()
1035 Full MatrixSpace of 0 by 0 dense matrices over Algebraic
1038 The matrix space for complex/quaternion Hermitian matrix EJA
1039 is the space in which their real-embeddings live, not the
1040 original complex/quaternion matrix space::
1042 sage: J = ComplexHermitianEJA(2,field=QQ,orthonormalize=False)
1043 sage: J.matrix_space()
1044 Full MatrixSpace of 4 by 4 dense matrices over Rational Field
1045 sage: J = QuaternionHermitianEJA(1,field=QQ,orthonormalize=False)
1046 sage: J.matrix_space()
1047 Full MatrixSpace of 4 by 4 dense matrices over Rational Field
1050 if self
.is_trivial():
1051 return MatrixSpace(self
.base_ring(), 0)
1053 return self
.matrix_basis()[0].parent()
1059 Return the unit element of this algebra.
1063 sage: from mjo.eja.eja_algebra import (HadamardEJA,
1068 We can compute unit element in the Hadamard EJA::
1070 sage: J = HadamardEJA(5)
1072 e0 + e1 + e2 + e3 + e4
1074 The unit element in the Hadamard EJA is inherited in the
1075 subalgebras generated by its elements::
1077 sage: J = HadamardEJA(5)
1079 e0 + e1 + e2 + e3 + e4
1080 sage: x = sum(J.gens())
1081 sage: A = x.subalgebra_generated_by(orthonormalize=False)
1084 sage: A.one().superalgebra_element()
1085 e0 + e1 + e2 + e3 + e4
1089 The identity element acts like the identity, regardless of
1090 whether or not we orthonormalize::
1092 sage: set_random_seed()
1093 sage: J = random_eja()
1094 sage: x = J.random_element()
1095 sage: J.one()*x == x and x*J.one() == x
1097 sage: A = x.subalgebra_generated_by()
1098 sage: y = A.random_element()
1099 sage: A.one()*y == y and y*A.one() == y
1104 sage: set_random_seed()
1105 sage: J = random_eja(field=QQ, orthonormalize=False)
1106 sage: x = J.random_element()
1107 sage: J.one()*x == x and x*J.one() == x
1109 sage: A = x.subalgebra_generated_by(orthonormalize=False)
1110 sage: y = A.random_element()
1111 sage: A.one()*y == y and y*A.one() == y
1114 The matrix of the unit element's operator is the identity,
1115 regardless of the base field and whether or not we
1118 sage: set_random_seed()
1119 sage: J = random_eja()
1120 sage: actual = J.one().operator().matrix()
1121 sage: expected = matrix.identity(J.base_ring(), J.dimension())
1122 sage: actual == expected
1124 sage: x = J.random_element()
1125 sage: A = x.subalgebra_generated_by()
1126 sage: actual = A.one().operator().matrix()
1127 sage: expected = matrix.identity(A.base_ring(), A.dimension())
1128 sage: actual == expected
1133 sage: set_random_seed()
1134 sage: J = random_eja(field=QQ, orthonormalize=False)
1135 sage: actual = J.one().operator().matrix()
1136 sage: expected = matrix.identity(J.base_ring(), J.dimension())
1137 sage: actual == expected
1139 sage: x = J.random_element()
1140 sage: A = x.subalgebra_generated_by(orthonormalize=False)
1141 sage: actual = A.one().operator().matrix()
1142 sage: expected = matrix.identity(A.base_ring(), A.dimension())
1143 sage: actual == expected
1146 Ensure that the cached unit element (often precomputed by
1147 hand) agrees with the computed one::
1149 sage: set_random_seed()
1150 sage: J = random_eja()
1151 sage: cached = J.one()
1152 sage: J.one.clear_cache()
1153 sage: J.one() == cached
1158 sage: set_random_seed()
1159 sage: J = random_eja(field=QQ, orthonormalize=False)
1160 sage: cached = J.one()
1161 sage: J.one.clear_cache()
1162 sage: J.one() == cached
1166 # We can brute-force compute the matrices of the operators
1167 # that correspond to the basis elements of this algebra.
1168 # If some linear combination of those basis elements is the
1169 # algebra identity, then the same linear combination of
1170 # their matrices has to be the identity matrix.
1172 # Of course, matrices aren't vectors in sage, so we have to
1173 # appeal to the "long vectors" isometry.
1174 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
1176 # Now we use basic linear algebra to find the coefficients,
1177 # of the matrices-as-vectors-linear-combination, which should
1178 # work for the original algebra basis too.
1179 A
= matrix(self
.base_ring(), oper_vecs
)
1181 # We used the isometry on the left-hand side already, but we
1182 # still need to do it for the right-hand side. Recall that we
1183 # wanted something that summed to the identity matrix.
1184 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
1186 # Now if there's an identity element in the algebra, this
1187 # should work. We solve on the left to avoid having to
1188 # transpose the matrix "A".
1189 return self
.from_vector(A
.solve_left(b
))
1192 def peirce_decomposition(self
, c
):
1194 The Peirce decomposition of this algebra relative to the
1197 In the future, this can be extended to a complete system of
1198 orthogonal idempotents.
1202 - ``c`` -- an idempotent of this algebra.
1206 A triple (J0, J5, J1) containing two subalgebras and one subspace
1209 - ``J0`` -- the algebra on the eigenspace of ``c.operator()``
1210 corresponding to the eigenvalue zero.
1212 - ``J5`` -- the eigenspace (NOT a subalgebra) of ``c.operator()``
1213 corresponding to the eigenvalue one-half.
1215 - ``J1`` -- the algebra on the eigenspace of ``c.operator()``
1216 corresponding to the eigenvalue one.
1218 These are the only possible eigenspaces for that operator, and this
1219 algebra is a direct sum of them. The spaces ``J0`` and ``J1`` are
1220 orthogonal, and are subalgebras of this algebra with the appropriate
1225 sage: from mjo.eja.eja_algebra import random_eja, RealSymmetricEJA
1229 The canonical example comes from the symmetric matrices, which
1230 decompose into diagonal and off-diagonal parts::
1232 sage: J = RealSymmetricEJA(3)
1233 sage: C = matrix(QQ, [ [1,0,0],
1237 sage: J0,J5,J1 = J.peirce_decomposition(c)
1239 Euclidean Jordan algebra of dimension 1...
1241 Vector space of degree 6 and dimension 2...
1243 Euclidean Jordan algebra of dimension 3...
1244 sage: J0.one().to_matrix()
1248 sage: orig_df = AA.options.display_format
1249 sage: AA.options.display_format = 'radical'
1250 sage: J.from_vector(J5.basis()[0]).to_matrix()
1254 sage: J.from_vector(J5.basis()[1]).to_matrix()
1258 sage: AA.options.display_format = orig_df
1259 sage: J1.one().to_matrix()
1266 Every algebra decomposes trivially with respect to its identity
1269 sage: set_random_seed()
1270 sage: J = random_eja()
1271 sage: J0,J5,J1 = J.peirce_decomposition(J.one())
1272 sage: J0.dimension() == 0 and J5.dimension() == 0
1274 sage: J1.superalgebra() == J and J1.dimension() == J.dimension()
1277 The decomposition is into eigenspaces, and its components are
1278 therefore necessarily orthogonal. Moreover, the identity
1279 elements in the two subalgebras are the projections onto their
1280 respective subspaces of the superalgebra's identity element::
1282 sage: set_random_seed()
1283 sage: J = random_eja()
1284 sage: x = J.random_element()
1285 sage: if not J.is_trivial():
1286 ....: while x.is_nilpotent():
1287 ....: x = J.random_element()
1288 sage: c = x.subalgebra_idempotent()
1289 sage: J0,J5,J1 = J.peirce_decomposition(c)
1291 sage: for (w,y,z) in zip(J0.basis(), J5.basis(), J1.basis()):
1292 ....: w = w.superalgebra_element()
1293 ....: y = J.from_vector(y)
1294 ....: z = z.superalgebra_element()
1295 ....: ipsum += w.inner_product(y).abs()
1296 ....: ipsum += w.inner_product(z).abs()
1297 ....: ipsum += y.inner_product(z).abs()
1300 sage: J1(c) == J1.one()
1302 sage: J0(J.one() - c) == J0.one()
1306 if not c
.is_idempotent():
1307 raise ValueError("element is not idempotent: %s" % c
)
1309 # Default these to what they should be if they turn out to be
1310 # trivial, because eigenspaces_left() won't return eigenvalues
1311 # corresponding to trivial spaces (e.g. it returns only the
1312 # eigenspace corresponding to lambda=1 if you take the
1313 # decomposition relative to the identity element).
1314 trivial
= self
.subalgebra(())
1315 J0
= trivial
# eigenvalue zero
1316 J5
= VectorSpace(self
.base_ring(), 0) # eigenvalue one-half
1317 J1
= trivial
# eigenvalue one
1319 for (eigval
, eigspace
) in c
.operator().matrix().right_eigenspaces():
1320 if eigval
== ~
(self
.base_ring()(2)):
1323 gens
= tuple( self
.from_vector(b
) for b
in eigspace
.basis() )
1324 subalg
= self
.subalgebra(gens
, check_axioms
=False)
1330 raise ValueError("unexpected eigenvalue: %s" % eigval
)
1335 def random_element(self
, thorough
=False):
1337 Return a random element of this algebra.
1339 Our algebra superclass method only returns a linear
1340 combination of at most two basis elements. We instead
1341 want the vector space "random element" method that
1342 returns a more diverse selection.
1346 - ``thorough`` -- (boolean; default False) whether or not we
1347 should generate irrational coefficients for the random
1348 element when our base ring is irrational; this slows the
1349 algebra operations to a crawl, but any truly random method
1353 # For a general base ring... maybe we can trust this to do the
1354 # right thing? Unlikely, but.
1355 V
= self
.vector_space()
1356 v
= V
.random_element()
1358 if self
.base_ring() is AA
:
1359 # The "random element" method of the algebraic reals is
1360 # stupid at the moment, and only returns integers between
1361 # -2 and 2, inclusive:
1363 # https://trac.sagemath.org/ticket/30875
1365 # Instead, we implement our own "random vector" method,
1366 # and then coerce that into the algebra. We use the vector
1367 # space degree here instead of the dimension because a
1368 # subalgebra could (for example) be spanned by only two
1369 # vectors, each with five coordinates. We need to
1370 # generate all five coordinates.
1372 v
*= QQbar
.random_element().real()
1374 v
*= QQ
.random_element()
1376 return self
.from_vector(V
.coordinate_vector(v
))
1378 def random_elements(self
, count
, thorough
=False):
1380 Return ``count`` random elements as a tuple.
1384 - ``thorough`` -- (boolean; default False) whether or not we
1385 should generate irrational coefficients for the random
1386 elements when our base ring is irrational; this slows the
1387 algebra operations to a crawl, but any truly random method
1392 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1396 sage: J = JordanSpinEJA(3)
1397 sage: x,y,z = J.random_elements(3)
1398 sage: all( [ x in J, y in J, z in J ])
1400 sage: len( J.random_elements(10) ) == 10
1404 return tuple( self
.random_element(thorough
)
1405 for idx
in range(count
) )
1409 def _charpoly_coefficients(self
):
1411 The `r` polynomial coefficients of the "characteristic polynomial
1416 sage: from mjo.eja.eja_algebra import random_eja
1420 The theory shows that these are all homogeneous polynomials of
1423 sage: set_random_seed()
1424 sage: J = random_eja()
1425 sage: all(p.is_homogeneous() for p in J._charpoly_coefficients())
1429 n
= self
.dimension()
1430 R
= self
.coordinate_polynomial_ring()
1432 F
= R
.fraction_field()
1435 # From a result in my book, these are the entries of the
1436 # basis representation of L_x.
1437 return sum( vars[k
]*self
.monomial(k
).operator().matrix()[i
,j
]
1440 L_x
= matrix(F
, n
, n
, L_x_i_j
)
1443 if self
.rank
.is_in_cache():
1445 # There's no need to pad the system with redundant
1446 # columns if we *know* they'll be redundant.
1449 # Compute an extra power in case the rank is equal to
1450 # the dimension (otherwise, we would stop at x^(r-1)).
1451 x_powers
= [ (L_x
**k
)*self
.one().to_vector()
1452 for k
in range(n
+1) ]
1453 A
= matrix
.column(F
, x_powers
[:n
])
1454 AE
= A
.extended_echelon_form()
1461 # The theory says that only the first "r" coefficients are
1462 # nonzero, and they actually live in the original polynomial
1463 # ring and not the fraction field. We negate them because in
1464 # the actual characteristic polynomial, they get moved to the
1465 # other side where x^r lives. We don't bother to trim A_rref
1466 # down to a square matrix and solve the resulting system,
1467 # because the upper-left r-by-r portion of A_rref is
1468 # guaranteed to be the identity matrix, so e.g.
1470 # A_rref.solve_right(Y)
1472 # would just be returning Y.
1473 return (-E
*b
)[:r
].change_ring(R
)
1478 Return the rank of this EJA.
1480 This is a cached method because we know the rank a priori for
1481 all of the algebras we can construct. Thus we can avoid the
1482 expensive ``_charpoly_coefficients()`` call unless we truly
1483 need to compute the whole characteristic polynomial.
1487 sage: from mjo.eja.eja_algebra import (HadamardEJA,
1488 ....: JordanSpinEJA,
1489 ....: RealSymmetricEJA,
1490 ....: ComplexHermitianEJA,
1491 ....: QuaternionHermitianEJA,
1496 The rank of the Jordan spin algebra is always two::
1498 sage: JordanSpinEJA(2).rank()
1500 sage: JordanSpinEJA(3).rank()
1502 sage: JordanSpinEJA(4).rank()
1505 The rank of the `n`-by-`n` Hermitian real, complex, or
1506 quaternion matrices is `n`::
1508 sage: RealSymmetricEJA(4).rank()
1510 sage: ComplexHermitianEJA(3).rank()
1512 sage: QuaternionHermitianEJA(2).rank()
1517 Ensure that every EJA that we know how to construct has a
1518 positive integer rank, unless the algebra is trivial in
1519 which case its rank will be zero::
1521 sage: set_random_seed()
1522 sage: J = random_eja()
1526 sage: r > 0 or (r == 0 and J.is_trivial())
1529 Ensure that computing the rank actually works, since the ranks
1530 of all simple algebras are known and will be cached by default::
1532 sage: set_random_seed() # long time
1533 sage: J = random_eja() # long time
1534 sage: cached = J.rank() # long time
1535 sage: J.rank.clear_cache() # long time
1536 sage: J.rank() == cached # long time
1540 return len(self
._charpoly
_coefficients
())
1543 def subalgebra(self
, basis
, **kwargs
):
1545 Create a subalgebra of this algebra from the given basis.
1547 from mjo
.eja
.eja_subalgebra
import FiniteDimensionalEJASubalgebra
1548 return FiniteDimensionalEJASubalgebra(self
, basis
, **kwargs
)
1551 def vector_space(self
):
1553 Return the vector space that underlies this algebra.
1557 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1561 sage: J = RealSymmetricEJA(2)
1562 sage: J.vector_space()
1563 Vector space of dimension 3 over...
1566 return self
.zero().to_vector().parent().ambient_vector_space()
1570 class RationalBasisEJA(FiniteDimensionalEJA
):
1572 New class for algebras whose supplied basis elements have all rational entries.
1576 sage: from mjo.eja.eja_algebra import BilinearFormEJA
1580 The supplied basis is orthonormalized by default::
1582 sage: B = matrix(QQ, [[1, 0, 0], [0, 25, -32], [0, -32, 41]])
1583 sage: J = BilinearFormEJA(B)
1584 sage: J.matrix_basis()
1601 # Abuse the check_field parameter to check that the entries of
1602 # out basis (in ambient coordinates) are in the field QQ.
1603 if not all( all(b_i
in QQ
for b_i
in b
.list()) for b
in basis
):
1604 raise TypeError("basis not rational")
1606 super().__init
__(basis
,
1610 check_field
=check_field
,
1613 self
._rational
_algebra
= None
1615 # There's no point in constructing the extra algebra if this
1616 # one is already rational.
1618 # Note: the same Jordan and inner-products work here,
1619 # because they are necessarily defined with respect to
1620 # ambient coordinates and not any particular basis.
1621 self
._rational
_algebra
= FiniteDimensionalEJA(
1626 associative
=self
.is_associative(),
1627 orthonormalize
=False,
1632 def _charpoly_coefficients(self
):
1636 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
1637 ....: JordanSpinEJA)
1641 The base ring of the resulting polynomial coefficients is what
1642 it should be, and not the rationals (unless the algebra was
1643 already over the rationals)::
1645 sage: J = JordanSpinEJA(3)
1646 sage: J._charpoly_coefficients()
1647 (X1^2 - X2^2 - X3^2, -2*X1)
1648 sage: a0 = J._charpoly_coefficients()[0]
1650 Algebraic Real Field
1651 sage: a0.base_ring()
1652 Algebraic Real Field
1655 if self
._rational
_algebra
is None:
1656 # There's no need to construct *another* algebra over the
1657 # rationals if this one is already over the
1658 # rationals. Likewise, if we never orthonormalized our
1659 # basis, we might as well just use the given one.
1660 return super()._charpoly
_coefficients
()
1662 # Do the computation over the rationals. The answer will be
1663 # the same, because all we've done is a change of basis.
1664 # Then, change back from QQ to our real base ring
1665 a
= ( a_i
.change_ring(self
.base_ring())
1666 for a_i
in self
._rational
_algebra
._charpoly
_coefficients
() )
1668 if self
._deortho
_matrix
is None:
1669 # This can happen if our base ring was, say, AA and we
1670 # chose not to (or didn't need to) orthonormalize. It's
1671 # still faster to do the computations over QQ even if
1672 # the numbers in the boxes stay the same.
1675 # Otherwise, convert the coordinate variables back to the
1676 # deorthonormalized ones.
1677 R
= self
.coordinate_polynomial_ring()
1678 from sage
.modules
.free_module_element
import vector
1679 X
= vector(R
, R
.gens())
1680 BX
= self
._deortho
_matrix
*X
1682 subs_dict
= { X[i]: BX[i] for i in range(len(X)) }
1683 return tuple( a_i
.subs(subs_dict
) for a_i
in a
)
1685 class ConcreteEJA(RationalBasisEJA
):
1687 A class for the Euclidean Jordan algebras that we know by name.
1689 These are the Jordan algebras whose basis, multiplication table,
1690 rank, and so on are known a priori. More to the point, they are
1691 the Euclidean Jordan algebras for which we are able to conjure up
1692 a "random instance."
1696 sage: from mjo.eja.eja_algebra import ConcreteEJA
1700 Our basis is normalized with respect to the algebra's inner
1701 product, unless we specify otherwise::
1703 sage: set_random_seed()
1704 sage: J = ConcreteEJA.random_instance()
1705 sage: all( b.norm() == 1 for b in J.gens() )
1708 Since our basis is orthonormal with respect to the algebra's inner
1709 product, and since we know that this algebra is an EJA, any
1710 left-multiplication operator's matrix will be symmetric because
1711 natural->EJA basis representation is an isometry and within the
1712 EJA the operator is self-adjoint by the Jordan axiom::
1714 sage: set_random_seed()
1715 sage: J = ConcreteEJA.random_instance()
1716 sage: x = J.random_element()
1717 sage: x.operator().is_self_adjoint()
1722 def _max_random_instance_size():
1724 Return an integer "size" that is an upper bound on the size of
1725 this algebra when it is used in a random test
1726 case. Unfortunately, the term "size" is ambiguous -- when
1727 dealing with `R^n` under either the Hadamard or Jordan spin
1728 product, the "size" refers to the dimension `n`. When dealing
1729 with a matrix algebra (real symmetric or complex/quaternion
1730 Hermitian), it refers to the size of the matrix, which is far
1731 less than the dimension of the underlying vector space.
1733 This method must be implemented in each subclass.
1735 raise NotImplementedError
1738 def random_instance(cls
, *args
, **kwargs
):
1740 Return a random instance of this type of algebra.
1742 This method should be implemented in each subclass.
1744 from sage
.misc
.prandom
import choice
1745 eja_class
= choice(cls
.__subclasses
__())
1747 # These all bubble up to the RationalBasisEJA superclass
1748 # constructor, so any (kw)args valid there are also valid
1750 return eja_class
.random_instance(*args
, **kwargs
)
1755 def dimension_over_reals():
1757 The dimension of this matrix's base ring over the reals.
1759 The reals are dimension one over themselves, obviously; that's
1760 just `\mathbb{R}^{1}`. Likewise, the complex numbers `a + bi`
1761 have dimension two. Finally, the quaternions have dimension
1762 four over the reals.
1764 This is used to determine the size of the matrix returned from
1765 :meth:`real_embed`, among other things.
1767 raise NotImplementedError
1770 def real_embed(cls
,M
):
1772 Embed the matrix ``M`` into a space of real matrices.
1774 The matrix ``M`` can have entries in any field at the moment:
1775 the real numbers, complex numbers, or quaternions. And although
1776 they are not a field, we can probably support octonions at some
1777 point, too. This function returns a real matrix that "acts like"
1778 the original with respect to matrix multiplication; i.e.
1780 real_embed(M*N) = real_embed(M)*real_embed(N)
1783 if M
.ncols() != M
.nrows():
1784 raise ValueError("the matrix 'M' must be square")
1789 def real_unembed(cls
,M
):
1791 The inverse of :meth:`real_embed`.
1793 if M
.ncols() != M
.nrows():
1794 raise ValueError("the matrix 'M' must be square")
1795 if not ZZ(M
.nrows()).mod(cls
.dimension_over_reals()).is_zero():
1796 raise ValueError("the matrix 'M' must be a real embedding")
1800 def jordan_product(X
,Y
):
1801 return (X
*Y
+ Y
*X
)/2
1804 def trace_inner_product(cls
,X
,Y
):
1806 Compute the trace inner-product of two real-embeddings.
1810 sage: from mjo.eja.eja_algebra import (RealSymmetricEJA,
1811 ....: ComplexHermitianEJA,
1812 ....: QuaternionHermitianEJA)
1816 This gives the same answer as it would if we computed the trace
1817 from the unembedded (original) matrices::
1819 sage: set_random_seed()
1820 sage: J = RealSymmetricEJA.random_instance()
1821 sage: x,y = J.random_elements(2)
1822 sage: Xe = x.to_matrix()
1823 sage: Ye = y.to_matrix()
1824 sage: X = J.real_unembed(Xe)
1825 sage: Y = J.real_unembed(Ye)
1826 sage: expected = (X*Y).trace()
1827 sage: actual = J.trace_inner_product(Xe,Ye)
1828 sage: actual == expected
1833 sage: set_random_seed()
1834 sage: J = ComplexHermitianEJA.random_instance()
1835 sage: x,y = J.random_elements(2)
1836 sage: Xe = x.to_matrix()
1837 sage: Ye = y.to_matrix()
1838 sage: X = J.real_unembed(Xe)
1839 sage: Y = J.real_unembed(Ye)
1840 sage: expected = (X*Y).trace().real()
1841 sage: actual = J.trace_inner_product(Xe,Ye)
1842 sage: actual == expected
1847 sage: set_random_seed()
1848 sage: J = QuaternionHermitianEJA.random_instance()
1849 sage: x,y = J.random_elements(2)
1850 sage: Xe = x.to_matrix()
1851 sage: Ye = y.to_matrix()
1852 sage: X = J.real_unembed(Xe)
1853 sage: Y = J.real_unembed(Ye)
1854 sage: expected = (X*Y).trace().coefficient_tuple()[0]
1855 sage: actual = J.trace_inner_product(Xe,Ye)
1856 sage: actual == expected
1860 Xu
= cls
.real_unembed(X
)
1861 Yu
= cls
.real_unembed(Y
)
1862 tr
= (Xu
*Yu
).trace()
1865 # Works in QQ, AA, RDF, et cetera.
1867 except AttributeError:
1868 # A quaternion doesn't have a real() method, but does
1869 # have coefficient_tuple() method that returns the
1870 # coefficients of 1, i, j, and k -- in that order.
1871 return tr
.coefficient_tuple()[0]
1874 class RealMatrixEJA(MatrixEJA
):
1876 def dimension_over_reals():
1880 class RealSymmetricEJA(ConcreteEJA
, RealMatrixEJA
):
1882 The rank-n simple EJA consisting of real symmetric n-by-n
1883 matrices, the usual symmetric Jordan product, and the trace inner
1884 product. It has dimension `(n^2 + n)/2` over the reals.
1888 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1892 sage: J = RealSymmetricEJA(2)
1893 sage: e0, e1, e2 = J.gens()
1901 In theory, our "field" can be any subfield of the reals::
1903 sage: RealSymmetricEJA(2, field=RDF, check_axioms=True)
1904 Euclidean Jordan algebra of dimension 3 over Real Double Field
1905 sage: RealSymmetricEJA(2, field=RR, check_axioms=True)
1906 Euclidean Jordan algebra of dimension 3 over Real Field with
1907 53 bits of precision
1911 The dimension of this algebra is `(n^2 + n) / 2`::
1913 sage: set_random_seed()
1914 sage: n_max = RealSymmetricEJA._max_random_instance_size()
1915 sage: n = ZZ.random_element(1, n_max)
1916 sage: J = RealSymmetricEJA(n)
1917 sage: J.dimension() == (n^2 + n)/2
1920 The Jordan multiplication is what we think it is::
1922 sage: set_random_seed()
1923 sage: J = RealSymmetricEJA.random_instance()
1924 sage: x,y = J.random_elements(2)
1925 sage: actual = (x*y).to_matrix()
1926 sage: X = x.to_matrix()
1927 sage: Y = y.to_matrix()
1928 sage: expected = (X*Y + Y*X)/2
1929 sage: actual == expected
1931 sage: J(expected) == x*y
1934 We can change the generator prefix::
1936 sage: RealSymmetricEJA(3, prefix='q').gens()
1937 (q0, q1, q2, q3, q4, q5)
1939 We can construct the (trivial) algebra of rank zero::
1941 sage: RealSymmetricEJA(0)
1942 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1946 def _denormalized_basis(cls
, n
):
1948 Return a basis for the space of real symmetric n-by-n matrices.
1952 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1956 sage: set_random_seed()
1957 sage: n = ZZ.random_element(1,5)
1958 sage: B = RealSymmetricEJA._denormalized_basis(n)
1959 sage: all( M.is_symmetric() for M in B)
1963 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1967 for j
in range(i
+1):
1968 Eij
= matrix(ZZ
, n
, lambda k
,l
: k
==i
and l
==j
)
1972 Sij
= Eij
+ Eij
.transpose()
1978 def _max_random_instance_size():
1979 return 4 # Dimension 10
1982 def random_instance(cls
, **kwargs
):
1984 Return a random instance of this type of algebra.
1986 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1987 return cls(n
, **kwargs
)
1989 def __init__(self
, n
, **kwargs
):
1990 # We know this is a valid EJA, but will double-check
1991 # if the user passes check_axioms=True.
1992 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
1998 super().__init
__(self
._denormalized
_basis
(n
),
1999 self
.jordan_product
,
2000 self
.trace_inner_product
,
2001 associative
=associative
,
2004 # TODO: this could be factored out somehow, but is left here
2005 # because the MatrixEJA is not presently a subclass of the
2006 # FDEJA class that defines rank() and one().
2007 self
.rank
.set_cache(n
)
2008 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
2009 self
.one
.set_cache(self(idV
))
2013 class ComplexMatrixEJA(MatrixEJA
):
2014 # A manual dictionary-cache for the complex_extension() method,
2015 # since apparently @classmethods can't also be @cached_methods.
2016 _complex_extension
= {}
2019 def complex_extension(cls
,field
):
2021 The complex field that we embed/unembed, as an extension
2022 of the given ``field``.
2024 if field
in cls
._complex
_extension
:
2025 return cls
._complex
_extension
[field
]
2027 # Sage doesn't know how to adjoin the complex "i" (the root of
2028 # x^2 + 1) to a field in a general way. Here, we just enumerate
2029 # all of the cases that I have cared to support so far.
2031 # Sage doesn't know how to embed AA into QQbar, i.e. how
2032 # to adjoin sqrt(-1) to AA.
2034 elif not field
.is_exact():
2036 F
= field
.complex_field()
2038 # Works for QQ and... maybe some other fields.
2039 R
= PolynomialRing(field
, 'z')
2041 F
= field
.extension(z
**2 + 1, 'I', embedding
=CLF(-1).sqrt())
2043 cls
._complex
_extension
[field
] = F
2047 def dimension_over_reals():
2051 def real_embed(cls
,M
):
2053 Embed the n-by-n complex matrix ``M`` into the space of real
2054 matrices of size 2n-by-2n via the map the sends each entry `z = a +
2055 bi` to the block matrix ``[[a,b],[-b,a]]``.
2059 sage: from mjo.eja.eja_algebra import ComplexMatrixEJA
2063 sage: F = QuadraticField(-1, 'I')
2064 sage: x1 = F(4 - 2*i)
2065 sage: x2 = F(1 + 2*i)
2068 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
2069 sage: ComplexMatrixEJA.real_embed(M)
2078 Embedding is a homomorphism (isomorphism, in fact)::
2080 sage: set_random_seed()
2081 sage: n = ZZ.random_element(3)
2082 sage: F = QuadraticField(-1, 'I')
2083 sage: X = random_matrix(F, n)
2084 sage: Y = random_matrix(F, n)
2085 sage: Xe = ComplexMatrixEJA.real_embed(X)
2086 sage: Ye = ComplexMatrixEJA.real_embed(Y)
2087 sage: XYe = ComplexMatrixEJA.real_embed(X*Y)
2092 super().real_embed(M
)
2095 # We don't need any adjoined elements...
2096 field
= M
.base_ring().base_ring()
2102 blocks
.append(matrix(field
, 2, [ [ a
, b
],
2105 return matrix
.block(field
, n
, blocks
)
2109 def real_unembed(cls
,M
):
2111 The inverse of _embed_complex_matrix().
2115 sage: from mjo.eja.eja_algebra import ComplexMatrixEJA
2119 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
2120 ....: [-2, 1, -4, 3],
2121 ....: [ 9, 10, 11, 12],
2122 ....: [-10, 9, -12, 11] ])
2123 sage: ComplexMatrixEJA.real_unembed(A)
2125 [ 10*I + 9 12*I + 11]
2129 Unembedding is the inverse of embedding::
2131 sage: set_random_seed()
2132 sage: F = QuadraticField(-1, 'I')
2133 sage: M = random_matrix(F, 3)
2134 sage: Me = ComplexMatrixEJA.real_embed(M)
2135 sage: ComplexMatrixEJA.real_unembed(Me) == M
2139 super().real_unembed(M
)
2141 d
= cls
.dimension_over_reals()
2142 F
= cls
.complex_extension(M
.base_ring())
2145 # Go top-left to bottom-right (reading order), converting every
2146 # 2-by-2 block we see to a single complex element.
2148 for k
in range(n
/d
):
2149 for j
in range(n
/d
):
2150 submat
= M
[d
*k
:d
*k
+d
,d
*j
:d
*j
+d
]
2151 if submat
[0,0] != submat
[1,1]:
2152 raise ValueError('bad on-diagonal submatrix')
2153 if submat
[0,1] != -submat
[1,0]:
2154 raise ValueError('bad off-diagonal submatrix')
2155 z
= submat
[0,0] + submat
[0,1]*i
2158 return matrix(F
, n
/d
, elements
)
2161 class ComplexHermitianEJA(ConcreteEJA
, ComplexMatrixEJA
):
2163 The rank-n simple EJA consisting of complex Hermitian n-by-n
2164 matrices over the real numbers, the usual symmetric Jordan product,
2165 and the real-part-of-trace inner product. It has dimension `n^2` over
2170 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
2174 In theory, our "field" can be any subfield of the reals::
2176 sage: ComplexHermitianEJA(2, field=RDF, check_axioms=True)
2177 Euclidean Jordan algebra of dimension 4 over Real Double Field
2178 sage: ComplexHermitianEJA(2, field=RR, check_axioms=True)
2179 Euclidean Jordan algebra of dimension 4 over Real Field with
2180 53 bits of precision
2184 The dimension of this algebra is `n^2`::
2186 sage: set_random_seed()
2187 sage: n_max = ComplexHermitianEJA._max_random_instance_size()
2188 sage: n = ZZ.random_element(1, n_max)
2189 sage: J = ComplexHermitianEJA(n)
2190 sage: J.dimension() == n^2
2193 The Jordan multiplication is what we think it is::
2195 sage: set_random_seed()
2196 sage: J = ComplexHermitianEJA.random_instance()
2197 sage: x,y = J.random_elements(2)
2198 sage: actual = (x*y).to_matrix()
2199 sage: X = x.to_matrix()
2200 sage: Y = y.to_matrix()
2201 sage: expected = (X*Y + Y*X)/2
2202 sage: actual == expected
2204 sage: J(expected) == x*y
2207 We can change the generator prefix::
2209 sage: ComplexHermitianEJA(2, prefix='z').gens()
2212 We can construct the (trivial) algebra of rank zero::
2214 sage: ComplexHermitianEJA(0)
2215 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2220 def _denormalized_basis(cls
, n
):
2222 Returns a basis for the space of complex Hermitian n-by-n matrices.
2224 Why do we embed these? Basically, because all of numerical linear
2225 algebra assumes that you're working with vectors consisting of `n`
2226 entries from a field and scalars from the same field. There's no way
2227 to tell SageMath that (for example) the vectors contain complex
2228 numbers, while the scalar field is real.
2232 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
2236 sage: set_random_seed()
2237 sage: n = ZZ.random_element(1,5)
2238 sage: B = ComplexHermitianEJA._denormalized_basis(n)
2239 sage: all( M.is_symmetric() for M in B)
2244 R
= PolynomialRing(field
, 'z')
2246 F
= field
.extension(z
**2 + 1, 'I')
2249 # This is like the symmetric case, but we need to be careful:
2251 # * We want conjugate-symmetry, not just symmetry.
2252 # * The diagonal will (as a result) be real.
2255 Eij
= matrix
.zero(F
,n
)
2257 for j
in range(i
+1):
2261 Sij
= cls
.real_embed(Eij
)
2264 # The second one has a minus because it's conjugated.
2265 Eij
[j
,i
] = 1 # Eij = Eij + Eij.transpose()
2266 Sij_real
= cls
.real_embed(Eij
)
2268 # Eij = I*Eij - I*Eij.transpose()
2271 Sij_imag
= cls
.real_embed(Eij
)
2277 # Since we embedded these, we can drop back to the "field" that we
2278 # started with instead of the complex extension "F".
2279 return tuple( s
.change_ring(field
) for s
in S
)
2282 def __init__(self
, n
, **kwargs
):
2283 # We know this is a valid EJA, but will double-check
2284 # if the user passes check_axioms=True.
2285 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2291 super().__init
__(self
._denormalized
_basis
(n
),
2292 self
.jordan_product
,
2293 self
.trace_inner_product
,
2294 associative
=associative
,
2296 # TODO: this could be factored out somehow, but is left here
2297 # because the MatrixEJA is not presently a subclass of the
2298 # FDEJA class that defines rank() and one().
2299 self
.rank
.set_cache(n
)
2300 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
2301 self
.one
.set_cache(self(idV
))
2304 def _max_random_instance_size():
2305 return 3 # Dimension 9
2308 def random_instance(cls
, **kwargs
):
2310 Return a random instance of this type of algebra.
2312 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2313 return cls(n
, **kwargs
)
2315 class QuaternionMatrixEJA(MatrixEJA
):
2317 # A manual dictionary-cache for the quaternion_extension() method,
2318 # since apparently @classmethods can't also be @cached_methods.
2319 _quaternion_extension
= {}
2322 def quaternion_extension(cls
,field
):
2324 The quaternion field that we embed/unembed, as an extension
2325 of the given ``field``.
2327 if field
in cls
._quaternion
_extension
:
2328 return cls
._quaternion
_extension
[field
]
2330 Q
= QuaternionAlgebra(field
,-1,-1)
2332 cls
._quaternion
_extension
[field
] = Q
2336 def dimension_over_reals():
2340 def real_embed(cls
,M
):
2342 Embed the n-by-n quaternion matrix ``M`` into the space of real
2343 matrices of size 4n-by-4n by first sending each quaternion entry `z
2344 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
2345 c+di],[-c + di, a-bi]]`, and then embedding those into a real
2350 sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA
2354 sage: Q = QuaternionAlgebra(QQ,-1,-1)
2355 sage: i,j,k = Q.gens()
2356 sage: x = 1 + 2*i + 3*j + 4*k
2357 sage: M = matrix(Q, 1, [[x]])
2358 sage: QuaternionMatrixEJA.real_embed(M)
2364 Embedding is a homomorphism (isomorphism, in fact)::
2366 sage: set_random_seed()
2367 sage: n = ZZ.random_element(2)
2368 sage: Q = QuaternionAlgebra(QQ,-1,-1)
2369 sage: X = random_matrix(Q, n)
2370 sage: Y = random_matrix(Q, n)
2371 sage: Xe = QuaternionMatrixEJA.real_embed(X)
2372 sage: Ye = QuaternionMatrixEJA.real_embed(Y)
2373 sage: XYe = QuaternionMatrixEJA.real_embed(X*Y)
2378 super().real_embed(M
)
2379 quaternions
= M
.base_ring()
2382 F
= QuadraticField(-1, 'I')
2387 t
= z
.coefficient_tuple()
2392 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
2393 [-c
+ d
*i
, a
- b
*i
]])
2394 realM
= ComplexMatrixEJA
.real_embed(cplxM
)
2395 blocks
.append(realM
)
2397 # We should have real entries by now, so use the realest field
2398 # we've got for the return value.
2399 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
2404 def real_unembed(cls
,M
):
2406 The inverse of _embed_quaternion_matrix().
2410 sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA
2414 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
2415 ....: [-2, 1, -4, 3],
2416 ....: [-3, 4, 1, -2],
2417 ....: [-4, -3, 2, 1]])
2418 sage: QuaternionMatrixEJA.real_unembed(M)
2419 [1 + 2*i + 3*j + 4*k]
2423 Unembedding is the inverse of embedding::
2425 sage: set_random_seed()
2426 sage: Q = QuaternionAlgebra(QQ, -1, -1)
2427 sage: M = random_matrix(Q, 3)
2428 sage: Me = QuaternionMatrixEJA.real_embed(M)
2429 sage: QuaternionMatrixEJA.real_unembed(Me) == M
2433 super().real_unembed(M
)
2435 d
= cls
.dimension_over_reals()
2437 # Use the base ring of the matrix to ensure that its entries can be
2438 # multiplied by elements of the quaternion algebra.
2439 Q
= cls
.quaternion_extension(M
.base_ring())
2442 # Go top-left to bottom-right (reading order), converting every
2443 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
2446 for l
in range(n
/d
):
2447 for m
in range(n
/d
):
2448 submat
= ComplexMatrixEJA
.real_unembed(
2449 M
[d
*l
:d
*l
+d
,d
*m
:d
*m
+d
] )
2450 if submat
[0,0] != submat
[1,1].conjugate():
2451 raise ValueError('bad on-diagonal submatrix')
2452 if submat
[0,1] != -submat
[1,0].conjugate():
2453 raise ValueError('bad off-diagonal submatrix')
2454 z
= submat
[0,0].real()
2455 z
+= submat
[0,0].imag()*i
2456 z
+= submat
[0,1].real()*j
2457 z
+= submat
[0,1].imag()*k
2460 return matrix(Q
, n
/d
, elements
)
2463 class QuaternionHermitianEJA(ConcreteEJA
, QuaternionMatrixEJA
):
2465 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
2466 matrices, the usual symmetric Jordan product, and the
2467 real-part-of-trace inner product. It has dimension `2n^2 - n` over
2472 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2476 In theory, our "field" can be any subfield of the reals::
2478 sage: QuaternionHermitianEJA(2, field=RDF, check_axioms=True)
2479 Euclidean Jordan algebra of dimension 6 over Real Double Field
2480 sage: QuaternionHermitianEJA(2, field=RR, check_axioms=True)
2481 Euclidean Jordan algebra of dimension 6 over Real Field with
2482 53 bits of precision
2486 The dimension of this algebra is `2*n^2 - n`::
2488 sage: set_random_seed()
2489 sage: n_max = QuaternionHermitianEJA._max_random_instance_size()
2490 sage: n = ZZ.random_element(1, n_max)
2491 sage: J = QuaternionHermitianEJA(n)
2492 sage: J.dimension() == 2*(n^2) - n
2495 The Jordan multiplication is what we think it is::
2497 sage: set_random_seed()
2498 sage: J = QuaternionHermitianEJA.random_instance()
2499 sage: x,y = J.random_elements(2)
2500 sage: actual = (x*y).to_matrix()
2501 sage: X = x.to_matrix()
2502 sage: Y = y.to_matrix()
2503 sage: expected = (X*Y + Y*X)/2
2504 sage: actual == expected
2506 sage: J(expected) == x*y
2509 We can change the generator prefix::
2511 sage: QuaternionHermitianEJA(2, prefix='a').gens()
2512 (a0, a1, a2, a3, a4, a5)
2514 We can construct the (trivial) algebra of rank zero::
2516 sage: QuaternionHermitianEJA(0)
2517 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2521 def _denormalized_basis(cls
, n
):
2523 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
2525 Why do we embed these? Basically, because all of numerical
2526 linear algebra assumes that you're working with vectors consisting
2527 of `n` entries from a field and scalars from the same field. There's
2528 no way to tell SageMath that (for example) the vectors contain
2529 complex numbers, while the scalar field is real.
2533 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2537 sage: set_random_seed()
2538 sage: n = ZZ.random_element(1,5)
2539 sage: B = QuaternionHermitianEJA._denormalized_basis(n)
2540 sage: all( M.is_symmetric() for M in B )
2545 Q
= QuaternionAlgebra(QQ
,-1,-1)
2548 # This is like the symmetric case, but we need to be careful:
2550 # * We want conjugate-symmetry, not just symmetry.
2551 # * The diagonal will (as a result) be real.
2554 Eij
= matrix
.zero(Q
,n
)
2556 for j
in range(i
+1):
2560 Sij
= cls
.real_embed(Eij
)
2563 # The second, third, and fourth ones have a minus
2564 # because they're conjugated.
2565 # Eij = Eij + Eij.transpose()
2567 Sij_real
= cls
.real_embed(Eij
)
2569 # Eij = I*(Eij - Eij.transpose())
2572 Sij_I
= cls
.real_embed(Eij
)
2574 # Eij = J*(Eij - Eij.transpose())
2577 Sij_J
= cls
.real_embed(Eij
)
2579 # Eij = K*(Eij - Eij.transpose())
2582 Sij_K
= cls
.real_embed(Eij
)
2588 # Since we embedded these, we can drop back to the "field" that we
2589 # started with instead of the quaternion algebra "Q".
2590 return tuple( s
.change_ring(field
) for s
in S
)
2593 def __init__(self
, n
, **kwargs
):
2594 # We know this is a valid EJA, but will double-check
2595 # if the user passes check_axioms=True.
2596 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2602 super().__init
__(self
._denormalized
_basis
(n
),
2603 self
.jordan_product
,
2604 self
.trace_inner_product
,
2605 associative
=associative
,
2608 # TODO: this could be factored out somehow, but is left here
2609 # because the MatrixEJA is not presently a subclass of the
2610 # FDEJA class that defines rank() and one().
2611 self
.rank
.set_cache(n
)
2612 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
2613 self
.one
.set_cache(self(idV
))
2617 def _max_random_instance_size():
2619 The maximum rank of a random QuaternionHermitianEJA.
2621 return 2 # Dimension 6
2624 def random_instance(cls
, **kwargs
):
2626 Return a random instance of this type of algebra.
2628 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2629 return cls(n
, **kwargs
)
2632 class HadamardEJA(ConcreteEJA
):
2634 Return the Euclidean Jordan Algebra corresponding to the set
2635 `R^n` under the Hadamard product.
2637 Note: this is nothing more than the Cartesian product of ``n``
2638 copies of the spin algebra. Once Cartesian product algebras
2639 are implemented, this can go.
2643 sage: from mjo.eja.eja_algebra import HadamardEJA
2647 This multiplication table can be verified by hand::
2649 sage: J = HadamardEJA(3)
2650 sage: e0,e1,e2 = J.gens()
2666 We can change the generator prefix::
2668 sage: HadamardEJA(3, prefix='r').gens()
2672 def __init__(self
, n
, **kwargs
):
2674 jordan_product
= lambda x
,y
: x
2675 inner_product
= lambda x
,y
: x
2677 def jordan_product(x
,y
):
2679 return P( xi
*yi
for (xi
,yi
) in zip(x
,y
) )
2681 def inner_product(x
,y
):
2684 # New defaults for keyword arguments. Don't orthonormalize
2685 # because our basis is already orthonormal with respect to our
2686 # inner-product. Don't check the axioms, because we know this
2687 # is a valid EJA... but do double-check if the user passes
2688 # check_axioms=True. Note: we DON'T override the "check_field"
2689 # default here, because the user can pass in a field!
2690 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2691 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2693 column_basis
= tuple( b
.column() for b
in FreeModule(ZZ
, n
).basis() )
2694 super().__init
__(column_basis
,
2699 self
.rank
.set_cache(n
)
2702 self
.one
.set_cache( self
.zero() )
2704 self
.one
.set_cache( sum(self
.gens()) )
2707 def _max_random_instance_size():
2709 The maximum dimension of a random HadamardEJA.
2714 def random_instance(cls
, **kwargs
):
2716 Return a random instance of this type of algebra.
2718 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2719 return cls(n
, **kwargs
)
2722 class BilinearFormEJA(ConcreteEJA
):
2724 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2725 with the half-trace inner product and jordan product ``x*y =
2726 (<Bx,y>,y_bar>, x0*y_bar + y0*x_bar)`` where `B = 1 \times B22` is
2727 a symmetric positive-definite "bilinear form" matrix. Its
2728 dimension is the size of `B`, and it has rank two in dimensions
2729 larger than two. It reduces to the ``JordanSpinEJA`` when `B` is
2730 the identity matrix of order ``n``.
2732 We insist that the one-by-one upper-left identity block of `B` be
2733 passed in as well so that we can be passed a matrix of size zero
2734 to construct a trivial algebra.
2738 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
2739 ....: JordanSpinEJA)
2743 When no bilinear form is specified, the identity matrix is used,
2744 and the resulting algebra is the Jordan spin algebra::
2746 sage: B = matrix.identity(AA,3)
2747 sage: J0 = BilinearFormEJA(B)
2748 sage: J1 = JordanSpinEJA(3)
2749 sage: J0.multiplication_table() == J0.multiplication_table()
2752 An error is raised if the matrix `B` does not correspond to a
2753 positive-definite bilinear form::
2755 sage: B = matrix.random(QQ,2,3)
2756 sage: J = BilinearFormEJA(B)
2757 Traceback (most recent call last):
2759 ValueError: bilinear form is not positive-definite
2760 sage: B = matrix.zero(QQ,3)
2761 sage: J = BilinearFormEJA(B)
2762 Traceback (most recent call last):
2764 ValueError: bilinear form is not positive-definite
2768 We can create a zero-dimensional algebra::
2770 sage: B = matrix.identity(AA,0)
2771 sage: J = BilinearFormEJA(B)
2775 We can check the multiplication condition given in the Jordan, von
2776 Neumann, and Wigner paper (and also discussed on my "On the
2777 symmetry..." paper). Note that this relies heavily on the standard
2778 choice of basis, as does anything utilizing the bilinear form
2779 matrix. We opt not to orthonormalize the basis, because if we
2780 did, we would have to normalize the `s_{i}` in a similar manner::
2782 sage: set_random_seed()
2783 sage: n = ZZ.random_element(5)
2784 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2785 sage: B11 = matrix.identity(QQ,1)
2786 sage: B22 = M.transpose()*M
2787 sage: B = block_matrix(2,2,[ [B11,0 ],
2789 sage: J = BilinearFormEJA(B, orthonormalize=False)
2790 sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
2791 sage: V = J.vector_space()
2792 sage: sis = [ J( V([0] + (M.inverse()*ei).list()).column() )
2793 ....: for ei in eis ]
2794 sage: actual = [ sis[i]*sis[j]
2795 ....: for i in range(n-1)
2796 ....: for j in range(n-1) ]
2797 sage: expected = [ J.one() if i == j else J.zero()
2798 ....: for i in range(n-1)
2799 ....: for j in range(n-1) ]
2800 sage: actual == expected
2804 def __init__(self
, B
, **kwargs
):
2805 # The matrix "B" is supplied by the user in most cases,
2806 # so it makes sense to check whether or not its positive-
2807 # definite unless we are specifically asked not to...
2808 if ("check_axioms" not in kwargs
) or kwargs
["check_axioms"]:
2809 if not B
.is_positive_definite():
2810 raise ValueError("bilinear form is not positive-definite")
2812 # However, all of the other data for this EJA is computed
2813 # by us in manner that guarantees the axioms are
2814 # satisfied. So, again, unless we are specifically asked to
2815 # verify things, we'll skip the rest of the checks.
2816 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2818 def inner_product(x
,y
):
2819 return (y
.T
*B
*x
)[0,0]
2821 def jordan_product(x
,y
):
2827 z0
= inner_product(y
,x
)
2828 zbar
= y0
*xbar
+ x0
*ybar
2829 return P([z0
] + zbar
.list())
2832 column_basis
= tuple( b
.column() for b
in FreeModule(ZZ
, n
).basis() )
2834 # TODO: I haven't actually checked this, but it seems legit.
2839 super().__init
__(column_basis
,
2842 associative
=associative
,
2845 # The rank of this algebra is two, unless we're in a
2846 # one-dimensional ambient space (because the rank is bounded
2847 # by the ambient dimension).
2848 self
.rank
.set_cache(min(n
,2))
2851 self
.one
.set_cache( self
.zero() )
2853 self
.one
.set_cache( self
.monomial(0) )
2856 def _max_random_instance_size():
2858 The maximum dimension of a random BilinearFormEJA.
2863 def random_instance(cls
, **kwargs
):
2865 Return a random instance of this algebra.
2867 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2869 B
= matrix
.identity(ZZ
, n
)
2870 return cls(B
, **kwargs
)
2872 B11
= matrix
.identity(ZZ
, 1)
2873 M
= matrix
.random(ZZ
, n
-1)
2874 I
= matrix
.identity(ZZ
, n
-1)
2876 while alpha
.is_zero():
2877 alpha
= ZZ
.random_element().abs()
2878 B22
= M
.transpose()*M
+ alpha
*I
2880 from sage
.matrix
.special
import block_matrix
2881 B
= block_matrix(2,2, [ [B11
, ZZ(0) ],
2884 return cls(B
, **kwargs
)
2887 class JordanSpinEJA(BilinearFormEJA
):
2889 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2890 with the usual inner product and jordan product ``x*y =
2891 (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2896 sage: from mjo.eja.eja_algebra import JordanSpinEJA
2900 This multiplication table can be verified by hand::
2902 sage: J = JordanSpinEJA(4)
2903 sage: e0,e1,e2,e3 = J.gens()
2919 We can change the generator prefix::
2921 sage: JordanSpinEJA(2, prefix='B').gens()
2926 Ensure that we have the usual inner product on `R^n`::
2928 sage: set_random_seed()
2929 sage: J = JordanSpinEJA.random_instance()
2930 sage: x,y = J.random_elements(2)
2931 sage: actual = x.inner_product(y)
2932 sage: expected = x.to_vector().inner_product(y.to_vector())
2933 sage: actual == expected
2937 def __init__(self
, n
, **kwargs
):
2938 # This is a special case of the BilinearFormEJA with the
2939 # identity matrix as its bilinear form.
2940 B
= matrix
.identity(ZZ
, n
)
2942 # Don't orthonormalize because our basis is already
2943 # orthonormal with respect to our inner-product.
2944 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2946 # But also don't pass check_field=False here, because the user
2947 # can pass in a field!
2948 super().__init
__(B
, **kwargs
)
2951 def _max_random_instance_size():
2953 The maximum dimension of a random JordanSpinEJA.
2958 def random_instance(cls
, **kwargs
):
2960 Return a random instance of this type of algebra.
2962 Needed here to override the implementation for ``BilinearFormEJA``.
2964 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2965 return cls(n
, **kwargs
)
2968 class TrivialEJA(ConcreteEJA
):
2970 The trivial Euclidean Jordan algebra consisting of only a zero element.
2974 sage: from mjo.eja.eja_algebra import TrivialEJA
2978 sage: J = TrivialEJA()
2985 sage: 7*J.one()*12*J.one()
2987 sage: J.one().inner_product(J.one())
2989 sage: J.one().norm()
2991 sage: J.one().subalgebra_generated_by()
2992 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2997 def __init__(self
, **kwargs
):
2998 jordan_product
= lambda x
,y
: x
2999 inner_product
= lambda x
,y
: 0
3002 # New defaults for keyword arguments
3003 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
3004 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
3006 super().__init
__(basis
,
3012 # The rank is zero using my definition, namely the dimension of the
3013 # largest subalgebra generated by any element.
3014 self
.rank
.set_cache(0)
3015 self
.one
.set_cache( self
.zero() )
3018 def random_instance(cls
, **kwargs
):
3019 # We don't take a "size" argument so the superclass method is
3020 # inappropriate for us.
3021 return cls(**kwargs
)
3024 class CartesianProductEJA(FiniteDimensionalEJA
):
3026 The external (orthogonal) direct sum of two or more Euclidean
3027 Jordan algebras. Every Euclidean Jordan algebra decomposes into an
3028 orthogonal direct sum of simple Euclidean Jordan algebras which is
3029 then isometric to a Cartesian product, so no generality is lost by
3030 providing only this construction.
3034 sage: from mjo.eja.eja_algebra import (random_eja,
3035 ....: CartesianProductEJA,
3037 ....: JordanSpinEJA,
3038 ....: RealSymmetricEJA)
3042 The Jordan product is inherited from our factors and implemented by
3043 our CombinatorialFreeModule Cartesian product superclass::
3045 sage: set_random_seed()
3046 sage: J1 = HadamardEJA(2)
3047 sage: J2 = RealSymmetricEJA(2)
3048 sage: J = cartesian_product([J1,J2])
3049 sage: x,y = J.random_elements(2)
3053 The ability to retrieve the original factors is implemented by our
3054 CombinatorialFreeModule Cartesian product superclass::
3056 sage: J1 = HadamardEJA(2, field=QQ)
3057 sage: J2 = JordanSpinEJA(3, field=QQ)
3058 sage: J = cartesian_product([J1,J2])
3059 sage: J.cartesian_factors()
3060 (Euclidean Jordan algebra of dimension 2 over Rational Field,
3061 Euclidean Jordan algebra of dimension 3 over Rational Field)
3063 You can provide more than two factors::
3065 sage: J1 = HadamardEJA(2)
3066 sage: J2 = JordanSpinEJA(3)
3067 sage: J3 = RealSymmetricEJA(3)
3068 sage: cartesian_product([J1,J2,J3])
3069 Euclidean Jordan algebra of dimension 2 over Algebraic Real
3070 Field (+) Euclidean Jordan algebra of dimension 3 over Algebraic
3071 Real Field (+) Euclidean Jordan algebra of dimension 6 over
3072 Algebraic Real Field
3074 Rank is additive on a Cartesian product::
3076 sage: J1 = HadamardEJA(1)
3077 sage: J2 = RealSymmetricEJA(2)
3078 sage: J = cartesian_product([J1,J2])
3079 sage: J1.rank.clear_cache()
3080 sage: J2.rank.clear_cache()
3081 sage: J.rank.clear_cache()
3084 sage: J.rank() == J1.rank() + J2.rank()
3087 The same rank computation works over the rationals, with whatever
3090 sage: J1 = HadamardEJA(1, field=QQ, orthonormalize=False)
3091 sage: J2 = RealSymmetricEJA(2, field=QQ, orthonormalize=False)
3092 sage: J = cartesian_product([J1,J2])
3093 sage: J1.rank.clear_cache()
3094 sage: J2.rank.clear_cache()
3095 sage: J.rank.clear_cache()
3098 sage: J.rank() == J1.rank() + J2.rank()
3101 The product algebra will be associative if and only if all of its
3102 components are associative::
3104 sage: J1 = HadamardEJA(2)
3105 sage: J1.is_associative()
3107 sage: J2 = HadamardEJA(3)
3108 sage: J2.is_associative()
3110 sage: J3 = RealSymmetricEJA(3)
3111 sage: J3.is_associative()
3113 sage: CP1 = cartesian_product([J1,J2])
3114 sage: CP1.is_associative()
3116 sage: CP2 = cartesian_product([J1,J3])
3117 sage: CP2.is_associative()
3120 Cartesian products of Cartesian products work::
3122 sage: J1 = JordanSpinEJA(1)
3123 sage: J2 = JordanSpinEJA(1)
3124 sage: J3 = JordanSpinEJA(1)
3125 sage: J = cartesian_product([J1,cartesian_product([J2,J3])])
3126 sage: J.multiplication_table()
3127 +----++----+----+----+
3128 | * || e0 | e1 | e2 |
3129 +====++====+====+====+
3130 | e0 || e0 | 0 | 0 |
3131 +----++----+----+----+
3132 | e1 || 0 | e1 | 0 |
3133 +----++----+----+----+
3134 | e2 || 0 | 0 | e2 |
3135 +----++----+----+----+
3136 sage: HadamardEJA(3).multiplication_table()
3137 +----++----+----+----+
3138 | * || e0 | e1 | e2 |
3139 +====++====+====+====+
3140 | e0 || e0 | 0 | 0 |
3141 +----++----+----+----+
3142 | e1 || 0 | e1 | 0 |
3143 +----++----+----+----+
3144 | e2 || 0 | 0 | e2 |
3145 +----++----+----+----+
3149 All factors must share the same base field::
3151 sage: J1 = HadamardEJA(2, field=QQ)
3152 sage: J2 = RealSymmetricEJA(2)
3153 sage: CartesianProductEJA((J1,J2))
3154 Traceback (most recent call last):
3156 ValueError: all factors must share the same base field
3158 The cached unit element is the same one that would be computed::
3160 sage: set_random_seed() # long time
3161 sage: J1 = random_eja() # long time
3162 sage: J2 = random_eja() # long time
3163 sage: J = cartesian_product([J1,J2]) # long time
3164 sage: actual = J.one() # long time
3165 sage: J.one.clear_cache() # long time
3166 sage: expected = J.one() # long time
3167 sage: actual == expected # long time
3171 Element
= FiniteDimensionalEJAElement
3174 def __init__(self
, factors
, **kwargs
):
3179 self
._sets
= factors
3181 field
= factors
[0].base_ring()
3182 if not all( J
.base_ring() == field
for J
in factors
):
3183 raise ValueError("all factors must share the same base field")
3185 associative
= all( f
.is_associative() for f
in factors
)
3187 MS
= self
.matrix_space()
3191 for b
in factors
[i
].matrix_basis():
3196 basis
= tuple( MS(b
) for b
in basis
)
3198 # Define jordan/inner products that operate on that matrix_basis.
3199 def jordan_product(x
,y
):
3201 (factors
[i
](x
[i
])*factors
[i
](y
[i
])).to_matrix()
3205 def inner_product(x
, y
):
3207 factors
[i
](x
[i
]).inner_product(factors
[i
](y
[i
]))
3211 # There's no need to check the field since it already came
3212 # from an EJA. Likewise the axioms are guaranteed to be
3213 # satisfied, unless the guy writing this class sucks.
3215 # If you want the basis to be orthonormalized, orthonormalize
3217 FiniteDimensionalEJA
.__init
__(self
,
3222 orthonormalize
=False,
3223 associative
=associative
,
3224 cartesian_product
=True,
3228 ones
= tuple(J
.one().to_matrix() for J
in factors
)
3229 self
.one
.set_cache(self(ones
))
3230 self
.rank
.set_cache(sum(J
.rank() for J
in factors
))
3232 def cartesian_factors(self
):
3233 # Copy/pasted from CombinatorialFreeModule_CartesianProduct.
3236 def cartesian_factor(self
, i
):
3238 Return the ``i``th factor of this algebra.
3240 return self
._sets
[i
]
3243 # Copy/pasted from CombinatorialFreeModule_CartesianProduct.
3244 from sage
.categories
.cartesian_product
import cartesian_product
3245 return cartesian_product
.symbol
.join("%s" % factor
3246 for factor
in self
._sets
)
3248 def matrix_space(self
):
3250 Return the space that our matrix basis lives in as a Cartesian
3255 sage: from mjo.eja.eja_algebra import (HadamardEJA,
3256 ....: RealSymmetricEJA)
3260 sage: J1 = HadamardEJA(1)
3261 sage: J2 = RealSymmetricEJA(2)
3262 sage: J = cartesian_product([J1,J2])
3263 sage: J.matrix_space()
3264 The Cartesian product of (Full MatrixSpace of 1 by 1 dense
3265 matrices over Algebraic Real Field, Full MatrixSpace of 2
3266 by 2 dense matrices over Algebraic Real Field)
3269 from sage
.categories
.cartesian_product
import cartesian_product
3270 return cartesian_product( [J
.matrix_space()
3271 for J
in self
.cartesian_factors()] )
3274 def cartesian_projection(self
, i
):
3278 sage: from mjo.eja.eja_algebra import (random_eja,
3279 ....: JordanSpinEJA,
3281 ....: RealSymmetricEJA,
3282 ....: ComplexHermitianEJA)
3286 The projection morphisms are Euclidean Jordan algebra
3289 sage: J1 = HadamardEJA(2)
3290 sage: J2 = RealSymmetricEJA(2)
3291 sage: J = cartesian_product([J1,J2])
3292 sage: J.cartesian_projection(0)
3293 Linear operator between finite-dimensional Euclidean Jordan
3294 algebras represented by the matrix:
3297 Domain: Euclidean Jordan algebra of dimension 2 over Algebraic
3298 Real Field (+) Euclidean Jordan algebra of dimension 3 over
3299 Algebraic Real Field
3300 Codomain: Euclidean Jordan algebra of dimension 2 over Algebraic
3302 sage: J.cartesian_projection(1)
3303 Linear operator between finite-dimensional Euclidean Jordan
3304 algebras represented by the matrix:
3308 Domain: Euclidean Jordan algebra of dimension 2 over Algebraic
3309 Real Field (+) Euclidean Jordan algebra of dimension 3 over
3310 Algebraic Real Field
3311 Codomain: Euclidean Jordan algebra of dimension 3 over Algebraic
3314 The projections work the way you'd expect on the vector
3315 representation of an element::
3317 sage: J1 = JordanSpinEJA(2)
3318 sage: J2 = ComplexHermitianEJA(2)
3319 sage: J = cartesian_product([J1,J2])
3320 sage: pi_left = J.cartesian_projection(0)
3321 sage: pi_right = J.cartesian_projection(1)
3322 sage: pi_left(J.one()).to_vector()
3324 sage: pi_right(J.one()).to_vector()
3326 sage: J.one().to_vector()
3331 The answer never changes::
3333 sage: set_random_seed()
3334 sage: J1 = random_eja()
3335 sage: J2 = random_eja()
3336 sage: J = cartesian_product([J1,J2])
3337 sage: P0 = J.cartesian_projection(0)
3338 sage: P1 = J.cartesian_projection(0)
3343 offset
= sum( self
.cartesian_factor(k
).dimension()
3345 Ji
= self
.cartesian_factor(i
)
3346 Pi
= self
._module
_morphism
(lambda j
: Ji
.monomial(j
- offset
),
3349 return FiniteDimensionalEJAOperator(self
,Ji
,Pi
.matrix())
3352 def cartesian_embedding(self
, i
):
3356 sage: from mjo.eja.eja_algebra import (random_eja,
3357 ....: JordanSpinEJA,
3359 ....: RealSymmetricEJA)
3363 The embedding morphisms are Euclidean Jordan algebra
3366 sage: J1 = HadamardEJA(2)
3367 sage: J2 = RealSymmetricEJA(2)
3368 sage: J = cartesian_product([J1,J2])
3369 sage: J.cartesian_embedding(0)
3370 Linear operator between finite-dimensional Euclidean Jordan
3371 algebras represented by the matrix:
3377 Domain: Euclidean Jordan algebra of dimension 2 over
3378 Algebraic Real Field
3379 Codomain: Euclidean Jordan algebra of dimension 2 over
3380 Algebraic Real Field (+) Euclidean Jordan algebra of
3381 dimension 3 over Algebraic Real Field
3382 sage: J.cartesian_embedding(1)
3383 Linear operator between finite-dimensional Euclidean Jordan
3384 algebras represented by the matrix:
3390 Domain: Euclidean Jordan algebra of dimension 3 over
3391 Algebraic Real Field
3392 Codomain: Euclidean Jordan algebra of dimension 2 over
3393 Algebraic Real Field (+) Euclidean Jordan algebra of
3394 dimension 3 over Algebraic Real Field
3396 The embeddings work the way you'd expect on the vector
3397 representation of an element::
3399 sage: J1 = JordanSpinEJA(3)
3400 sage: J2 = RealSymmetricEJA(2)
3401 sage: J = cartesian_product([J1,J2])
3402 sage: iota_left = J.cartesian_embedding(0)
3403 sage: iota_right = J.cartesian_embedding(1)
3404 sage: iota_left(J1.zero()) == J.zero()
3406 sage: iota_right(J2.zero()) == J.zero()
3408 sage: J1.one().to_vector()
3410 sage: iota_left(J1.one()).to_vector()
3412 sage: J2.one().to_vector()
3414 sage: iota_right(J2.one()).to_vector()
3416 sage: J.one().to_vector()
3421 The answer never changes::
3423 sage: set_random_seed()
3424 sage: J1 = random_eja()
3425 sage: J2 = random_eja()
3426 sage: J = cartesian_product([J1,J2])
3427 sage: E0 = J.cartesian_embedding(0)
3428 sage: E1 = J.cartesian_embedding(0)
3432 Composing a projection with the corresponding inclusion should
3433 produce the identity map, and mismatching them should produce
3436 sage: set_random_seed()
3437 sage: J1 = random_eja()
3438 sage: J2 = random_eja()
3439 sage: J = cartesian_product([J1,J2])
3440 sage: iota_left = J.cartesian_embedding(0)
3441 sage: iota_right = J.cartesian_embedding(1)
3442 sage: pi_left = J.cartesian_projection(0)
3443 sage: pi_right = J.cartesian_projection(1)
3444 sage: pi_left*iota_left == J1.one().operator()
3446 sage: pi_right*iota_right == J2.one().operator()
3448 sage: (pi_left*iota_right).is_zero()
3450 sage: (pi_right*iota_left).is_zero()
3454 offset
= sum( self
.cartesian_factor(k
).dimension()
3456 Ji
= self
.cartesian_factor(i
)
3457 Ei
= Ji
._module
_morphism
(lambda j
: self
.monomial(j
+ offset
),
3459 return FiniteDimensionalEJAOperator(Ji
,self
,Ei
.matrix())
3463 FiniteDimensionalEJA
.CartesianProduct
= CartesianProductEJA
3465 class RationalBasisCartesianProductEJA(CartesianProductEJA
,
3468 A separate class for products of algebras for which we know a
3473 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
3474 ....: RealSymmetricEJA)
3478 This gives us fast characteristic polynomial computations in
3479 product algebras, too::
3482 sage: J1 = JordanSpinEJA(2)
3483 sage: J2 = RealSymmetricEJA(3)
3484 sage: J = cartesian_product([J1,J2])
3485 sage: J.characteristic_polynomial_of().degree()
3491 def __init__(self
, algebras
, **kwargs
):
3492 CartesianProductEJA
.__init
__(self
, algebras
, **kwargs
)
3494 self
._rational
_algebra
= None
3495 if self
.vector_space().base_field() is not QQ
:
3496 self
._rational
_algebra
= cartesian_product([
3497 r
._rational
_algebra
for r
in algebras
3501 RationalBasisEJA
.CartesianProduct
= RationalBasisCartesianProductEJA
3503 random_eja
= ConcreteEJA
.random_instance
3505 # def random_eja(*args, **kwargs):
3506 # J1 = ConcreteEJA.random_instance(*args, **kwargs)
3508 # # This might make Cartesian products appear roughly as often as
3509 # # any other ConcreteEJA.
3510 # if ZZ.random_element(len(ConcreteEJA.__subclasses__()) + 1) == 0:
3511 # # Use random_eja() again so we can get more than two factors.
3512 # J2 = random_eja(*args, **kwargs)
3513 # J = cartesian_product([J1,J2])