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,
150 if not field
.is_subring(RR
):
151 # Note: this does return true for the real algebraic
152 # field, the rationals, and any quadratic field where
153 # we've specified a real embedding.
154 raise ValueError("scalar field is not real")
156 from mjo
.eja
.eja_utils
import _change_ring
157 # If the basis given to us wasn't over the field that it's
158 # supposed to be over, fix that. Or, you know, crash.
159 basis
= tuple( _change_ring(b
, field
) for b
in basis
)
162 # Check commutativity of the Jordan and inner-products.
163 # This has to be done before we build the multiplication
164 # and inner-product tables/matrices, because we take
165 # advantage of symmetry in the process.
166 if not all( jordan_product(bi
,bj
) == jordan_product(bj
,bi
)
169 raise ValueError("Jordan product is not commutative")
171 if not all( inner_product(bi
,bj
) == inner_product(bj
,bi
)
174 raise ValueError("inner-product is not commutative")
177 category
= MagmaticAlgebras(field
).FiniteDimensional()
178 category
= category
.WithBasis().Unital().Commutative()
180 if associative
is None:
181 # We should figure it out. As with check_axioms, we have to do
182 # this without the help of the _jordan_product_is_associative()
183 # method because we need to know the category before we
184 # initialize the algebra.
185 associative
= all( jordan_product(jordan_product(bi
,bj
),bk
)
187 jordan_product(bi
,jordan_product(bj
,bk
))
193 # Element subalgebras can take advantage of this.
194 category
= category
.Associative()
195 if cartesian_product
:
196 # Use join() here because otherwise we only get the
197 # "Cartesian product of..." and not the things themselves.
198 category
= category
.join([category
,
199 category
.CartesianProducts()])
201 # Call the superclass constructor so that we can use its from_vector()
202 # method to build our multiplication table.
203 CombinatorialFreeModule
.__init
__(self
,
210 # Now comes all of the hard work. We'll be constructing an
211 # ambient vector space V that our (vectorized) basis lives in,
212 # as well as a subspace W of V spanned by those (vectorized)
213 # basis elements. The W-coordinates are the coefficients that
214 # we see in things like x = 1*e1 + 2*e2.
219 degree
= len(_all2list(basis
[0]))
221 # Build an ambient space that fits our matrix basis when
222 # written out as "long vectors."
223 V
= VectorSpace(field
, degree
)
225 # The matrix that will hole the orthonormal -> unorthonormal
226 # coordinate transformation.
227 self
._deortho
_matrix
= None
230 # Save a copy of the un-orthonormalized basis for later.
231 # Convert it to ambient V (vector) coordinates while we're
232 # at it, because we'd have to do it later anyway.
233 deortho_vector_basis
= tuple( V(_all2list(b
)) for b
in basis
)
235 from mjo
.eja
.eja_utils
import gram_schmidt
236 basis
= tuple(gram_schmidt(basis
, inner_product
))
238 # Save the (possibly orthonormalized) matrix basis for
240 self
._matrix
_basis
= basis
242 # Now create the vector space for the algebra, which will have
243 # its own set of non-ambient coordinates (in terms of the
245 vector_basis
= tuple( V(_all2list(b
)) for b
in basis
)
246 W
= V
.span_of_basis( vector_basis
, check
=check_axioms
)
249 # Now "W" is the vector space of our algebra coordinates. The
250 # variables "X1", "X2",... refer to the entries of vectors in
251 # W. Thus to convert back and forth between the orthonormal
252 # coordinates and the given ones, we need to stick the original
254 U
= V
.span_of_basis( deortho_vector_basis
, check
=check_axioms
)
255 self
._deortho
_matrix
= matrix( U
.coordinate_vector(q
)
256 for q
in vector_basis
)
259 # Now we actually compute the multiplication and inner-product
260 # tables/matrices using the possibly-orthonormalized basis.
261 self
._inner
_product
_matrix
= matrix
.identity(field
, n
)
262 self
._multiplication
_table
= [ [0 for j
in range(i
+1)]
265 # Note: the Jordan and inner-products are defined in terms
266 # of the ambient basis. It's important that their arguments
267 # are in ambient coordinates as well.
270 # ortho basis w.r.t. ambient coords
274 # The jordan product returns a matrixy answer, so we
275 # have to convert it to the algebra coordinates.
276 elt
= jordan_product(q_i
, q_j
)
277 elt
= W
.coordinate_vector(V(_all2list(elt
)))
278 self
._multiplication
_table
[i
][j
] = self
.from_vector(elt
)
280 if not orthonormalize
:
281 # If we're orthonormalizing the basis with respect
282 # to an inner-product, then the inner-product
283 # matrix with respect to the resulting basis is
284 # just going to be the identity.
285 ip
= inner_product(q_i
, q_j
)
286 self
._inner
_product
_matrix
[i
,j
] = ip
287 self
._inner
_product
_matrix
[j
,i
] = ip
289 self
._inner
_product
_matrix
._cache
= {'hermitian': True}
290 self
._inner
_product
_matrix
.set_immutable()
293 if not self
._is
_jordanian
():
294 raise ValueError("Jordan identity does not hold")
295 if not self
._inner
_product
_is
_associative
():
296 raise ValueError("inner product is not associative")
299 def _coerce_map_from_base_ring(self
):
301 Disable the map from the base ring into the algebra.
303 Performing a nonsense conversion like this automatically
304 is counterpedagogical. The fallback is to try the usual
305 element constructor, which should also fail.
309 sage: from mjo.eja.eja_algebra import random_eja
313 sage: set_random_seed()
314 sage: J = random_eja()
316 Traceback (most recent call last):
318 ValueError: not an element of this algebra
324 def product_on_basis(self
, i
, j
):
326 Returns the Jordan product of the `i` and `j`th basis elements.
328 This completely defines the Jordan product on the algebra, and
329 is used direclty by our superclass machinery to implement
334 sage: from mjo.eja.eja_algebra import random_eja
338 sage: set_random_seed()
339 sage: J = random_eja()
340 sage: n = J.dimension()
343 sage: ei_ej = J.zero()*J.zero()
345 ....: i = ZZ.random_element(n)
346 ....: j = ZZ.random_element(n)
347 ....: ei = J.monomial(i)
348 ....: ej = J.monomial(j)
349 ....: ei_ej = J.product_on_basis(i,j)
354 # We only stored the lower-triangular portion of the
355 # multiplication table.
357 return self
._multiplication
_table
[i
][j
]
359 return self
._multiplication
_table
[j
][i
]
361 def inner_product(self
, x
, y
):
363 The inner product associated with this Euclidean Jordan algebra.
365 Defaults to the trace inner product, but can be overridden by
366 subclasses if they are sure that the necessary properties are
371 sage: from mjo.eja.eja_algebra import (random_eja,
373 ....: BilinearFormEJA)
377 Our inner product is "associative," which means the following for
378 a symmetric bilinear form::
380 sage: set_random_seed()
381 sage: J = random_eja()
382 sage: x,y,z = J.random_elements(3)
383 sage: (x*y).inner_product(z) == y.inner_product(x*z)
388 Ensure that this is the usual inner product for the algebras
391 sage: set_random_seed()
392 sage: J = HadamardEJA.random_instance()
393 sage: x,y = J.random_elements(2)
394 sage: actual = x.inner_product(y)
395 sage: expected = x.to_vector().inner_product(y.to_vector())
396 sage: actual == expected
399 Ensure that this is one-half of the trace inner-product in a
400 BilinearFormEJA that isn't just the reals (when ``n`` isn't
401 one). This is in Faraut and Koranyi, and also my "On the
404 sage: set_random_seed()
405 sage: J = BilinearFormEJA.random_instance()
406 sage: n = J.dimension()
407 sage: x = J.random_element()
408 sage: y = J.random_element()
409 sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
413 B
= self
._inner
_product
_matrix
414 return (B
*x
.to_vector()).inner_product(y
.to_vector())
417 def is_associative(self
):
419 Return whether or not this algebra's Jordan product is associative.
423 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
427 sage: J = ComplexHermitianEJA(3, field=QQ, orthonormalize=False)
428 sage: J.is_associative()
430 sage: x = sum(J.gens())
431 sage: A = x.subalgebra_generated_by(orthonormalize=False)
432 sage: A.is_associative()
436 return "Associative" in self
.category().axioms()
438 def _is_commutative(self
):
440 Whether or not this algebra's multiplication table is commutative.
442 This method should of course always return ``True``, unless
443 this algebra was constructed with ``check_axioms=False`` and
444 passed an invalid multiplication table.
446 return all( x
*y
== y
*x
for x
in self
.gens() for y
in self
.gens() )
448 def _is_jordanian(self
):
450 Whether or not this algebra's multiplication table respects the
451 Jordan identity `(x^{2})(xy) = x(x^{2}y)`.
453 We only check one arrangement of `x` and `y`, so for a
454 ``True`` result to be truly true, you should also check
455 :meth:`_is_commutative`. This method should of course always
456 return ``True``, unless this algebra was constructed with
457 ``check_axioms=False`` and passed an invalid multiplication table.
459 return all( (self
.monomial(i
)**2)*(self
.monomial(i
)*self
.monomial(j
))
461 (self
.monomial(i
))*((self
.monomial(i
)**2)*self
.monomial(j
))
462 for i
in range(self
.dimension())
463 for j
in range(self
.dimension()) )
465 def _jordan_product_is_associative(self
):
467 Return whether or not this algebra's Jordan product is
468 associative; that is, whether or not `x*(y*z) = (x*y)*z`
471 This method should agree with :meth:`is_associative` unless
472 you lied about the value of the ``associative`` parameter
473 when you constructed the algebra.
477 sage: from mjo.eja.eja_algebra import (random_eja,
478 ....: RealSymmetricEJA,
479 ....: ComplexHermitianEJA,
480 ....: QuaternionHermitianEJA)
484 sage: J = RealSymmetricEJA(4, orthonormalize=False)
485 sage: J._jordan_product_is_associative()
487 sage: x = sum(J.gens())
488 sage: A = x.subalgebra_generated_by()
489 sage: A._jordan_product_is_associative()
494 sage: J = ComplexHermitianEJA(2,field=QQ,orthonormalize=False)
495 sage: J._jordan_product_is_associative()
497 sage: x = sum(J.gens())
498 sage: A = x.subalgebra_generated_by(orthonormalize=False)
499 sage: A._jordan_product_is_associative()
504 sage: J = QuaternionHermitianEJA(2)
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 The values we've presupplied to the constructors agree with
517 sage: set_random_seed()
518 sage: J = random_eja()
519 sage: J.is_associative() == J._jordan_product_is_associative()
525 # Used to check whether or not something is zero.
528 # I don't know of any examples that make this magnitude
529 # necessary because I don't know how to make an
530 # associative algebra when the element subalgebra
531 # construction is unreliable (as it is over RDF; we can't
532 # find the degree of an element because we can't compute
533 # the rank of a matrix). But even multiplication of floats
534 # is non-associative, so *some* epsilon is needed... let's
535 # just take the one from _inner_product_is_associative?
538 for i
in range(self
.dimension()):
539 for j
in range(self
.dimension()):
540 for k
in range(self
.dimension()):
544 diff
= (x
*y
)*z
- x
*(y
*z
)
546 if diff
.norm() > epsilon
:
551 def _inner_product_is_associative(self
):
553 Return whether or not this algebra's inner product `B` is
554 associative; that is, whether or not `B(xy,z) = B(x,yz)`.
556 This method should of course always return ``True``, unless
557 this algebra was constructed with ``check_axioms=False`` and
558 passed an invalid Jordan or inner-product.
562 # Used to check whether or not something is zero.
565 # This choice is sufficient to allow the construction of
566 # QuaternionHermitianEJA(2, field=RDF) with check_axioms=True.
569 for i
in range(self
.dimension()):
570 for j
in range(self
.dimension()):
571 for k
in range(self
.dimension()):
575 diff
= (x
*y
).inner_product(z
) - x
.inner_product(y
*z
)
577 if diff
.abs() > epsilon
:
582 def _element_constructor_(self
, elt
):
584 Construct an element of this algebra from its vector or matrix
587 This gets called only after the parent element _call_ method
588 fails to find a coercion for the argument.
592 sage: from mjo.eja.eja_algebra import (random_eja,
595 ....: RealSymmetricEJA)
599 The identity in `S^n` is converted to the identity in the EJA::
601 sage: J = RealSymmetricEJA(3)
602 sage: I = matrix.identity(QQ,3)
603 sage: J(I) == J.one()
606 This skew-symmetric matrix can't be represented in the EJA::
608 sage: J = RealSymmetricEJA(3)
609 sage: A = matrix(QQ,3, lambda i,j: i-j)
611 Traceback (most recent call last):
613 ValueError: not an element of this algebra
615 Tuples work as well, provided that the matrix basis for the
616 algebra consists of them::
618 sage: J1 = HadamardEJA(3)
619 sage: J2 = RealSymmetricEJA(2)
620 sage: J = cartesian_product([J1,J2])
621 sage: J( (J1.matrix_basis()[1], J2.matrix_basis()[2]) )
626 Ensure that we can convert any element back and forth
627 faithfully between its matrix and algebra representations::
629 sage: set_random_seed()
630 sage: J = random_eja()
631 sage: x = J.random_element()
632 sage: J(x.to_matrix()) == x
635 We cannot coerce elements between algebras just because their
636 matrix representations are compatible::
638 sage: J1 = HadamardEJA(3)
639 sage: J2 = JordanSpinEJA(3)
641 Traceback (most recent call last):
643 ValueError: not an element of this algebra
645 Traceback (most recent call last):
647 ValueError: not an element of this algebra
649 msg
= "not an element of this algebra"
650 if elt
in self
.base_ring():
651 # Ensure that no base ring -> algebra coercion is performed
652 # by this method. There's some stupidity in sage that would
653 # otherwise propagate to this method; for example, sage thinks
654 # that the integer 3 belongs to the space of 2-by-2 matrices.
655 raise ValueError(msg
)
658 # Try to convert a vector into a column-matrix...
660 except (AttributeError, TypeError):
661 # and ignore failure, because we weren't really expecting
662 # a vector as an argument anyway.
665 if elt
not in self
.matrix_space():
666 raise ValueError(msg
)
668 # Thanks for nothing! Matrix spaces aren't vector spaces in
669 # Sage, so we have to figure out its matrix-basis coordinates
670 # ourselves. We use the basis space's ring instead of the
671 # element's ring because the basis space might be an algebraic
672 # closure whereas the base ring of the 3-by-3 identity matrix
673 # could be QQ instead of QQbar.
675 # And, we also have to handle Cartesian product bases (when
676 # the matrix basis consists of tuples) here. The "good news"
677 # is that we're already converting everything to long vectors,
678 # and that strategy works for tuples as well.
680 # We pass check=False because the matrix basis is "guaranteed"
681 # to be linearly independent... right? Ha ha.
683 V
= VectorSpace(self
.base_ring(), len(elt
))
684 W
= V
.span_of_basis( (V(_all2list(s
)) for s
in self
.matrix_basis()),
688 coords
= W
.coordinate_vector(V(elt
))
689 except ArithmeticError: # vector is not in free module
690 raise ValueError(msg
)
692 return self
.from_vector(coords
)
696 Return a string representation of ``self``.
700 sage: from mjo.eja.eja_algebra import JordanSpinEJA
704 Ensure that it says what we think it says::
706 sage: JordanSpinEJA(2, field=AA)
707 Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
708 sage: JordanSpinEJA(3, field=RDF)
709 Euclidean Jordan algebra of dimension 3 over Real Double Field
712 fmt
= "Euclidean Jordan algebra of dimension {} over {}"
713 return fmt
.format(self
.dimension(), self
.base_ring())
717 def characteristic_polynomial_of(self
):
719 Return the algebra's "characteristic polynomial of" function,
720 which is itself a multivariate polynomial that, when evaluated
721 at the coordinates of some algebra element, returns that
722 element's characteristic polynomial.
724 The resulting polynomial has `n+1` variables, where `n` is the
725 dimension of this algebra. The first `n` variables correspond to
726 the coordinates of an algebra element: when evaluated at the
727 coordinates of an algebra element with respect to a certain
728 basis, the result is a univariate polynomial (in the one
729 remaining variable ``t``), namely the characteristic polynomial
734 sage: from mjo.eja.eja_algebra import JordanSpinEJA, TrivialEJA
738 The characteristic polynomial in the spin algebra is given in
739 Alizadeh, Example 11.11::
741 sage: J = JordanSpinEJA(3)
742 sage: p = J.characteristic_polynomial_of(); p
743 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
744 sage: xvec = J.one().to_vector()
748 By definition, the characteristic polynomial is a monic
749 degree-zero polynomial in a rank-zero algebra. Note that
750 Cayley-Hamilton is indeed satisfied since the polynomial
751 ``1`` evaluates to the identity element of the algebra on
754 sage: J = TrivialEJA()
755 sage: J.characteristic_polynomial_of()
762 # The list of coefficient polynomials a_0, a_1, a_2, ..., a_(r-1).
763 a
= self
._charpoly
_coefficients
()
765 # We go to a bit of trouble here to reorder the
766 # indeterminates, so that it's easier to evaluate the
767 # characteristic polynomial at x's coordinates and get back
768 # something in terms of t, which is what we want.
769 S
= PolynomialRing(self
.base_ring(),'t')
773 S
= PolynomialRing(S
, R
.variable_names())
776 return (t
**r
+ sum( a
[k
]*(t
**k
) for k
in range(r
) ))
778 def coordinate_polynomial_ring(self
):
780 The multivariate polynomial ring in which this algebra's
781 :meth:`characteristic_polynomial_of` lives.
785 sage: from mjo.eja.eja_algebra import (HadamardEJA,
786 ....: RealSymmetricEJA)
790 sage: J = HadamardEJA(2)
791 sage: J.coordinate_polynomial_ring()
792 Multivariate Polynomial Ring in X1, X2...
793 sage: J = RealSymmetricEJA(3,field=QQ,orthonormalize=False)
794 sage: J.coordinate_polynomial_ring()
795 Multivariate Polynomial Ring in X1, X2, X3, X4, X5, X6...
798 var_names
= tuple( "X%d" % z
for z
in range(1, self
.dimension()+1) )
799 return PolynomialRing(self
.base_ring(), var_names
)
801 def inner_product(self
, x
, y
):
803 The inner product associated with this Euclidean Jordan algebra.
805 Defaults to the trace inner product, but can be overridden by
806 subclasses if they are sure that the necessary properties are
811 sage: from mjo.eja.eja_algebra import (random_eja,
813 ....: BilinearFormEJA)
817 Our inner product is "associative," which means the following for
818 a symmetric bilinear form::
820 sage: set_random_seed()
821 sage: J = random_eja()
822 sage: x,y,z = J.random_elements(3)
823 sage: (x*y).inner_product(z) == y.inner_product(x*z)
828 Ensure that this is the usual inner product for the algebras
831 sage: set_random_seed()
832 sage: J = HadamardEJA.random_instance()
833 sage: x,y = J.random_elements(2)
834 sage: actual = x.inner_product(y)
835 sage: expected = x.to_vector().inner_product(y.to_vector())
836 sage: actual == expected
839 Ensure that this is one-half of the trace inner-product in a
840 BilinearFormEJA that isn't just the reals (when ``n`` isn't
841 one). This is in Faraut and Koranyi, and also my "On the
844 sage: set_random_seed()
845 sage: J = BilinearFormEJA.random_instance()
846 sage: n = J.dimension()
847 sage: x = J.random_element()
848 sage: y = J.random_element()
849 sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
852 B
= self
._inner
_product
_matrix
853 return (B
*x
.to_vector()).inner_product(y
.to_vector())
856 def is_trivial(self
):
858 Return whether or not this algebra is trivial.
860 A trivial algebra contains only the zero element.
864 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
869 sage: J = ComplexHermitianEJA(3)
875 sage: J = TrivialEJA()
880 return self
.dimension() == 0
883 def multiplication_table(self
):
885 Return a visual representation of this algebra's multiplication
886 table (on basis elements).
890 sage: from mjo.eja.eja_algebra import JordanSpinEJA
894 sage: J = JordanSpinEJA(4)
895 sage: J.multiplication_table()
896 +----++----+----+----+----+
897 | * || e0 | e1 | e2 | e3 |
898 +====++====+====+====+====+
899 | e0 || e0 | e1 | e2 | e3 |
900 +----++----+----+----+----+
901 | e1 || e1 | e0 | 0 | 0 |
902 +----++----+----+----+----+
903 | e2 || e2 | 0 | e0 | 0 |
904 +----++----+----+----+----+
905 | e3 || e3 | 0 | 0 | e0 |
906 +----++----+----+----+----+
910 # Prepend the header row.
911 M
= [["*"] + list(self
.gens())]
913 # And to each subsequent row, prepend an entry that belongs to
914 # the left-side "header column."
915 M
+= [ [self
.monomial(i
)] + [ self
.monomial(i
)*self
.monomial(j
)
919 return table(M
, header_row
=True, header_column
=True, frame
=True)
922 def matrix_basis(self
):
924 Return an (often more natural) representation of this algebras
925 basis as an ordered tuple of matrices.
927 Every finite-dimensional Euclidean Jordan Algebra is a, up to
928 Jordan isomorphism, a direct sum of five simple
929 algebras---four of which comprise Hermitian matrices. And the
930 last type of algebra can of course be thought of as `n`-by-`1`
931 column matrices (ambiguusly called column vectors) to avoid
932 special cases. As a result, matrices (and column vectors) are
933 a natural representation format for Euclidean Jordan algebra
936 But, when we construct an algebra from a basis of matrices,
937 those matrix representations are lost in favor of coordinate
938 vectors *with respect to* that basis. We could eventually
939 convert back if we tried hard enough, but having the original
940 representations handy is valuable enough that we simply store
941 them and return them from this method.
943 Why implement this for non-matrix algebras? Avoiding special
944 cases for the :class:`BilinearFormEJA` pays with simplicity in
945 its own right. But mainly, we would like to be able to assume
946 that elements of a :class:`CartesianProductEJA` can be displayed
947 nicely, without having to have special classes for direct sums
948 one of whose components was a matrix algebra.
952 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
953 ....: RealSymmetricEJA)
957 sage: J = RealSymmetricEJA(2)
959 Finite family {0: e0, 1: e1, 2: e2}
960 sage: J.matrix_basis()
962 [1 0] [ 0 0.7071067811865475?] [0 0]
963 [0 0], [0.7071067811865475? 0], [0 1]
968 sage: J = JordanSpinEJA(2)
970 Finite family {0: e0, 1: e1}
971 sage: J.matrix_basis()
977 return self
._matrix
_basis
980 def matrix_space(self
):
982 Return the matrix space in which this algebra's elements live, if
983 we think of them as matrices (including column vectors of the
986 "By default" this will be an `n`-by-`1` column-matrix space,
987 except when the algebra is trivial. There it's `n`-by-`n`
988 (where `n` is zero), to ensure that two elements of the matrix
989 space (empty matrices) can be multiplied. For algebras of
990 matrices, this returns the space in which their
991 real embeddings live.
995 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
997 ....: QuaternionHermitianEJA,
1002 By default, the matrix representation is just a column-matrix
1003 equivalent to the vector representation::
1005 sage: J = JordanSpinEJA(3)
1006 sage: J.matrix_space()
1007 Full MatrixSpace of 3 by 1 dense matrices over Algebraic
1010 The matrix representation in the trivial algebra is
1011 zero-by-zero instead of the usual `n`-by-one::
1013 sage: J = TrivialEJA()
1014 sage: J.matrix_space()
1015 Full MatrixSpace of 0 by 0 dense matrices over Algebraic
1018 The matrix space for complex/quaternion Hermitian matrix EJA
1019 is the space in which their real-embeddings live, not the
1020 original complex/quaternion matrix space::
1022 sage: J = ComplexHermitianEJA(2,field=QQ,orthonormalize=False)
1023 sage: J.matrix_space()
1024 Full MatrixSpace of 4 by 4 dense matrices over Rational Field
1025 sage: J = QuaternionHermitianEJA(1,field=QQ,orthonormalize=False)
1026 sage: J.matrix_space()
1027 Full MatrixSpace of 4 by 4 dense matrices over Rational Field
1030 if self
.is_trivial():
1031 return MatrixSpace(self
.base_ring(), 0)
1033 return self
.matrix_basis()[0].parent()
1039 Return the unit element of this algebra.
1043 sage: from mjo.eja.eja_algebra import (HadamardEJA,
1048 We can compute unit element in the Hadamard EJA::
1050 sage: J = HadamardEJA(5)
1052 e0 + e1 + e2 + e3 + e4
1054 The unit element in the Hadamard EJA is inherited in the
1055 subalgebras generated by its elements::
1057 sage: J = HadamardEJA(5)
1059 e0 + e1 + e2 + e3 + e4
1060 sage: x = sum(J.gens())
1061 sage: A = x.subalgebra_generated_by(orthonormalize=False)
1064 sage: A.one().superalgebra_element()
1065 e0 + e1 + e2 + e3 + e4
1069 The identity element acts like the identity, regardless of
1070 whether or not we orthonormalize::
1072 sage: set_random_seed()
1073 sage: J = random_eja()
1074 sage: x = J.random_element()
1075 sage: J.one()*x == x and x*J.one() == x
1077 sage: A = x.subalgebra_generated_by()
1078 sage: y = A.random_element()
1079 sage: A.one()*y == y and y*A.one() == y
1084 sage: set_random_seed()
1085 sage: J = random_eja(field=QQ, orthonormalize=False)
1086 sage: x = J.random_element()
1087 sage: J.one()*x == x and x*J.one() == x
1089 sage: A = x.subalgebra_generated_by(orthonormalize=False)
1090 sage: y = A.random_element()
1091 sage: A.one()*y == y and y*A.one() == y
1094 The matrix of the unit element's operator is the identity,
1095 regardless of the base field and whether or not we
1098 sage: set_random_seed()
1099 sage: J = random_eja()
1100 sage: actual = J.one().operator().matrix()
1101 sage: expected = matrix.identity(J.base_ring(), J.dimension())
1102 sage: actual == expected
1104 sage: x = J.random_element()
1105 sage: A = x.subalgebra_generated_by()
1106 sage: actual = A.one().operator().matrix()
1107 sage: expected = matrix.identity(A.base_ring(), A.dimension())
1108 sage: actual == expected
1113 sage: set_random_seed()
1114 sage: J = random_eja(field=QQ, orthonormalize=False)
1115 sage: actual = J.one().operator().matrix()
1116 sage: expected = matrix.identity(J.base_ring(), J.dimension())
1117 sage: actual == expected
1119 sage: x = J.random_element()
1120 sage: A = x.subalgebra_generated_by(orthonormalize=False)
1121 sage: actual = A.one().operator().matrix()
1122 sage: expected = matrix.identity(A.base_ring(), A.dimension())
1123 sage: actual == expected
1126 Ensure that the cached unit element (often precomputed by
1127 hand) agrees with the computed one::
1129 sage: set_random_seed()
1130 sage: J = random_eja()
1131 sage: cached = J.one()
1132 sage: J.one.clear_cache()
1133 sage: J.one() == cached
1138 sage: set_random_seed()
1139 sage: J = random_eja(field=QQ, orthonormalize=False)
1140 sage: cached = J.one()
1141 sage: J.one.clear_cache()
1142 sage: J.one() == cached
1146 # We can brute-force compute the matrices of the operators
1147 # that correspond to the basis elements of this algebra.
1148 # If some linear combination of those basis elements is the
1149 # algebra identity, then the same linear combination of
1150 # their matrices has to be the identity matrix.
1152 # Of course, matrices aren't vectors in sage, so we have to
1153 # appeal to the "long vectors" isometry.
1154 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
1156 # Now we use basic linear algebra to find the coefficients,
1157 # of the matrices-as-vectors-linear-combination, which should
1158 # work for the original algebra basis too.
1159 A
= matrix(self
.base_ring(), oper_vecs
)
1161 # We used the isometry on the left-hand side already, but we
1162 # still need to do it for the right-hand side. Recall that we
1163 # wanted something that summed to the identity matrix.
1164 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
1166 # Now if there's an identity element in the algebra, this
1167 # should work. We solve on the left to avoid having to
1168 # transpose the matrix "A".
1169 return self
.from_vector(A
.solve_left(b
))
1172 def peirce_decomposition(self
, c
):
1174 The Peirce decomposition of this algebra relative to the
1177 In the future, this can be extended to a complete system of
1178 orthogonal idempotents.
1182 - ``c`` -- an idempotent of this algebra.
1186 A triple (J0, J5, J1) containing two subalgebras and one subspace
1189 - ``J0`` -- the algebra on the eigenspace of ``c.operator()``
1190 corresponding to the eigenvalue zero.
1192 - ``J5`` -- the eigenspace (NOT a subalgebra) of ``c.operator()``
1193 corresponding to the eigenvalue one-half.
1195 - ``J1`` -- the algebra on the eigenspace of ``c.operator()``
1196 corresponding to the eigenvalue one.
1198 These are the only possible eigenspaces for that operator, and this
1199 algebra is a direct sum of them. The spaces ``J0`` and ``J1`` are
1200 orthogonal, and are subalgebras of this algebra with the appropriate
1205 sage: from mjo.eja.eja_algebra import random_eja, RealSymmetricEJA
1209 The canonical example comes from the symmetric matrices, which
1210 decompose into diagonal and off-diagonal parts::
1212 sage: J = RealSymmetricEJA(3)
1213 sage: C = matrix(QQ, [ [1,0,0],
1217 sage: J0,J5,J1 = J.peirce_decomposition(c)
1219 Euclidean Jordan algebra of dimension 1...
1221 Vector space of degree 6 and dimension 2...
1223 Euclidean Jordan algebra of dimension 3...
1224 sage: J0.one().to_matrix()
1228 sage: orig_df = AA.options.display_format
1229 sage: AA.options.display_format = 'radical'
1230 sage: J.from_vector(J5.basis()[0]).to_matrix()
1234 sage: J.from_vector(J5.basis()[1]).to_matrix()
1238 sage: AA.options.display_format = orig_df
1239 sage: J1.one().to_matrix()
1246 Every algebra decomposes trivially with respect to its identity
1249 sage: set_random_seed()
1250 sage: J = random_eja()
1251 sage: J0,J5,J1 = J.peirce_decomposition(J.one())
1252 sage: J0.dimension() == 0 and J5.dimension() == 0
1254 sage: J1.superalgebra() == J and J1.dimension() == J.dimension()
1257 The decomposition is into eigenspaces, and its components are
1258 therefore necessarily orthogonal. Moreover, the identity
1259 elements in the two subalgebras are the projections onto their
1260 respective subspaces of the superalgebra's identity element::
1262 sage: set_random_seed()
1263 sage: J = random_eja()
1264 sage: x = J.random_element()
1265 sage: if not J.is_trivial():
1266 ....: while x.is_nilpotent():
1267 ....: x = J.random_element()
1268 sage: c = x.subalgebra_idempotent()
1269 sage: J0,J5,J1 = J.peirce_decomposition(c)
1271 sage: for (w,y,z) in zip(J0.basis(), J5.basis(), J1.basis()):
1272 ....: w = w.superalgebra_element()
1273 ....: y = J.from_vector(y)
1274 ....: z = z.superalgebra_element()
1275 ....: ipsum += w.inner_product(y).abs()
1276 ....: ipsum += w.inner_product(z).abs()
1277 ....: ipsum += y.inner_product(z).abs()
1280 sage: J1(c) == J1.one()
1282 sage: J0(J.one() - c) == J0.one()
1286 if not c
.is_idempotent():
1287 raise ValueError("element is not idempotent: %s" % c
)
1289 # Default these to what they should be if they turn out to be
1290 # trivial, because eigenspaces_left() won't return eigenvalues
1291 # corresponding to trivial spaces (e.g. it returns only the
1292 # eigenspace corresponding to lambda=1 if you take the
1293 # decomposition relative to the identity element).
1294 trivial
= self
.subalgebra(())
1295 J0
= trivial
# eigenvalue zero
1296 J5
= VectorSpace(self
.base_ring(), 0) # eigenvalue one-half
1297 J1
= trivial
# eigenvalue one
1299 for (eigval
, eigspace
) in c
.operator().matrix().right_eigenspaces():
1300 if eigval
== ~
(self
.base_ring()(2)):
1303 gens
= tuple( self
.from_vector(b
) for b
in eigspace
.basis() )
1304 subalg
= self
.subalgebra(gens
, check_axioms
=False)
1310 raise ValueError("unexpected eigenvalue: %s" % eigval
)
1315 def random_element(self
, thorough
=False):
1317 Return a random element of this algebra.
1319 Our algebra superclass method only returns a linear
1320 combination of at most two basis elements. We instead
1321 want the vector space "random element" method that
1322 returns a more diverse selection.
1326 - ``thorough`` -- (boolean; default False) whether or not we
1327 should generate irrational coefficients for the random
1328 element when our base ring is irrational; this slows the
1329 algebra operations to a crawl, but any truly random method
1333 # For a general base ring... maybe we can trust this to do the
1334 # right thing? Unlikely, but.
1335 V
= self
.vector_space()
1336 v
= V
.random_element()
1338 if self
.base_ring() is AA
:
1339 # The "random element" method of the algebraic reals is
1340 # stupid at the moment, and only returns integers between
1341 # -2 and 2, inclusive:
1343 # https://trac.sagemath.org/ticket/30875
1345 # Instead, we implement our own "random vector" method,
1346 # and then coerce that into the algebra. We use the vector
1347 # space degree here instead of the dimension because a
1348 # subalgebra could (for example) be spanned by only two
1349 # vectors, each with five coordinates. We need to
1350 # generate all five coordinates.
1352 v
*= QQbar
.random_element().real()
1354 v
*= QQ
.random_element()
1356 return self
.from_vector(V
.coordinate_vector(v
))
1358 def random_elements(self
, count
, thorough
=False):
1360 Return ``count`` random elements as a tuple.
1364 - ``thorough`` -- (boolean; default False) whether or not we
1365 should generate irrational coefficients for the random
1366 elements when our base ring is irrational; this slows the
1367 algebra operations to a crawl, but any truly random method
1372 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1376 sage: J = JordanSpinEJA(3)
1377 sage: x,y,z = J.random_elements(3)
1378 sage: all( [ x in J, y in J, z in J ])
1380 sage: len( J.random_elements(10) ) == 10
1384 return tuple( self
.random_element(thorough
)
1385 for idx
in range(count
) )
1389 def _charpoly_coefficients(self
):
1391 The `r` polynomial coefficients of the "characteristic polynomial
1396 sage: from mjo.eja.eja_algebra import random_eja
1400 The theory shows that these are all homogeneous polynomials of
1403 sage: set_random_seed()
1404 sage: J = random_eja()
1405 sage: all(p.is_homogeneous() for p in J._charpoly_coefficients())
1409 n
= self
.dimension()
1410 R
= self
.coordinate_polynomial_ring()
1412 F
= R
.fraction_field()
1415 # From a result in my book, these are the entries of the
1416 # basis representation of L_x.
1417 return sum( vars[k
]*self
.monomial(k
).operator().matrix()[i
,j
]
1420 L_x
= matrix(F
, n
, n
, L_x_i_j
)
1423 if self
.rank
.is_in_cache():
1425 # There's no need to pad the system with redundant
1426 # columns if we *know* they'll be redundant.
1429 # Compute an extra power in case the rank is equal to
1430 # the dimension (otherwise, we would stop at x^(r-1)).
1431 x_powers
= [ (L_x
**k
)*self
.one().to_vector()
1432 for k
in range(n
+1) ]
1433 A
= matrix
.column(F
, x_powers
[:n
])
1434 AE
= A
.extended_echelon_form()
1441 # The theory says that only the first "r" coefficients are
1442 # nonzero, and they actually live in the original polynomial
1443 # ring and not the fraction field. We negate them because in
1444 # the actual characteristic polynomial, they get moved to the
1445 # other side where x^r lives. We don't bother to trim A_rref
1446 # down to a square matrix and solve the resulting system,
1447 # because the upper-left r-by-r portion of A_rref is
1448 # guaranteed to be the identity matrix, so e.g.
1450 # A_rref.solve_right(Y)
1452 # would just be returning Y.
1453 return (-E
*b
)[:r
].change_ring(R
)
1458 Return the rank of this EJA.
1460 This is a cached method because we know the rank a priori for
1461 all of the algebras we can construct. Thus we can avoid the
1462 expensive ``_charpoly_coefficients()`` call unless we truly
1463 need to compute the whole characteristic polynomial.
1467 sage: from mjo.eja.eja_algebra import (HadamardEJA,
1468 ....: JordanSpinEJA,
1469 ....: RealSymmetricEJA,
1470 ....: ComplexHermitianEJA,
1471 ....: QuaternionHermitianEJA,
1476 The rank of the Jordan spin algebra is always two::
1478 sage: JordanSpinEJA(2).rank()
1480 sage: JordanSpinEJA(3).rank()
1482 sage: JordanSpinEJA(4).rank()
1485 The rank of the `n`-by-`n` Hermitian real, complex, or
1486 quaternion matrices is `n`::
1488 sage: RealSymmetricEJA(4).rank()
1490 sage: ComplexHermitianEJA(3).rank()
1492 sage: QuaternionHermitianEJA(2).rank()
1497 Ensure that every EJA that we know how to construct has a
1498 positive integer rank, unless the algebra is trivial in
1499 which case its rank will be zero::
1501 sage: set_random_seed()
1502 sage: J = random_eja()
1506 sage: r > 0 or (r == 0 and J.is_trivial())
1509 Ensure that computing the rank actually works, since the ranks
1510 of all simple algebras are known and will be cached by default::
1512 sage: set_random_seed() # long time
1513 sage: J = random_eja() # long time
1514 sage: cached = J.rank() # long time
1515 sage: J.rank.clear_cache() # long time
1516 sage: J.rank() == cached # long time
1520 return len(self
._charpoly
_coefficients
())
1523 def subalgebra(self
, basis
, **kwargs
):
1525 Create a subalgebra of this algebra from the given basis.
1527 from mjo
.eja
.eja_subalgebra
import FiniteDimensionalEJASubalgebra
1528 return FiniteDimensionalEJASubalgebra(self
, basis
, **kwargs
)
1531 def vector_space(self
):
1533 Return the vector space that underlies this algebra.
1537 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1541 sage: J = RealSymmetricEJA(2)
1542 sage: J.vector_space()
1543 Vector space of dimension 3 over...
1546 return self
.zero().to_vector().parent().ambient_vector_space()
1550 class RationalBasisEJA(FiniteDimensionalEJA
):
1552 New class for algebras whose supplied basis elements have all rational entries.
1556 sage: from mjo.eja.eja_algebra import BilinearFormEJA
1560 The supplied basis is orthonormalized by default::
1562 sage: B = matrix(QQ, [[1, 0, 0], [0, 25, -32], [0, -32, 41]])
1563 sage: J = BilinearFormEJA(B)
1564 sage: J.matrix_basis()
1581 # Abuse the check_field parameter to check that the entries of
1582 # out basis (in ambient coordinates) are in the field QQ.
1583 if not all( all(b_i
in QQ
for b_i
in b
.list()) for b
in basis
):
1584 raise TypeError("basis not rational")
1586 super().__init
__(basis
,
1590 check_field
=check_field
,
1593 self
._rational
_algebra
= None
1595 # There's no point in constructing the extra algebra if this
1596 # one is already rational.
1598 # Note: the same Jordan and inner-products work here,
1599 # because they are necessarily defined with respect to
1600 # ambient coordinates and not any particular basis.
1601 self
._rational
_algebra
= FiniteDimensionalEJA(
1606 associative
=self
.is_associative(),
1607 orthonormalize
=False,
1612 def _charpoly_coefficients(self
):
1616 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
1617 ....: JordanSpinEJA)
1621 The base ring of the resulting polynomial coefficients is what
1622 it should be, and not the rationals (unless the algebra was
1623 already over the rationals)::
1625 sage: J = JordanSpinEJA(3)
1626 sage: J._charpoly_coefficients()
1627 (X1^2 - X2^2 - X3^2, -2*X1)
1628 sage: a0 = J._charpoly_coefficients()[0]
1630 Algebraic Real Field
1631 sage: a0.base_ring()
1632 Algebraic Real Field
1635 if self
._rational
_algebra
is None:
1636 # There's no need to construct *another* algebra over the
1637 # rationals if this one is already over the
1638 # rationals. Likewise, if we never orthonormalized our
1639 # basis, we might as well just use the given one.
1640 return super()._charpoly
_coefficients
()
1642 # Do the computation over the rationals. The answer will be
1643 # the same, because all we've done is a change of basis.
1644 # Then, change back from QQ to our real base ring
1645 a
= ( a_i
.change_ring(self
.base_ring())
1646 for a_i
in self
._rational
_algebra
._charpoly
_coefficients
() )
1648 if self
._deortho
_matrix
is None:
1649 # This can happen if our base ring was, say, AA and we
1650 # chose not to (or didn't need to) orthonormalize. It's
1651 # still faster to do the computations over QQ even if
1652 # the numbers in the boxes stay the same.
1655 # Otherwise, convert the coordinate variables back to the
1656 # deorthonormalized ones.
1657 R
= self
.coordinate_polynomial_ring()
1658 from sage
.modules
.free_module_element
import vector
1659 X
= vector(R
, R
.gens())
1660 BX
= self
._deortho
_matrix
*X
1662 subs_dict
= { X[i]: BX[i] for i in range(len(X)) }
1663 return tuple( a_i
.subs(subs_dict
) for a_i
in a
)
1665 class ConcreteEJA(RationalBasisEJA
):
1667 A class for the Euclidean Jordan algebras that we know by name.
1669 These are the Jordan algebras whose basis, multiplication table,
1670 rank, and so on are known a priori. More to the point, they are
1671 the Euclidean Jordan algebras for which we are able to conjure up
1672 a "random instance."
1676 sage: from mjo.eja.eja_algebra import ConcreteEJA
1680 Our basis is normalized with respect to the algebra's inner
1681 product, unless we specify otherwise::
1683 sage: set_random_seed()
1684 sage: J = ConcreteEJA.random_instance()
1685 sage: all( b.norm() == 1 for b in J.gens() )
1688 Since our basis is orthonormal with respect to the algebra's inner
1689 product, and since we know that this algebra is an EJA, any
1690 left-multiplication operator's matrix will be symmetric because
1691 natural->EJA basis representation is an isometry and within the
1692 EJA the operator is self-adjoint by the Jordan axiom::
1694 sage: set_random_seed()
1695 sage: J = ConcreteEJA.random_instance()
1696 sage: x = J.random_element()
1697 sage: x.operator().is_self_adjoint()
1702 def _max_random_instance_size():
1704 Return an integer "size" that is an upper bound on the size of
1705 this algebra when it is used in a random test
1706 case. Unfortunately, the term "size" is ambiguous -- when
1707 dealing with `R^n` under either the Hadamard or Jordan spin
1708 product, the "size" refers to the dimension `n`. When dealing
1709 with a matrix algebra (real symmetric or complex/quaternion
1710 Hermitian), it refers to the size of the matrix, which is far
1711 less than the dimension of the underlying vector space.
1713 This method must be implemented in each subclass.
1715 raise NotImplementedError
1718 def random_instance(cls
, *args
, **kwargs
):
1720 Return a random instance of this type of algebra.
1722 This method should be implemented in each subclass.
1724 from sage
.misc
.prandom
import choice
1725 eja_class
= choice(cls
.__subclasses
__())
1727 # These all bubble up to the RationalBasisEJA superclass
1728 # constructor, so any (kw)args valid there are also valid
1730 return eja_class
.random_instance(*args
, **kwargs
)
1735 def dimension_over_reals():
1737 The dimension of this matrix's base ring over the reals.
1739 The reals are dimension one over themselves, obviously; that's
1740 just `\mathbb{R}^{1}`. Likewise, the complex numbers `a + bi`
1741 have dimension two. Finally, the quaternions have dimension
1742 four over the reals.
1744 This is used to determine the size of the matrix returned from
1745 :meth:`real_embed`, among other things.
1747 raise NotImplementedError
1750 def real_embed(cls
,M
):
1752 Embed the matrix ``M`` into a space of real matrices.
1754 The matrix ``M`` can have entries in any field at the moment:
1755 the real numbers, complex numbers, or quaternions. And although
1756 they are not a field, we can probably support octonions at some
1757 point, too. This function returns a real matrix that "acts like"
1758 the original with respect to matrix multiplication; i.e.
1760 real_embed(M*N) = real_embed(M)*real_embed(N)
1763 if M
.ncols() != M
.nrows():
1764 raise ValueError("the matrix 'M' must be square")
1769 def real_unembed(cls
,M
):
1771 The inverse of :meth:`real_embed`.
1773 if M
.ncols() != M
.nrows():
1774 raise ValueError("the matrix 'M' must be square")
1775 if not ZZ(M
.nrows()).mod(cls
.dimension_over_reals()).is_zero():
1776 raise ValueError("the matrix 'M' must be a real embedding")
1780 def jordan_product(X
,Y
):
1781 return (X
*Y
+ Y
*X
)/2
1784 def trace_inner_product(cls
,X
,Y
):
1786 Compute the trace inner-product of two real-embeddings.
1790 sage: from mjo.eja.eja_algebra import (RealSymmetricEJA,
1791 ....: ComplexHermitianEJA,
1792 ....: QuaternionHermitianEJA)
1796 This gives the same answer as it would if we computed the trace
1797 from the unembedded (original) matrices::
1799 sage: set_random_seed()
1800 sage: J = RealSymmetricEJA.random_instance()
1801 sage: x,y = J.random_elements(2)
1802 sage: Xe = x.to_matrix()
1803 sage: Ye = y.to_matrix()
1804 sage: X = J.real_unembed(Xe)
1805 sage: Y = J.real_unembed(Ye)
1806 sage: expected = (X*Y).trace()
1807 sage: actual = J.trace_inner_product(Xe,Ye)
1808 sage: actual == expected
1813 sage: set_random_seed()
1814 sage: J = ComplexHermitianEJA.random_instance()
1815 sage: x,y = J.random_elements(2)
1816 sage: Xe = x.to_matrix()
1817 sage: Ye = y.to_matrix()
1818 sage: X = J.real_unembed(Xe)
1819 sage: Y = J.real_unembed(Ye)
1820 sage: expected = (X*Y).trace().real()
1821 sage: actual = J.trace_inner_product(Xe,Ye)
1822 sage: actual == expected
1827 sage: set_random_seed()
1828 sage: J = QuaternionHermitianEJA.random_instance()
1829 sage: x,y = J.random_elements(2)
1830 sage: Xe = x.to_matrix()
1831 sage: Ye = y.to_matrix()
1832 sage: X = J.real_unembed(Xe)
1833 sage: Y = J.real_unembed(Ye)
1834 sage: expected = (X*Y).trace().coefficient_tuple()[0]
1835 sage: actual = J.trace_inner_product(Xe,Ye)
1836 sage: actual == expected
1840 Xu
= cls
.real_unembed(X
)
1841 Yu
= cls
.real_unembed(Y
)
1842 tr
= (Xu
*Yu
).trace()
1845 # Works in QQ, AA, RDF, et cetera.
1847 except AttributeError:
1848 # A quaternion doesn't have a real() method, but does
1849 # have coefficient_tuple() method that returns the
1850 # coefficients of 1, i, j, and k -- in that order.
1851 return tr
.coefficient_tuple()[0]
1854 class RealMatrixEJA(MatrixEJA
):
1856 def dimension_over_reals():
1860 class RealSymmetricEJA(ConcreteEJA
, RealMatrixEJA
):
1862 The rank-n simple EJA consisting of real symmetric n-by-n
1863 matrices, the usual symmetric Jordan product, and the trace inner
1864 product. It has dimension `(n^2 + n)/2` over the reals.
1868 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1872 sage: J = RealSymmetricEJA(2)
1873 sage: e0, e1, e2 = J.gens()
1881 In theory, our "field" can be any subfield of the reals::
1883 sage: RealSymmetricEJA(2, field=RDF, check_axioms=True)
1884 Euclidean Jordan algebra of dimension 3 over Real Double Field
1885 sage: RealSymmetricEJA(2, field=RR, check_axioms=True)
1886 Euclidean Jordan algebra of dimension 3 over Real Field with
1887 53 bits of precision
1891 The dimension of this algebra is `(n^2 + n) / 2`::
1893 sage: set_random_seed()
1894 sage: n_max = RealSymmetricEJA._max_random_instance_size()
1895 sage: n = ZZ.random_element(1, n_max)
1896 sage: J = RealSymmetricEJA(n)
1897 sage: J.dimension() == (n^2 + n)/2
1900 The Jordan multiplication is what we think it is::
1902 sage: set_random_seed()
1903 sage: J = RealSymmetricEJA.random_instance()
1904 sage: x,y = J.random_elements(2)
1905 sage: actual = (x*y).to_matrix()
1906 sage: X = x.to_matrix()
1907 sage: Y = y.to_matrix()
1908 sage: expected = (X*Y + Y*X)/2
1909 sage: actual == expected
1911 sage: J(expected) == x*y
1914 We can change the generator prefix::
1916 sage: RealSymmetricEJA(3, prefix='q').gens()
1917 (q0, q1, q2, q3, q4, q5)
1919 We can construct the (trivial) algebra of rank zero::
1921 sage: RealSymmetricEJA(0)
1922 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1926 def _denormalized_basis(cls
, n
):
1928 Return a basis for the space of real symmetric n-by-n matrices.
1932 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1936 sage: set_random_seed()
1937 sage: n = ZZ.random_element(1,5)
1938 sage: B = RealSymmetricEJA._denormalized_basis(n)
1939 sage: all( M.is_symmetric() for M in B)
1943 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1947 for j
in range(i
+1):
1948 Eij
= matrix(ZZ
, n
, lambda k
,l
: k
==i
and l
==j
)
1952 Sij
= Eij
+ Eij
.transpose()
1958 def _max_random_instance_size():
1959 return 4 # Dimension 10
1962 def random_instance(cls
, **kwargs
):
1964 Return a random instance of this type of algebra.
1966 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1967 return cls(n
, **kwargs
)
1969 def __init__(self
, n
, **kwargs
):
1970 # We know this is a valid EJA, but will double-check
1971 # if the user passes check_axioms=True.
1972 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
1978 super().__init
__(self
._denormalized
_basis
(n
),
1979 self
.jordan_product
,
1980 self
.trace_inner_product
,
1981 associative
=associative
,
1984 # TODO: this could be factored out somehow, but is left here
1985 # because the MatrixEJA is not presently a subclass of the
1986 # FDEJA class that defines rank() and one().
1987 self
.rank
.set_cache(n
)
1988 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
1989 self
.one
.set_cache(self(idV
))
1993 class ComplexMatrixEJA(MatrixEJA
):
1994 # A manual dictionary-cache for the complex_extension() method,
1995 # since apparently @classmethods can't also be @cached_methods.
1996 _complex_extension
= {}
1999 def complex_extension(cls
,field
):
2001 The complex field that we embed/unembed, as an extension
2002 of the given ``field``.
2004 if field
in cls
._complex
_extension
:
2005 return cls
._complex
_extension
[field
]
2007 # Sage doesn't know how to adjoin the complex "i" (the root of
2008 # x^2 + 1) to a field in a general way. Here, we just enumerate
2009 # all of the cases that I have cared to support so far.
2011 # Sage doesn't know how to embed AA into QQbar, i.e. how
2012 # to adjoin sqrt(-1) to AA.
2014 elif not field
.is_exact():
2016 F
= field
.complex_field()
2018 # Works for QQ and... maybe some other fields.
2019 R
= PolynomialRing(field
, 'z')
2021 F
= field
.extension(z
**2 + 1, 'I', embedding
=CLF(-1).sqrt())
2023 cls
._complex
_extension
[field
] = F
2027 def dimension_over_reals():
2031 def real_embed(cls
,M
):
2033 Embed the n-by-n complex matrix ``M`` into the space of real
2034 matrices of size 2n-by-2n via the map the sends each entry `z = a +
2035 bi` to the block matrix ``[[a,b],[-b,a]]``.
2039 sage: from mjo.eja.eja_algebra import ComplexMatrixEJA
2043 sage: F = QuadraticField(-1, 'I')
2044 sage: x1 = F(4 - 2*i)
2045 sage: x2 = F(1 + 2*i)
2048 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
2049 sage: ComplexMatrixEJA.real_embed(M)
2058 Embedding is a homomorphism (isomorphism, in fact)::
2060 sage: set_random_seed()
2061 sage: n = ZZ.random_element(3)
2062 sage: F = QuadraticField(-1, 'I')
2063 sage: X = random_matrix(F, n)
2064 sage: Y = random_matrix(F, n)
2065 sage: Xe = ComplexMatrixEJA.real_embed(X)
2066 sage: Ye = ComplexMatrixEJA.real_embed(Y)
2067 sage: XYe = ComplexMatrixEJA.real_embed(X*Y)
2072 super().real_embed(M
)
2075 # We don't need any adjoined elements...
2076 field
= M
.base_ring().base_ring()
2082 blocks
.append(matrix(field
, 2, [ [ a
, b
],
2085 return matrix
.block(field
, n
, blocks
)
2089 def real_unembed(cls
,M
):
2091 The inverse of _embed_complex_matrix().
2095 sage: from mjo.eja.eja_algebra import ComplexMatrixEJA
2099 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
2100 ....: [-2, 1, -4, 3],
2101 ....: [ 9, 10, 11, 12],
2102 ....: [-10, 9, -12, 11] ])
2103 sage: ComplexMatrixEJA.real_unembed(A)
2105 [ 10*I + 9 12*I + 11]
2109 Unembedding is the inverse of embedding::
2111 sage: set_random_seed()
2112 sage: F = QuadraticField(-1, 'I')
2113 sage: M = random_matrix(F, 3)
2114 sage: Me = ComplexMatrixEJA.real_embed(M)
2115 sage: ComplexMatrixEJA.real_unembed(Me) == M
2119 super().real_unembed(M
)
2121 d
= cls
.dimension_over_reals()
2122 F
= cls
.complex_extension(M
.base_ring())
2125 # Go top-left to bottom-right (reading order), converting every
2126 # 2-by-2 block we see to a single complex element.
2128 for k
in range(n
/d
):
2129 for j
in range(n
/d
):
2130 submat
= M
[d
*k
:d
*k
+d
,d
*j
:d
*j
+d
]
2131 if submat
[0,0] != submat
[1,1]:
2132 raise ValueError('bad on-diagonal submatrix')
2133 if submat
[0,1] != -submat
[1,0]:
2134 raise ValueError('bad off-diagonal submatrix')
2135 z
= submat
[0,0] + submat
[0,1]*i
2138 return matrix(F
, n
/d
, elements
)
2141 class ComplexHermitianEJA(ConcreteEJA
, ComplexMatrixEJA
):
2143 The rank-n simple EJA consisting of complex Hermitian n-by-n
2144 matrices over the real numbers, the usual symmetric Jordan product,
2145 and the real-part-of-trace inner product. It has dimension `n^2` over
2150 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
2154 In theory, our "field" can be any subfield of the reals::
2156 sage: ComplexHermitianEJA(2, field=RDF, check_axioms=True)
2157 Euclidean Jordan algebra of dimension 4 over Real Double Field
2158 sage: ComplexHermitianEJA(2, field=RR, check_axioms=True)
2159 Euclidean Jordan algebra of dimension 4 over Real Field with
2160 53 bits of precision
2164 The dimension of this algebra is `n^2`::
2166 sage: set_random_seed()
2167 sage: n_max = ComplexHermitianEJA._max_random_instance_size()
2168 sage: n = ZZ.random_element(1, n_max)
2169 sage: J = ComplexHermitianEJA(n)
2170 sage: J.dimension() == n^2
2173 The Jordan multiplication is what we think it is::
2175 sage: set_random_seed()
2176 sage: J = ComplexHermitianEJA.random_instance()
2177 sage: x,y = J.random_elements(2)
2178 sage: actual = (x*y).to_matrix()
2179 sage: X = x.to_matrix()
2180 sage: Y = y.to_matrix()
2181 sage: expected = (X*Y + Y*X)/2
2182 sage: actual == expected
2184 sage: J(expected) == x*y
2187 We can change the generator prefix::
2189 sage: ComplexHermitianEJA(2, prefix='z').gens()
2192 We can construct the (trivial) algebra of rank zero::
2194 sage: ComplexHermitianEJA(0)
2195 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2200 def _denormalized_basis(cls
, n
):
2202 Returns a basis for the space of complex Hermitian n-by-n matrices.
2204 Why do we embed these? Basically, because all of numerical linear
2205 algebra assumes that you're working with vectors consisting of `n`
2206 entries from a field and scalars from the same field. There's no way
2207 to tell SageMath that (for example) the vectors contain complex
2208 numbers, while the scalar field is real.
2212 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
2216 sage: set_random_seed()
2217 sage: n = ZZ.random_element(1,5)
2218 sage: B = ComplexHermitianEJA._denormalized_basis(n)
2219 sage: all( M.is_symmetric() for M in B)
2224 R
= PolynomialRing(field
, 'z')
2226 F
= field
.extension(z
**2 + 1, 'I')
2229 # This is like the symmetric case, but we need to be careful:
2231 # * We want conjugate-symmetry, not just symmetry.
2232 # * The diagonal will (as a result) be real.
2235 Eij
= matrix
.zero(F
,n
)
2237 for j
in range(i
+1):
2241 Sij
= cls
.real_embed(Eij
)
2244 # The second one has a minus because it's conjugated.
2245 Eij
[j
,i
] = 1 # Eij = Eij + Eij.transpose()
2246 Sij_real
= cls
.real_embed(Eij
)
2248 # Eij = I*Eij - I*Eij.transpose()
2251 Sij_imag
= cls
.real_embed(Eij
)
2257 # Since we embedded these, we can drop back to the "field" that we
2258 # started with instead of the complex extension "F".
2259 return tuple( s
.change_ring(field
) for s
in S
)
2262 def __init__(self
, n
, **kwargs
):
2263 # We know this is a valid EJA, but will double-check
2264 # if the user passes check_axioms=True.
2265 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2271 super().__init
__(self
._denormalized
_basis
(n
),
2272 self
.jordan_product
,
2273 self
.trace_inner_product
,
2274 associative
=associative
,
2276 # TODO: this could be factored out somehow, but is left here
2277 # because the MatrixEJA is not presently a subclass of the
2278 # FDEJA class that defines rank() and one().
2279 self
.rank
.set_cache(n
)
2280 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
2281 self
.one
.set_cache(self(idV
))
2284 def _max_random_instance_size():
2285 return 3 # Dimension 9
2288 def random_instance(cls
, **kwargs
):
2290 Return a random instance of this type of algebra.
2292 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2293 return cls(n
, **kwargs
)
2295 class QuaternionMatrixEJA(MatrixEJA
):
2297 # A manual dictionary-cache for the quaternion_extension() method,
2298 # since apparently @classmethods can't also be @cached_methods.
2299 _quaternion_extension
= {}
2302 def quaternion_extension(cls
,field
):
2304 The quaternion field that we embed/unembed, as an extension
2305 of the given ``field``.
2307 if field
in cls
._quaternion
_extension
:
2308 return cls
._quaternion
_extension
[field
]
2310 Q
= QuaternionAlgebra(field
,-1,-1)
2312 cls
._quaternion
_extension
[field
] = Q
2316 def dimension_over_reals():
2320 def real_embed(cls
,M
):
2322 Embed the n-by-n quaternion matrix ``M`` into the space of real
2323 matrices of size 4n-by-4n by first sending each quaternion entry `z
2324 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
2325 c+di],[-c + di, a-bi]]`, and then embedding those into a real
2330 sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA
2334 sage: Q = QuaternionAlgebra(QQ,-1,-1)
2335 sage: i,j,k = Q.gens()
2336 sage: x = 1 + 2*i + 3*j + 4*k
2337 sage: M = matrix(Q, 1, [[x]])
2338 sage: QuaternionMatrixEJA.real_embed(M)
2344 Embedding is a homomorphism (isomorphism, in fact)::
2346 sage: set_random_seed()
2347 sage: n = ZZ.random_element(2)
2348 sage: Q = QuaternionAlgebra(QQ,-1,-1)
2349 sage: X = random_matrix(Q, n)
2350 sage: Y = random_matrix(Q, n)
2351 sage: Xe = QuaternionMatrixEJA.real_embed(X)
2352 sage: Ye = QuaternionMatrixEJA.real_embed(Y)
2353 sage: XYe = QuaternionMatrixEJA.real_embed(X*Y)
2358 super().real_embed(M
)
2359 quaternions
= M
.base_ring()
2362 F
= QuadraticField(-1, 'I')
2367 t
= z
.coefficient_tuple()
2372 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
2373 [-c
+ d
*i
, a
- b
*i
]])
2374 realM
= ComplexMatrixEJA
.real_embed(cplxM
)
2375 blocks
.append(realM
)
2377 # We should have real entries by now, so use the realest field
2378 # we've got for the return value.
2379 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
2384 def real_unembed(cls
,M
):
2386 The inverse of _embed_quaternion_matrix().
2390 sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA
2394 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
2395 ....: [-2, 1, -4, 3],
2396 ....: [-3, 4, 1, -2],
2397 ....: [-4, -3, 2, 1]])
2398 sage: QuaternionMatrixEJA.real_unembed(M)
2399 [1 + 2*i + 3*j + 4*k]
2403 Unembedding is the inverse of embedding::
2405 sage: set_random_seed()
2406 sage: Q = QuaternionAlgebra(QQ, -1, -1)
2407 sage: M = random_matrix(Q, 3)
2408 sage: Me = QuaternionMatrixEJA.real_embed(M)
2409 sage: QuaternionMatrixEJA.real_unembed(Me) == M
2413 super().real_unembed(M
)
2415 d
= cls
.dimension_over_reals()
2417 # Use the base ring of the matrix to ensure that its entries can be
2418 # multiplied by elements of the quaternion algebra.
2419 Q
= cls
.quaternion_extension(M
.base_ring())
2422 # Go top-left to bottom-right (reading order), converting every
2423 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
2426 for l
in range(n
/d
):
2427 for m
in range(n
/d
):
2428 submat
= ComplexMatrixEJA
.real_unembed(
2429 M
[d
*l
:d
*l
+d
,d
*m
:d
*m
+d
] )
2430 if submat
[0,0] != submat
[1,1].conjugate():
2431 raise ValueError('bad on-diagonal submatrix')
2432 if submat
[0,1] != -submat
[1,0].conjugate():
2433 raise ValueError('bad off-diagonal submatrix')
2434 z
= submat
[0,0].real()
2435 z
+= submat
[0,0].imag()*i
2436 z
+= submat
[0,1].real()*j
2437 z
+= submat
[0,1].imag()*k
2440 return matrix(Q
, n
/d
, elements
)
2443 class QuaternionHermitianEJA(ConcreteEJA
, QuaternionMatrixEJA
):
2445 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
2446 matrices, the usual symmetric Jordan product, and the
2447 real-part-of-trace inner product. It has dimension `2n^2 - n` over
2452 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2456 In theory, our "field" can be any subfield of the reals::
2458 sage: QuaternionHermitianEJA(2, field=RDF, check_axioms=True)
2459 Euclidean Jordan algebra of dimension 6 over Real Double Field
2460 sage: QuaternionHermitianEJA(2, field=RR, check_axioms=True)
2461 Euclidean Jordan algebra of dimension 6 over Real Field with
2462 53 bits of precision
2466 The dimension of this algebra is `2*n^2 - n`::
2468 sage: set_random_seed()
2469 sage: n_max = QuaternionHermitianEJA._max_random_instance_size()
2470 sage: n = ZZ.random_element(1, n_max)
2471 sage: J = QuaternionHermitianEJA(n)
2472 sage: J.dimension() == 2*(n^2) - n
2475 The Jordan multiplication is what we think it is::
2477 sage: set_random_seed()
2478 sage: J = QuaternionHermitianEJA.random_instance()
2479 sage: x,y = J.random_elements(2)
2480 sage: actual = (x*y).to_matrix()
2481 sage: X = x.to_matrix()
2482 sage: Y = y.to_matrix()
2483 sage: expected = (X*Y + Y*X)/2
2484 sage: actual == expected
2486 sage: J(expected) == x*y
2489 We can change the generator prefix::
2491 sage: QuaternionHermitianEJA(2, prefix='a').gens()
2492 (a0, a1, a2, a3, a4, a5)
2494 We can construct the (trivial) algebra of rank zero::
2496 sage: QuaternionHermitianEJA(0)
2497 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2501 def _denormalized_basis(cls
, n
):
2503 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
2505 Why do we embed these? Basically, because all of numerical
2506 linear algebra assumes that you're working with vectors consisting
2507 of `n` entries from a field and scalars from the same field. There's
2508 no way to tell SageMath that (for example) the vectors contain
2509 complex numbers, while the scalar field is real.
2513 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2517 sage: set_random_seed()
2518 sage: n = ZZ.random_element(1,5)
2519 sage: B = QuaternionHermitianEJA._denormalized_basis(n)
2520 sage: all( M.is_symmetric() for M in B )
2525 Q
= QuaternionAlgebra(QQ
,-1,-1)
2528 # This is like the symmetric case, but we need to be careful:
2530 # * We want conjugate-symmetry, not just symmetry.
2531 # * The diagonal will (as a result) be real.
2534 Eij
= matrix
.zero(Q
,n
)
2536 for j
in range(i
+1):
2540 Sij
= cls
.real_embed(Eij
)
2543 # The second, third, and fourth ones have a minus
2544 # because they're conjugated.
2545 # Eij = Eij + Eij.transpose()
2547 Sij_real
= cls
.real_embed(Eij
)
2549 # Eij = I*(Eij - Eij.transpose())
2552 Sij_I
= cls
.real_embed(Eij
)
2554 # Eij = J*(Eij - Eij.transpose())
2557 Sij_J
= cls
.real_embed(Eij
)
2559 # Eij = K*(Eij - Eij.transpose())
2562 Sij_K
= cls
.real_embed(Eij
)
2568 # Since we embedded these, we can drop back to the "field" that we
2569 # started with instead of the quaternion algebra "Q".
2570 return tuple( s
.change_ring(field
) for s
in S
)
2573 def __init__(self
, n
, **kwargs
):
2574 # We know this is a valid EJA, but will double-check
2575 # if the user passes check_axioms=True.
2576 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2582 super().__init
__(self
._denormalized
_basis
(n
),
2583 self
.jordan_product
,
2584 self
.trace_inner_product
,
2585 associative
=associative
,
2588 # TODO: this could be factored out somehow, but is left here
2589 # because the MatrixEJA is not presently a subclass of the
2590 # FDEJA class that defines rank() and one().
2591 self
.rank
.set_cache(n
)
2592 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
2593 self
.one
.set_cache(self(idV
))
2597 def _max_random_instance_size():
2599 The maximum rank of a random QuaternionHermitianEJA.
2601 return 2 # Dimension 6
2604 def random_instance(cls
, **kwargs
):
2606 Return a random instance of this type of algebra.
2608 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2609 return cls(n
, **kwargs
)
2612 class HadamardEJA(ConcreteEJA
):
2614 Return the Euclidean Jordan Algebra corresponding to the set
2615 `R^n` under the Hadamard product.
2617 Note: this is nothing more than the Cartesian product of ``n``
2618 copies of the spin algebra. Once Cartesian product algebras
2619 are implemented, this can go.
2623 sage: from mjo.eja.eja_algebra import HadamardEJA
2627 This multiplication table can be verified by hand::
2629 sage: J = HadamardEJA(3)
2630 sage: e0,e1,e2 = J.gens()
2646 We can change the generator prefix::
2648 sage: HadamardEJA(3, prefix='r').gens()
2652 def __init__(self
, n
, **kwargs
):
2654 jordan_product
= lambda x
,y
: x
2655 inner_product
= lambda x
,y
: x
2657 def jordan_product(x
,y
):
2659 return P( xi
*yi
for (xi
,yi
) in zip(x
,y
) )
2661 def inner_product(x
,y
):
2664 # New defaults for keyword arguments. Don't orthonormalize
2665 # because our basis is already orthonormal with respect to our
2666 # inner-product. Don't check the axioms, because we know this
2667 # is a valid EJA... but do double-check if the user passes
2668 # check_axioms=True. Note: we DON'T override the "check_field"
2669 # default here, because the user can pass in a field!
2670 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2671 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2673 column_basis
= tuple( b
.column() for b
in FreeModule(ZZ
, n
).basis() )
2674 super().__init
__(column_basis
,
2679 self
.rank
.set_cache(n
)
2682 self
.one
.set_cache( self
.zero() )
2684 self
.one
.set_cache( sum(self
.gens()) )
2687 def _max_random_instance_size():
2689 The maximum dimension of a random HadamardEJA.
2694 def random_instance(cls
, **kwargs
):
2696 Return a random instance of this type of algebra.
2698 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2699 return cls(n
, **kwargs
)
2702 class BilinearFormEJA(ConcreteEJA
):
2704 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2705 with the half-trace inner product and jordan product ``x*y =
2706 (<Bx,y>,y_bar>, x0*y_bar + y0*x_bar)`` where `B = 1 \times B22` is
2707 a symmetric positive-definite "bilinear form" matrix. Its
2708 dimension is the size of `B`, and it has rank two in dimensions
2709 larger than two. It reduces to the ``JordanSpinEJA`` when `B` is
2710 the identity matrix of order ``n``.
2712 We insist that the one-by-one upper-left identity block of `B` be
2713 passed in as well so that we can be passed a matrix of size zero
2714 to construct a trivial algebra.
2718 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
2719 ....: JordanSpinEJA)
2723 When no bilinear form is specified, the identity matrix is used,
2724 and the resulting algebra is the Jordan spin algebra::
2726 sage: B = matrix.identity(AA,3)
2727 sage: J0 = BilinearFormEJA(B)
2728 sage: J1 = JordanSpinEJA(3)
2729 sage: J0.multiplication_table() == J0.multiplication_table()
2732 An error is raised if the matrix `B` does not correspond to a
2733 positive-definite bilinear form::
2735 sage: B = matrix.random(QQ,2,3)
2736 sage: J = BilinearFormEJA(B)
2737 Traceback (most recent call last):
2739 ValueError: bilinear form is not positive-definite
2740 sage: B = matrix.zero(QQ,3)
2741 sage: J = BilinearFormEJA(B)
2742 Traceback (most recent call last):
2744 ValueError: bilinear form is not positive-definite
2748 We can create a zero-dimensional algebra::
2750 sage: B = matrix.identity(AA,0)
2751 sage: J = BilinearFormEJA(B)
2755 We can check the multiplication condition given in the Jordan, von
2756 Neumann, and Wigner paper (and also discussed on my "On the
2757 symmetry..." paper). Note that this relies heavily on the standard
2758 choice of basis, as does anything utilizing the bilinear form
2759 matrix. We opt not to orthonormalize the basis, because if we
2760 did, we would have to normalize the `s_{i}` in a similar manner::
2762 sage: set_random_seed()
2763 sage: n = ZZ.random_element(5)
2764 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2765 sage: B11 = matrix.identity(QQ,1)
2766 sage: B22 = M.transpose()*M
2767 sage: B = block_matrix(2,2,[ [B11,0 ],
2769 sage: J = BilinearFormEJA(B, orthonormalize=False)
2770 sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
2771 sage: V = J.vector_space()
2772 sage: sis = [ J( V([0] + (M.inverse()*ei).list()).column() )
2773 ....: for ei in eis ]
2774 sage: actual = [ sis[i]*sis[j]
2775 ....: for i in range(n-1)
2776 ....: for j in range(n-1) ]
2777 sage: expected = [ J.one() if i == j else J.zero()
2778 ....: for i in range(n-1)
2779 ....: for j in range(n-1) ]
2780 sage: actual == expected
2784 def __init__(self
, B
, **kwargs
):
2785 # The matrix "B" is supplied by the user in most cases,
2786 # so it makes sense to check whether or not its positive-
2787 # definite unless we are specifically asked not to...
2788 if ("check_axioms" not in kwargs
) or kwargs
["check_axioms"]:
2789 if not B
.is_positive_definite():
2790 raise ValueError("bilinear form is not positive-definite")
2792 # However, all of the other data for this EJA is computed
2793 # by us in manner that guarantees the axioms are
2794 # satisfied. So, again, unless we are specifically asked to
2795 # verify things, we'll skip the rest of the checks.
2796 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2798 def inner_product(x
,y
):
2799 return (y
.T
*B
*x
)[0,0]
2801 def jordan_product(x
,y
):
2807 z0
= inner_product(y
,x
)
2808 zbar
= y0
*xbar
+ x0
*ybar
2809 return P([z0
] + zbar
.list())
2812 column_basis
= tuple( b
.column() for b
in FreeModule(ZZ
, n
).basis() )
2814 # TODO: I haven't actually checked this, but it seems legit.
2819 super().__init
__(column_basis
,
2822 associative
=associative
,
2825 # The rank of this algebra is two, unless we're in a
2826 # one-dimensional ambient space (because the rank is bounded
2827 # by the ambient dimension).
2828 self
.rank
.set_cache(min(n
,2))
2831 self
.one
.set_cache( self
.zero() )
2833 self
.one
.set_cache( self
.monomial(0) )
2836 def _max_random_instance_size():
2838 The maximum dimension of a random BilinearFormEJA.
2843 def random_instance(cls
, **kwargs
):
2845 Return a random instance of this algebra.
2847 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2849 B
= matrix
.identity(ZZ
, n
)
2850 return cls(B
, **kwargs
)
2852 B11
= matrix
.identity(ZZ
, 1)
2853 M
= matrix
.random(ZZ
, n
-1)
2854 I
= matrix
.identity(ZZ
, n
-1)
2856 while alpha
.is_zero():
2857 alpha
= ZZ
.random_element().abs()
2858 B22
= M
.transpose()*M
+ alpha
*I
2860 from sage
.matrix
.special
import block_matrix
2861 B
= block_matrix(2,2, [ [B11
, ZZ(0) ],
2864 return cls(B
, **kwargs
)
2867 class JordanSpinEJA(BilinearFormEJA
):
2869 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2870 with the usual inner product and jordan product ``x*y =
2871 (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2876 sage: from mjo.eja.eja_algebra import JordanSpinEJA
2880 This multiplication table can be verified by hand::
2882 sage: J = JordanSpinEJA(4)
2883 sage: e0,e1,e2,e3 = J.gens()
2899 We can change the generator prefix::
2901 sage: JordanSpinEJA(2, prefix='B').gens()
2906 Ensure that we have the usual inner product on `R^n`::
2908 sage: set_random_seed()
2909 sage: J = JordanSpinEJA.random_instance()
2910 sage: x,y = J.random_elements(2)
2911 sage: actual = x.inner_product(y)
2912 sage: expected = x.to_vector().inner_product(y.to_vector())
2913 sage: actual == expected
2917 def __init__(self
, n
, **kwargs
):
2918 # This is a special case of the BilinearFormEJA with the
2919 # identity matrix as its bilinear form.
2920 B
= matrix
.identity(ZZ
, n
)
2922 # Don't orthonormalize because our basis is already
2923 # orthonormal with respect to our inner-product.
2924 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2926 # But also don't pass check_field=False here, because the user
2927 # can pass in a field!
2928 super().__init
__(B
, **kwargs
)
2931 def _max_random_instance_size():
2933 The maximum dimension of a random JordanSpinEJA.
2938 def random_instance(cls
, **kwargs
):
2940 Return a random instance of this type of algebra.
2942 Needed here to override the implementation for ``BilinearFormEJA``.
2944 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2945 return cls(n
, **kwargs
)
2948 class TrivialEJA(ConcreteEJA
):
2950 The trivial Euclidean Jordan algebra consisting of only a zero element.
2954 sage: from mjo.eja.eja_algebra import TrivialEJA
2958 sage: J = TrivialEJA()
2965 sage: 7*J.one()*12*J.one()
2967 sage: J.one().inner_product(J.one())
2969 sage: J.one().norm()
2971 sage: J.one().subalgebra_generated_by()
2972 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2977 def __init__(self
, **kwargs
):
2978 jordan_product
= lambda x
,y
: x
2979 inner_product
= lambda x
,y
: 0
2982 # New defaults for keyword arguments
2983 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2984 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2986 super().__init
__(basis
,
2992 # The rank is zero using my definition, namely the dimension of the
2993 # largest subalgebra generated by any element.
2994 self
.rank
.set_cache(0)
2995 self
.one
.set_cache( self
.zero() )
2998 def random_instance(cls
, **kwargs
):
2999 # We don't take a "size" argument so the superclass method is
3000 # inappropriate for us.
3001 return cls(**kwargs
)
3004 class CartesianProductEJA(FiniteDimensionalEJA
):
3006 The external (orthogonal) direct sum of two or more Euclidean
3007 Jordan algebras. Every Euclidean Jordan algebra decomposes into an
3008 orthogonal direct sum of simple Euclidean Jordan algebras which is
3009 then isometric to a Cartesian product, so no generality is lost by
3010 providing only this construction.
3014 sage: from mjo.eja.eja_algebra import (random_eja,
3015 ....: CartesianProductEJA,
3017 ....: JordanSpinEJA,
3018 ....: RealSymmetricEJA)
3022 The Jordan product is inherited from our factors and implemented by
3023 our CombinatorialFreeModule Cartesian product superclass::
3025 sage: set_random_seed()
3026 sage: J1 = HadamardEJA(2)
3027 sage: J2 = RealSymmetricEJA(2)
3028 sage: J = cartesian_product([J1,J2])
3029 sage: x,y = J.random_elements(2)
3033 The ability to retrieve the original factors is implemented by our
3034 CombinatorialFreeModule Cartesian product superclass::
3036 sage: J1 = HadamardEJA(2, field=QQ)
3037 sage: J2 = JordanSpinEJA(3, field=QQ)
3038 sage: J = cartesian_product([J1,J2])
3039 sage: J.cartesian_factors()
3040 (Euclidean Jordan algebra of dimension 2 over Rational Field,
3041 Euclidean Jordan algebra of dimension 3 over Rational Field)
3043 You can provide more than two factors::
3045 sage: J1 = HadamardEJA(2)
3046 sage: J2 = JordanSpinEJA(3)
3047 sage: J3 = RealSymmetricEJA(3)
3048 sage: cartesian_product([J1,J2,J3])
3049 Euclidean Jordan algebra of dimension 2 over Algebraic Real
3050 Field (+) Euclidean Jordan algebra of dimension 3 over Algebraic
3051 Real Field (+) Euclidean Jordan algebra of dimension 6 over
3052 Algebraic Real Field
3054 Rank is additive on a Cartesian product::
3056 sage: J1 = HadamardEJA(1)
3057 sage: J2 = RealSymmetricEJA(2)
3058 sage: J = cartesian_product([J1,J2])
3059 sage: J1.rank.clear_cache()
3060 sage: J2.rank.clear_cache()
3061 sage: J.rank.clear_cache()
3064 sage: J.rank() == J1.rank() + J2.rank()
3067 The same rank computation works over the rationals, with whatever
3070 sage: J1 = HadamardEJA(1, field=QQ, orthonormalize=False)
3071 sage: J2 = RealSymmetricEJA(2, field=QQ, orthonormalize=False)
3072 sage: J = cartesian_product([J1,J2])
3073 sage: J1.rank.clear_cache()
3074 sage: J2.rank.clear_cache()
3075 sage: J.rank.clear_cache()
3078 sage: J.rank() == J1.rank() + J2.rank()
3081 The product algebra will be associative if and only if all of its
3082 components are associative::
3084 sage: J1 = HadamardEJA(2)
3085 sage: J1.is_associative()
3087 sage: J2 = HadamardEJA(3)
3088 sage: J2.is_associative()
3090 sage: J3 = RealSymmetricEJA(3)
3091 sage: J3.is_associative()
3093 sage: CP1 = cartesian_product([J1,J2])
3094 sage: CP1.is_associative()
3096 sage: CP2 = cartesian_product([J1,J3])
3097 sage: CP2.is_associative()
3100 Cartesian products of Cartesian products work::
3102 sage: J1 = JordanSpinEJA(1)
3103 sage: J2 = JordanSpinEJA(1)
3104 sage: J3 = JordanSpinEJA(1)
3105 sage: J = cartesian_product([J1,cartesian_product([J2,J3])])
3106 sage: J.multiplication_table()
3107 +----++----+----+----+
3108 | * || e0 | e1 | e2 |
3109 +====++====+====+====+
3110 | e0 || e0 | 0 | 0 |
3111 +----++----+----+----+
3112 | e1 || 0 | e1 | 0 |
3113 +----++----+----+----+
3114 | e2 || 0 | 0 | e2 |
3115 +----++----+----+----+
3116 sage: HadamardEJA(3).multiplication_table()
3117 +----++----+----+----+
3118 | * || e0 | e1 | e2 |
3119 +====++====+====+====+
3120 | e0 || e0 | 0 | 0 |
3121 +----++----+----+----+
3122 | e1 || 0 | e1 | 0 |
3123 +----++----+----+----+
3124 | e2 || 0 | 0 | e2 |
3125 +----++----+----+----+
3129 All factors must share the same base field::
3131 sage: J1 = HadamardEJA(2, field=QQ)
3132 sage: J2 = RealSymmetricEJA(2)
3133 sage: CartesianProductEJA((J1,J2))
3134 Traceback (most recent call last):
3136 ValueError: all factors must share the same base field
3138 The cached unit element is the same one that would be computed::
3140 sage: set_random_seed() # long time
3141 sage: J1 = random_eja() # long time
3142 sage: J2 = random_eja() # long time
3143 sage: J = cartesian_product([J1,J2]) # long time
3144 sage: actual = J.one() # long time
3145 sage: J.one.clear_cache() # long time
3146 sage: expected = J.one() # long time
3147 sage: actual == expected # long time
3151 Element
= FiniteDimensionalEJAElement
3154 def __init__(self
, factors
, **kwargs
):
3159 self
._sets
= factors
3161 field
= factors
[0].base_ring()
3162 if not all( J
.base_ring() == field
for J
in factors
):
3163 raise ValueError("all factors must share the same base field")
3165 associative
= all( f
.is_associative() for f
in factors
)
3167 MS
= self
.matrix_space()
3171 for b
in factors
[i
].matrix_basis():
3176 basis
= tuple( MS(b
) for b
in basis
)
3178 # Define jordan/inner products that operate on that matrix_basis.
3179 def jordan_product(x
,y
):
3181 (factors
[i
](x
[i
])*factors
[i
](y
[i
])).to_matrix()
3185 def inner_product(x
, y
):
3187 factors
[i
](x
[i
]).inner_product(factors
[i
](y
[i
]))
3191 # There's no need to check the field since it already came
3192 # from an EJA. Likewise the axioms are guaranteed to be
3193 # satisfied, unless the guy writing this class sucks.
3195 # If you want the basis to be orthonormalized, orthonormalize
3197 FiniteDimensionalEJA
.__init
__(self
,
3202 orthonormalize
=False,
3203 associative
=associative
,
3204 cartesian_product
=True,
3208 ones
= tuple(J
.one().to_matrix() for J
in factors
)
3209 self
.one
.set_cache(self(ones
))
3210 self
.rank
.set_cache(sum(J
.rank() for J
in factors
))
3212 def cartesian_factors(self
):
3213 # Copy/pasted from CombinatorialFreeModule_CartesianProduct.
3216 def cartesian_factor(self
, i
):
3218 Return the ``i``th factor of this algebra.
3220 return self
._sets
[i
]
3223 # Copy/pasted from CombinatorialFreeModule_CartesianProduct.
3224 from sage
.categories
.cartesian_product
import cartesian_product
3225 return cartesian_product
.symbol
.join("%s" % factor
3226 for factor
in self
._sets
)
3228 def matrix_space(self
):
3230 Return the space that our matrix basis lives in as a Cartesian
3235 sage: from mjo.eja.eja_algebra import (HadamardEJA,
3236 ....: RealSymmetricEJA)
3240 sage: J1 = HadamardEJA(1)
3241 sage: J2 = RealSymmetricEJA(2)
3242 sage: J = cartesian_product([J1,J2])
3243 sage: J.matrix_space()
3244 The Cartesian product of (Full MatrixSpace of 1 by 1 dense
3245 matrices over Algebraic Real Field, Full MatrixSpace of 2
3246 by 2 dense matrices over Algebraic Real Field)
3249 from sage
.categories
.cartesian_product
import cartesian_product
3250 return cartesian_product( [J
.matrix_space()
3251 for J
in self
.cartesian_factors()] )
3254 def cartesian_projection(self
, i
):
3258 sage: from mjo.eja.eja_algebra import (random_eja,
3259 ....: JordanSpinEJA,
3261 ....: RealSymmetricEJA,
3262 ....: ComplexHermitianEJA)
3266 The projection morphisms are Euclidean Jordan algebra
3269 sage: J1 = HadamardEJA(2)
3270 sage: J2 = RealSymmetricEJA(2)
3271 sage: J = cartesian_product([J1,J2])
3272 sage: J.cartesian_projection(0)
3273 Linear operator between finite-dimensional Euclidean Jordan
3274 algebras represented by the matrix:
3277 Domain: Euclidean Jordan algebra of dimension 2 over Algebraic
3278 Real Field (+) Euclidean Jordan algebra of dimension 3 over
3279 Algebraic Real Field
3280 Codomain: Euclidean Jordan algebra of dimension 2 over Algebraic
3282 sage: J.cartesian_projection(1)
3283 Linear operator between finite-dimensional Euclidean Jordan
3284 algebras represented by the matrix:
3288 Domain: Euclidean Jordan algebra of dimension 2 over Algebraic
3289 Real Field (+) Euclidean Jordan algebra of dimension 3 over
3290 Algebraic Real Field
3291 Codomain: Euclidean Jordan algebra of dimension 3 over Algebraic
3294 The projections work the way you'd expect on the vector
3295 representation of an element::
3297 sage: J1 = JordanSpinEJA(2)
3298 sage: J2 = ComplexHermitianEJA(2)
3299 sage: J = cartesian_product([J1,J2])
3300 sage: pi_left = J.cartesian_projection(0)
3301 sage: pi_right = J.cartesian_projection(1)
3302 sage: pi_left(J.one()).to_vector()
3304 sage: pi_right(J.one()).to_vector()
3306 sage: J.one().to_vector()
3311 The answer never changes::
3313 sage: set_random_seed()
3314 sage: J1 = random_eja()
3315 sage: J2 = random_eja()
3316 sage: J = cartesian_product([J1,J2])
3317 sage: P0 = J.cartesian_projection(0)
3318 sage: P1 = J.cartesian_projection(0)
3323 offset
= sum( self
.cartesian_factor(k
).dimension()
3325 Ji
= self
.cartesian_factor(i
)
3326 Pi
= self
._module
_morphism
(lambda j
: Ji
.monomial(j
- offset
),
3329 return FiniteDimensionalEJAOperator(self
,Ji
,Pi
.matrix())
3332 def cartesian_embedding(self
, i
):
3336 sage: from mjo.eja.eja_algebra import (random_eja,
3337 ....: JordanSpinEJA,
3339 ....: RealSymmetricEJA)
3343 The embedding morphisms are Euclidean Jordan algebra
3346 sage: J1 = HadamardEJA(2)
3347 sage: J2 = RealSymmetricEJA(2)
3348 sage: J = cartesian_product([J1,J2])
3349 sage: J.cartesian_embedding(0)
3350 Linear operator between finite-dimensional Euclidean Jordan
3351 algebras represented by the matrix:
3357 Domain: Euclidean Jordan algebra of dimension 2 over
3358 Algebraic Real Field
3359 Codomain: Euclidean Jordan algebra of dimension 2 over
3360 Algebraic Real Field (+) Euclidean Jordan algebra of
3361 dimension 3 over Algebraic Real Field
3362 sage: J.cartesian_embedding(1)
3363 Linear operator between finite-dimensional Euclidean Jordan
3364 algebras represented by the matrix:
3370 Domain: Euclidean Jordan algebra of dimension 3 over
3371 Algebraic Real Field
3372 Codomain: Euclidean Jordan algebra of dimension 2 over
3373 Algebraic Real Field (+) Euclidean Jordan algebra of
3374 dimension 3 over Algebraic Real Field
3376 The embeddings work the way you'd expect on the vector
3377 representation of an element::
3379 sage: J1 = JordanSpinEJA(3)
3380 sage: J2 = RealSymmetricEJA(2)
3381 sage: J = cartesian_product([J1,J2])
3382 sage: iota_left = J.cartesian_embedding(0)
3383 sage: iota_right = J.cartesian_embedding(1)
3384 sage: iota_left(J1.zero()) == J.zero()
3386 sage: iota_right(J2.zero()) == J.zero()
3388 sage: J1.one().to_vector()
3390 sage: iota_left(J1.one()).to_vector()
3392 sage: J2.one().to_vector()
3394 sage: iota_right(J2.one()).to_vector()
3396 sage: J.one().to_vector()
3401 The answer never changes::
3403 sage: set_random_seed()
3404 sage: J1 = random_eja()
3405 sage: J2 = random_eja()
3406 sage: J = cartesian_product([J1,J2])
3407 sage: E0 = J.cartesian_embedding(0)
3408 sage: E1 = J.cartesian_embedding(0)
3412 Composing a projection with the corresponding inclusion should
3413 produce the identity map, and mismatching them should produce
3416 sage: set_random_seed()
3417 sage: J1 = random_eja()
3418 sage: J2 = random_eja()
3419 sage: J = cartesian_product([J1,J2])
3420 sage: iota_left = J.cartesian_embedding(0)
3421 sage: iota_right = J.cartesian_embedding(1)
3422 sage: pi_left = J.cartesian_projection(0)
3423 sage: pi_right = J.cartesian_projection(1)
3424 sage: pi_left*iota_left == J1.one().operator()
3426 sage: pi_right*iota_right == J2.one().operator()
3428 sage: (pi_left*iota_right).is_zero()
3430 sage: (pi_right*iota_left).is_zero()
3434 offset
= sum( self
.cartesian_factor(k
).dimension()
3436 Ji
= self
.cartesian_factor(i
)
3437 Ei
= Ji
._module
_morphism
(lambda j
: self
.monomial(j
+ offset
),
3439 return FiniteDimensionalEJAOperator(Ji
,self
,Ei
.matrix())
3443 FiniteDimensionalEJA
.CartesianProduct
= CartesianProductEJA
3445 class RationalBasisCartesianProductEJA(CartesianProductEJA
,
3448 A separate class for products of algebras for which we know a
3453 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
3454 ....: RealSymmetricEJA)
3458 This gives us fast characteristic polynomial computations in
3459 product algebras, too::
3462 sage: J1 = JordanSpinEJA(2)
3463 sage: J2 = RealSymmetricEJA(3)
3464 sage: J = cartesian_product([J1,J2])
3465 sage: J.characteristic_polynomial_of().degree()
3471 def __init__(self
, algebras
, **kwargs
):
3472 CartesianProductEJA
.__init
__(self
, algebras
, **kwargs
)
3474 self
._rational
_algebra
= None
3475 if self
.vector_space().base_field() is not QQ
:
3476 self
._rational
_algebra
= cartesian_product([
3477 r
._rational
_algebra
for r
in algebras
3481 RationalBasisEJA
.CartesianProduct
= RationalBasisCartesianProductEJA
3483 def random_eja(*args
, **kwargs
):
3484 J1
= ConcreteEJA
.random_instance(*args
, **kwargs
)
3486 # This might make Cartesian products appear roughly as often as
3487 # any other ConcreteEJA.
3488 if ZZ
.random_element(len(ConcreteEJA
.__subclasses
__()) + 1) == 0:
3489 # Use random_eja() again so we can get more than two factors.
3490 J2
= random_eja(*args
, **kwargs
)
3491 J
= cartesian_product([J1
,J2
])