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 # This does in fact compute the real part of the trace.
1841 # If we compute the trace of e.g. a complex matrix M,
1842 # then we do so by adding up its diagonal entries --
1843 # call them z_1 through z_n. The real embedding of z_1
1844 # will be a 2-by-2 REAL matrix [a, b; -b, a] whose trace
1845 # as a REAL matrix will be 2*a = 2*Re(z_1). And so forth.
1846 return (X
*Y
).trace()/cls
.dimension_over_reals()
1849 class RealMatrixEJA(MatrixEJA
):
1851 def dimension_over_reals():
1855 class RealSymmetricEJA(ConcreteEJA
, RealMatrixEJA
):
1857 The rank-n simple EJA consisting of real symmetric n-by-n
1858 matrices, the usual symmetric Jordan product, and the trace inner
1859 product. It has dimension `(n^2 + n)/2` over the reals.
1863 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1867 sage: J = RealSymmetricEJA(2)
1868 sage: e0, e1, e2 = J.gens()
1876 In theory, our "field" can be any subfield of the reals::
1878 sage: RealSymmetricEJA(2, field=RDF, check_axioms=True)
1879 Euclidean Jordan algebra of dimension 3 over Real Double Field
1880 sage: RealSymmetricEJA(2, field=RR, check_axioms=True)
1881 Euclidean Jordan algebra of dimension 3 over Real Field with
1882 53 bits of precision
1886 The dimension of this algebra is `(n^2 + n) / 2`::
1888 sage: set_random_seed()
1889 sage: n_max = RealSymmetricEJA._max_random_instance_size()
1890 sage: n = ZZ.random_element(1, n_max)
1891 sage: J = RealSymmetricEJA(n)
1892 sage: J.dimension() == (n^2 + n)/2
1895 The Jordan multiplication is what we think it is::
1897 sage: set_random_seed()
1898 sage: J = RealSymmetricEJA.random_instance()
1899 sage: x,y = J.random_elements(2)
1900 sage: actual = (x*y).to_matrix()
1901 sage: X = x.to_matrix()
1902 sage: Y = y.to_matrix()
1903 sage: expected = (X*Y + Y*X)/2
1904 sage: actual == expected
1906 sage: J(expected) == x*y
1909 We can change the generator prefix::
1911 sage: RealSymmetricEJA(3, prefix='q').gens()
1912 (q0, q1, q2, q3, q4, q5)
1914 We can construct the (trivial) algebra of rank zero::
1916 sage: RealSymmetricEJA(0)
1917 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1921 def _denormalized_basis(cls
, n
):
1923 Return a basis for the space of real symmetric n-by-n matrices.
1927 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1931 sage: set_random_seed()
1932 sage: n = ZZ.random_element(1,5)
1933 sage: B = RealSymmetricEJA._denormalized_basis(n)
1934 sage: all( M.is_symmetric() for M in B)
1938 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1942 for j
in range(i
+1):
1943 Eij
= matrix(ZZ
, n
, lambda k
,l
: k
==i
and l
==j
)
1947 Sij
= Eij
+ Eij
.transpose()
1953 def _max_random_instance_size():
1954 return 4 # Dimension 10
1957 def random_instance(cls
, **kwargs
):
1959 Return a random instance of this type of algebra.
1961 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1962 return cls(n
, **kwargs
)
1964 def __init__(self
, n
, **kwargs
):
1965 # We know this is a valid EJA, but will double-check
1966 # if the user passes check_axioms=True.
1967 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
1973 super().__init
__(self
._denormalized
_basis
(n
),
1974 self
.jordan_product
,
1975 self
.trace_inner_product
,
1976 associative
=associative
,
1979 # TODO: this could be factored out somehow, but is left here
1980 # because the MatrixEJA is not presently a subclass of the
1981 # FDEJA class that defines rank() and one().
1982 self
.rank
.set_cache(n
)
1983 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
1984 self
.one
.set_cache(self(idV
))
1988 class ComplexMatrixEJA(MatrixEJA
):
1989 # A manual dictionary-cache for the complex_extension() method,
1990 # since apparently @classmethods can't also be @cached_methods.
1991 _complex_extension
= {}
1994 def complex_extension(cls
,field
):
1996 The complex field that we embed/unembed, as an extension
1997 of the given ``field``.
1999 if field
in cls
._complex
_extension
:
2000 return cls
._complex
_extension
[field
]
2002 # Sage doesn't know how to adjoin the complex "i" (the root of
2003 # x^2 + 1) to a field in a general way. Here, we just enumerate
2004 # all of the cases that I have cared to support so far.
2006 # Sage doesn't know how to embed AA into QQbar, i.e. how
2007 # to adjoin sqrt(-1) to AA.
2009 elif not field
.is_exact():
2011 F
= field
.complex_field()
2013 # Works for QQ and... maybe some other fields.
2014 R
= PolynomialRing(field
, 'z')
2016 F
= field
.extension(z
**2 + 1, 'I', embedding
=CLF(-1).sqrt())
2018 cls
._complex
_extension
[field
] = F
2022 def dimension_over_reals():
2026 def real_embed(cls
,M
):
2028 Embed the n-by-n complex matrix ``M`` into the space of real
2029 matrices of size 2n-by-2n via the map the sends each entry `z = a +
2030 bi` to the block matrix ``[[a,b],[-b,a]]``.
2034 sage: from mjo.eja.eja_algebra import ComplexMatrixEJA
2038 sage: F = QuadraticField(-1, 'I')
2039 sage: x1 = F(4 - 2*i)
2040 sage: x2 = F(1 + 2*i)
2043 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
2044 sage: ComplexMatrixEJA.real_embed(M)
2053 Embedding is a homomorphism (isomorphism, in fact)::
2055 sage: set_random_seed()
2056 sage: n = ZZ.random_element(3)
2057 sage: F = QuadraticField(-1, 'I')
2058 sage: X = random_matrix(F, n)
2059 sage: Y = random_matrix(F, n)
2060 sage: Xe = ComplexMatrixEJA.real_embed(X)
2061 sage: Ye = ComplexMatrixEJA.real_embed(Y)
2062 sage: XYe = ComplexMatrixEJA.real_embed(X*Y)
2067 super().real_embed(M
)
2070 # We don't need any adjoined elements...
2071 field
= M
.base_ring().base_ring()
2077 blocks
.append(matrix(field
, 2, [ [ a
, b
],
2080 return matrix
.block(field
, n
, blocks
)
2084 def real_unembed(cls
,M
):
2086 The inverse of _embed_complex_matrix().
2090 sage: from mjo.eja.eja_algebra import ComplexMatrixEJA
2094 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
2095 ....: [-2, 1, -4, 3],
2096 ....: [ 9, 10, 11, 12],
2097 ....: [-10, 9, -12, 11] ])
2098 sage: ComplexMatrixEJA.real_unembed(A)
2100 [ 10*I + 9 12*I + 11]
2104 Unembedding is the inverse of embedding::
2106 sage: set_random_seed()
2107 sage: F = QuadraticField(-1, 'I')
2108 sage: M = random_matrix(F, 3)
2109 sage: Me = ComplexMatrixEJA.real_embed(M)
2110 sage: ComplexMatrixEJA.real_unembed(Me) == M
2114 super().real_unembed(M
)
2116 d
= cls
.dimension_over_reals()
2117 F
= cls
.complex_extension(M
.base_ring())
2120 # Go top-left to bottom-right (reading order), converting every
2121 # 2-by-2 block we see to a single complex element.
2123 for k
in range(n
/d
):
2124 for j
in range(n
/d
):
2125 submat
= M
[d
*k
:d
*k
+d
,d
*j
:d
*j
+d
]
2126 if submat
[0,0] != submat
[1,1]:
2127 raise ValueError('bad on-diagonal submatrix')
2128 if submat
[0,1] != -submat
[1,0]:
2129 raise ValueError('bad off-diagonal submatrix')
2130 z
= submat
[0,0] + submat
[0,1]*i
2133 return matrix(F
, n
/d
, elements
)
2136 class ComplexHermitianEJA(ConcreteEJA
, ComplexMatrixEJA
):
2138 The rank-n simple EJA consisting of complex Hermitian n-by-n
2139 matrices over the real numbers, the usual symmetric Jordan product,
2140 and the real-part-of-trace inner product. It has dimension `n^2` over
2145 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
2149 In theory, our "field" can be any subfield of the reals::
2151 sage: ComplexHermitianEJA(2, field=RDF, check_axioms=True)
2152 Euclidean Jordan algebra of dimension 4 over Real Double Field
2153 sage: ComplexHermitianEJA(2, field=RR, check_axioms=True)
2154 Euclidean Jordan algebra of dimension 4 over Real Field with
2155 53 bits of precision
2159 The dimension of this algebra is `n^2`::
2161 sage: set_random_seed()
2162 sage: n_max = ComplexHermitianEJA._max_random_instance_size()
2163 sage: n = ZZ.random_element(1, n_max)
2164 sage: J = ComplexHermitianEJA(n)
2165 sage: J.dimension() == n^2
2168 The Jordan multiplication is what we think it is::
2170 sage: set_random_seed()
2171 sage: J = ComplexHermitianEJA.random_instance()
2172 sage: x,y = J.random_elements(2)
2173 sage: actual = (x*y).to_matrix()
2174 sage: X = x.to_matrix()
2175 sage: Y = y.to_matrix()
2176 sage: expected = (X*Y + Y*X)/2
2177 sage: actual == expected
2179 sage: J(expected) == x*y
2182 We can change the generator prefix::
2184 sage: ComplexHermitianEJA(2, prefix='z').gens()
2187 We can construct the (trivial) algebra of rank zero::
2189 sage: ComplexHermitianEJA(0)
2190 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2195 def _denormalized_basis(cls
, n
):
2197 Returns a basis for the space of complex Hermitian n-by-n matrices.
2199 Why do we embed these? Basically, because all of numerical linear
2200 algebra assumes that you're working with vectors consisting of `n`
2201 entries from a field and scalars from the same field. There's no way
2202 to tell SageMath that (for example) the vectors contain complex
2203 numbers, while the scalar field is real.
2207 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
2211 sage: set_random_seed()
2212 sage: n = ZZ.random_element(1,5)
2213 sage: B = ComplexHermitianEJA._denormalized_basis(n)
2214 sage: all( M.is_symmetric() for M in B)
2219 R
= PolynomialRing(field
, 'z')
2221 F
= field
.extension(z
**2 + 1, 'I')
2224 # This is like the symmetric case, but we need to be careful:
2226 # * We want conjugate-symmetry, not just symmetry.
2227 # * The diagonal will (as a result) be real.
2230 Eij
= matrix
.zero(F
,n
)
2232 for j
in range(i
+1):
2236 Sij
= cls
.real_embed(Eij
)
2239 # The second one has a minus because it's conjugated.
2240 Eij
[j
,i
] = 1 # Eij = Eij + Eij.transpose()
2241 Sij_real
= cls
.real_embed(Eij
)
2243 # Eij = I*Eij - I*Eij.transpose()
2246 Sij_imag
= cls
.real_embed(Eij
)
2252 # Since we embedded these, we can drop back to the "field" that we
2253 # started with instead of the complex extension "F".
2254 return tuple( s
.change_ring(field
) for s
in S
)
2257 def __init__(self
, n
, **kwargs
):
2258 # We know this is a valid EJA, but will double-check
2259 # if the user passes check_axioms=True.
2260 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2266 super().__init
__(self
._denormalized
_basis
(n
),
2267 self
.jordan_product
,
2268 self
.trace_inner_product
,
2269 associative
=associative
,
2271 # TODO: this could be factored out somehow, but is left here
2272 # because the MatrixEJA is not presently a subclass of the
2273 # FDEJA class that defines rank() and one().
2274 self
.rank
.set_cache(n
)
2275 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
2276 self
.one
.set_cache(self(idV
))
2279 def _max_random_instance_size():
2280 return 3 # Dimension 9
2283 def random_instance(cls
, **kwargs
):
2285 Return a random instance of this type of algebra.
2287 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2288 return cls(n
, **kwargs
)
2290 class QuaternionMatrixEJA(MatrixEJA
):
2292 # A manual dictionary-cache for the quaternion_extension() method,
2293 # since apparently @classmethods can't also be @cached_methods.
2294 _quaternion_extension
= {}
2297 def quaternion_extension(cls
,field
):
2299 The quaternion field that we embed/unembed, as an extension
2300 of the given ``field``.
2302 if field
in cls
._quaternion
_extension
:
2303 return cls
._quaternion
_extension
[field
]
2305 Q
= QuaternionAlgebra(field
,-1,-1)
2307 cls
._quaternion
_extension
[field
] = Q
2311 def dimension_over_reals():
2315 def real_embed(cls
,M
):
2317 Embed the n-by-n quaternion matrix ``M`` into the space of real
2318 matrices of size 4n-by-4n by first sending each quaternion entry `z
2319 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
2320 c+di],[-c + di, a-bi]]`, and then embedding those into a real
2325 sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA
2329 sage: Q = QuaternionAlgebra(QQ,-1,-1)
2330 sage: i,j,k = Q.gens()
2331 sage: x = 1 + 2*i + 3*j + 4*k
2332 sage: M = matrix(Q, 1, [[x]])
2333 sage: QuaternionMatrixEJA.real_embed(M)
2339 Embedding is a homomorphism (isomorphism, in fact)::
2341 sage: set_random_seed()
2342 sage: n = ZZ.random_element(2)
2343 sage: Q = QuaternionAlgebra(QQ,-1,-1)
2344 sage: X = random_matrix(Q, n)
2345 sage: Y = random_matrix(Q, n)
2346 sage: Xe = QuaternionMatrixEJA.real_embed(X)
2347 sage: Ye = QuaternionMatrixEJA.real_embed(Y)
2348 sage: XYe = QuaternionMatrixEJA.real_embed(X*Y)
2353 super().real_embed(M
)
2354 quaternions
= M
.base_ring()
2357 F
= QuadraticField(-1, 'I')
2362 t
= z
.coefficient_tuple()
2367 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
2368 [-c
+ d
*i
, a
- b
*i
]])
2369 realM
= ComplexMatrixEJA
.real_embed(cplxM
)
2370 blocks
.append(realM
)
2372 # We should have real entries by now, so use the realest field
2373 # we've got for the return value.
2374 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
2379 def real_unembed(cls
,M
):
2381 The inverse of _embed_quaternion_matrix().
2385 sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA
2389 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
2390 ....: [-2, 1, -4, 3],
2391 ....: [-3, 4, 1, -2],
2392 ....: [-4, -3, 2, 1]])
2393 sage: QuaternionMatrixEJA.real_unembed(M)
2394 [1 + 2*i + 3*j + 4*k]
2398 Unembedding is the inverse of embedding::
2400 sage: set_random_seed()
2401 sage: Q = QuaternionAlgebra(QQ, -1, -1)
2402 sage: M = random_matrix(Q, 3)
2403 sage: Me = QuaternionMatrixEJA.real_embed(M)
2404 sage: QuaternionMatrixEJA.real_unembed(Me) == M
2408 super().real_unembed(M
)
2410 d
= cls
.dimension_over_reals()
2412 # Use the base ring of the matrix to ensure that its entries can be
2413 # multiplied by elements of the quaternion algebra.
2414 Q
= cls
.quaternion_extension(M
.base_ring())
2417 # Go top-left to bottom-right (reading order), converting every
2418 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
2421 for l
in range(n
/d
):
2422 for m
in range(n
/d
):
2423 submat
= ComplexMatrixEJA
.real_unembed(
2424 M
[d
*l
:d
*l
+d
,d
*m
:d
*m
+d
] )
2425 if submat
[0,0] != submat
[1,1].conjugate():
2426 raise ValueError('bad on-diagonal submatrix')
2427 if submat
[0,1] != -submat
[1,0].conjugate():
2428 raise ValueError('bad off-diagonal submatrix')
2429 z
= submat
[0,0].real()
2430 z
+= submat
[0,0].imag()*i
2431 z
+= submat
[0,1].real()*j
2432 z
+= submat
[0,1].imag()*k
2435 return matrix(Q
, n
/d
, elements
)
2438 class QuaternionHermitianEJA(ConcreteEJA
, QuaternionMatrixEJA
):
2440 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
2441 matrices, the usual symmetric Jordan product, and the
2442 real-part-of-trace inner product. It has dimension `2n^2 - n` over
2447 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2451 In theory, our "field" can be any subfield of the reals::
2453 sage: QuaternionHermitianEJA(2, field=RDF, check_axioms=True)
2454 Euclidean Jordan algebra of dimension 6 over Real Double Field
2455 sage: QuaternionHermitianEJA(2, field=RR, check_axioms=True)
2456 Euclidean Jordan algebra of dimension 6 over Real Field with
2457 53 bits of precision
2461 The dimension of this algebra is `2*n^2 - n`::
2463 sage: set_random_seed()
2464 sage: n_max = QuaternionHermitianEJA._max_random_instance_size()
2465 sage: n = ZZ.random_element(1, n_max)
2466 sage: J = QuaternionHermitianEJA(n)
2467 sage: J.dimension() == 2*(n^2) - n
2470 The Jordan multiplication is what we think it is::
2472 sage: set_random_seed()
2473 sage: J = QuaternionHermitianEJA.random_instance()
2474 sage: x,y = J.random_elements(2)
2475 sage: actual = (x*y).to_matrix()
2476 sage: X = x.to_matrix()
2477 sage: Y = y.to_matrix()
2478 sage: expected = (X*Y + Y*X)/2
2479 sage: actual == expected
2481 sage: J(expected) == x*y
2484 We can change the generator prefix::
2486 sage: QuaternionHermitianEJA(2, prefix='a').gens()
2487 (a0, a1, a2, a3, a4, a5)
2489 We can construct the (trivial) algebra of rank zero::
2491 sage: QuaternionHermitianEJA(0)
2492 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2496 def _denormalized_basis(cls
, n
):
2498 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
2500 Why do we embed these? Basically, because all of numerical
2501 linear algebra assumes that you're working with vectors consisting
2502 of `n` entries from a field and scalars from the same field. There's
2503 no way to tell SageMath that (for example) the vectors contain
2504 complex numbers, while the scalar field is real.
2508 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2512 sage: set_random_seed()
2513 sage: n = ZZ.random_element(1,5)
2514 sage: B = QuaternionHermitianEJA._denormalized_basis(n)
2515 sage: all( M.is_symmetric() for M in B )
2520 Q
= QuaternionAlgebra(QQ
,-1,-1)
2523 # This is like the symmetric case, but we need to be careful:
2525 # * We want conjugate-symmetry, not just symmetry.
2526 # * The diagonal will (as a result) be real.
2529 Eij
= matrix
.zero(Q
,n
)
2531 for j
in range(i
+1):
2535 Sij
= cls
.real_embed(Eij
)
2538 # The second, third, and fourth ones have a minus
2539 # because they're conjugated.
2540 # Eij = Eij + Eij.transpose()
2542 Sij_real
= cls
.real_embed(Eij
)
2544 # Eij = I*(Eij - Eij.transpose())
2547 Sij_I
= cls
.real_embed(Eij
)
2549 # Eij = J*(Eij - Eij.transpose())
2552 Sij_J
= cls
.real_embed(Eij
)
2554 # Eij = K*(Eij - Eij.transpose())
2557 Sij_K
= cls
.real_embed(Eij
)
2563 # Since we embedded these, we can drop back to the "field" that we
2564 # started with instead of the quaternion algebra "Q".
2565 return tuple( s
.change_ring(field
) for s
in S
)
2568 def __init__(self
, n
, **kwargs
):
2569 # We know this is a valid EJA, but will double-check
2570 # if the user passes check_axioms=True.
2571 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2577 super().__init
__(self
._denormalized
_basis
(n
),
2578 self
.jordan_product
,
2579 self
.trace_inner_product
,
2580 associative
=associative
,
2583 # TODO: this could be factored out somehow, but is left here
2584 # because the MatrixEJA is not presently a subclass of the
2585 # FDEJA class that defines rank() and one().
2586 self
.rank
.set_cache(n
)
2587 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
2588 self
.one
.set_cache(self(idV
))
2592 def _max_random_instance_size():
2594 The maximum rank of a random QuaternionHermitianEJA.
2596 return 2 # Dimension 6
2599 def random_instance(cls
, **kwargs
):
2601 Return a random instance of this type of algebra.
2603 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2604 return cls(n
, **kwargs
)
2607 class HadamardEJA(ConcreteEJA
):
2609 Return the Euclidean Jordan Algebra corresponding to the set
2610 `R^n` under the Hadamard product.
2612 Note: this is nothing more than the Cartesian product of ``n``
2613 copies of the spin algebra. Once Cartesian product algebras
2614 are implemented, this can go.
2618 sage: from mjo.eja.eja_algebra import HadamardEJA
2622 This multiplication table can be verified by hand::
2624 sage: J = HadamardEJA(3)
2625 sage: e0,e1,e2 = J.gens()
2641 We can change the generator prefix::
2643 sage: HadamardEJA(3, prefix='r').gens()
2647 def __init__(self
, n
, **kwargs
):
2649 jordan_product
= lambda x
,y
: x
2650 inner_product
= lambda x
,y
: x
2652 def jordan_product(x
,y
):
2654 return P( xi
*yi
for (xi
,yi
) in zip(x
,y
) )
2656 def inner_product(x
,y
):
2659 # New defaults for keyword arguments. Don't orthonormalize
2660 # because our basis is already orthonormal with respect to our
2661 # inner-product. Don't check the axioms, because we know this
2662 # is a valid EJA... but do double-check if the user passes
2663 # check_axioms=True. Note: we DON'T override the "check_field"
2664 # default here, because the user can pass in a field!
2665 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2666 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2668 column_basis
= tuple( b
.column() for b
in FreeModule(ZZ
, n
).basis() )
2669 super().__init
__(column_basis
,
2674 self
.rank
.set_cache(n
)
2677 self
.one
.set_cache( self
.zero() )
2679 self
.one
.set_cache( sum(self
.gens()) )
2682 def _max_random_instance_size():
2684 The maximum dimension of a random HadamardEJA.
2689 def random_instance(cls
, **kwargs
):
2691 Return a random instance of this type of algebra.
2693 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2694 return cls(n
, **kwargs
)
2697 class BilinearFormEJA(ConcreteEJA
):
2699 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2700 with the half-trace inner product and jordan product ``x*y =
2701 (<Bx,y>,y_bar>, x0*y_bar + y0*x_bar)`` where `B = 1 \times B22` is
2702 a symmetric positive-definite "bilinear form" matrix. Its
2703 dimension is the size of `B`, and it has rank two in dimensions
2704 larger than two. It reduces to the ``JordanSpinEJA`` when `B` is
2705 the identity matrix of order ``n``.
2707 We insist that the one-by-one upper-left identity block of `B` be
2708 passed in as well so that we can be passed a matrix of size zero
2709 to construct a trivial algebra.
2713 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
2714 ....: JordanSpinEJA)
2718 When no bilinear form is specified, the identity matrix is used,
2719 and the resulting algebra is the Jordan spin algebra::
2721 sage: B = matrix.identity(AA,3)
2722 sage: J0 = BilinearFormEJA(B)
2723 sage: J1 = JordanSpinEJA(3)
2724 sage: J0.multiplication_table() == J0.multiplication_table()
2727 An error is raised if the matrix `B` does not correspond to a
2728 positive-definite bilinear form::
2730 sage: B = matrix.random(QQ,2,3)
2731 sage: J = BilinearFormEJA(B)
2732 Traceback (most recent call last):
2734 ValueError: bilinear form is not positive-definite
2735 sage: B = matrix.zero(QQ,3)
2736 sage: J = BilinearFormEJA(B)
2737 Traceback (most recent call last):
2739 ValueError: bilinear form is not positive-definite
2743 We can create a zero-dimensional algebra::
2745 sage: B = matrix.identity(AA,0)
2746 sage: J = BilinearFormEJA(B)
2750 We can check the multiplication condition given in the Jordan, von
2751 Neumann, and Wigner paper (and also discussed on my "On the
2752 symmetry..." paper). Note that this relies heavily on the standard
2753 choice of basis, as does anything utilizing the bilinear form
2754 matrix. We opt not to orthonormalize the basis, because if we
2755 did, we would have to normalize the `s_{i}` in a similar manner::
2757 sage: set_random_seed()
2758 sage: n = ZZ.random_element(5)
2759 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2760 sage: B11 = matrix.identity(QQ,1)
2761 sage: B22 = M.transpose()*M
2762 sage: B = block_matrix(2,2,[ [B11,0 ],
2764 sage: J = BilinearFormEJA(B, orthonormalize=False)
2765 sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
2766 sage: V = J.vector_space()
2767 sage: sis = [ J( V([0] + (M.inverse()*ei).list()).column() )
2768 ....: for ei in eis ]
2769 sage: actual = [ sis[i]*sis[j]
2770 ....: for i in range(n-1)
2771 ....: for j in range(n-1) ]
2772 sage: expected = [ J.one() if i == j else J.zero()
2773 ....: for i in range(n-1)
2774 ....: for j in range(n-1) ]
2775 sage: actual == expected
2779 def __init__(self
, B
, **kwargs
):
2780 # The matrix "B" is supplied by the user in most cases,
2781 # so it makes sense to check whether or not its positive-
2782 # definite unless we are specifically asked not to...
2783 if ("check_axioms" not in kwargs
) or kwargs
["check_axioms"]:
2784 if not B
.is_positive_definite():
2785 raise ValueError("bilinear form is not positive-definite")
2787 # However, all of the other data for this EJA is computed
2788 # by us in manner that guarantees the axioms are
2789 # satisfied. So, again, unless we are specifically asked to
2790 # verify things, we'll skip the rest of the checks.
2791 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2793 def inner_product(x
,y
):
2794 return (y
.T
*B
*x
)[0,0]
2796 def jordan_product(x
,y
):
2802 z0
= inner_product(y
,x
)
2803 zbar
= y0
*xbar
+ x0
*ybar
2804 return P([z0
] + zbar
.list())
2807 column_basis
= tuple( b
.column() for b
in FreeModule(ZZ
, n
).basis() )
2809 # TODO: I haven't actually checked this, but it seems legit.
2814 super().__init
__(column_basis
,
2817 associative
=associative
,
2820 # The rank of this algebra is two, unless we're in a
2821 # one-dimensional ambient space (because the rank is bounded
2822 # by the ambient dimension).
2823 self
.rank
.set_cache(min(n
,2))
2826 self
.one
.set_cache( self
.zero() )
2828 self
.one
.set_cache( self
.monomial(0) )
2831 def _max_random_instance_size():
2833 The maximum dimension of a random BilinearFormEJA.
2838 def random_instance(cls
, **kwargs
):
2840 Return a random instance of this algebra.
2842 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2844 B
= matrix
.identity(ZZ
, n
)
2845 return cls(B
, **kwargs
)
2847 B11
= matrix
.identity(ZZ
, 1)
2848 M
= matrix
.random(ZZ
, n
-1)
2849 I
= matrix
.identity(ZZ
, n
-1)
2851 while alpha
.is_zero():
2852 alpha
= ZZ
.random_element().abs()
2853 B22
= M
.transpose()*M
+ alpha
*I
2855 from sage
.matrix
.special
import block_matrix
2856 B
= block_matrix(2,2, [ [B11
, ZZ(0) ],
2859 return cls(B
, **kwargs
)
2862 class JordanSpinEJA(BilinearFormEJA
):
2864 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2865 with the usual inner product and jordan product ``x*y =
2866 (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2871 sage: from mjo.eja.eja_algebra import JordanSpinEJA
2875 This multiplication table can be verified by hand::
2877 sage: J = JordanSpinEJA(4)
2878 sage: e0,e1,e2,e3 = J.gens()
2894 We can change the generator prefix::
2896 sage: JordanSpinEJA(2, prefix='B').gens()
2901 Ensure that we have the usual inner product on `R^n`::
2903 sage: set_random_seed()
2904 sage: J = JordanSpinEJA.random_instance()
2905 sage: x,y = J.random_elements(2)
2906 sage: actual = x.inner_product(y)
2907 sage: expected = x.to_vector().inner_product(y.to_vector())
2908 sage: actual == expected
2912 def __init__(self
, n
, **kwargs
):
2913 # This is a special case of the BilinearFormEJA with the
2914 # identity matrix as its bilinear form.
2915 B
= matrix
.identity(ZZ
, n
)
2917 # Don't orthonormalize because our basis is already
2918 # orthonormal with respect to our inner-product.
2919 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2921 # But also don't pass check_field=False here, because the user
2922 # can pass in a field!
2923 super().__init
__(B
, **kwargs
)
2926 def _max_random_instance_size():
2928 The maximum dimension of a random JordanSpinEJA.
2933 def random_instance(cls
, **kwargs
):
2935 Return a random instance of this type of algebra.
2937 Needed here to override the implementation for ``BilinearFormEJA``.
2939 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2940 return cls(n
, **kwargs
)
2943 class TrivialEJA(ConcreteEJA
):
2945 The trivial Euclidean Jordan algebra consisting of only a zero element.
2949 sage: from mjo.eja.eja_algebra import TrivialEJA
2953 sage: J = TrivialEJA()
2960 sage: 7*J.one()*12*J.one()
2962 sage: J.one().inner_product(J.one())
2964 sage: J.one().norm()
2966 sage: J.one().subalgebra_generated_by()
2967 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2972 def __init__(self
, **kwargs
):
2973 jordan_product
= lambda x
,y
: x
2974 inner_product
= lambda x
,y
: 0
2977 # New defaults for keyword arguments
2978 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2979 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2981 super().__init
__(basis
,
2987 # The rank is zero using my definition, namely the dimension of the
2988 # largest subalgebra generated by any element.
2989 self
.rank
.set_cache(0)
2990 self
.one
.set_cache( self
.zero() )
2993 def random_instance(cls
, **kwargs
):
2994 # We don't take a "size" argument so the superclass method is
2995 # inappropriate for us.
2996 return cls(**kwargs
)
2999 class CartesianProductEJA(FiniteDimensionalEJA
):
3001 The external (orthogonal) direct sum of two or more Euclidean
3002 Jordan algebras. Every Euclidean Jordan algebra decomposes into an
3003 orthogonal direct sum of simple Euclidean Jordan algebras which is
3004 then isometric to a Cartesian product, so no generality is lost by
3005 providing only this construction.
3009 sage: from mjo.eja.eja_algebra import (random_eja,
3010 ....: CartesianProductEJA,
3012 ....: JordanSpinEJA,
3013 ....: RealSymmetricEJA)
3017 The Jordan product is inherited from our factors and implemented by
3018 our CombinatorialFreeModule Cartesian product superclass::
3020 sage: set_random_seed()
3021 sage: J1 = HadamardEJA(2)
3022 sage: J2 = RealSymmetricEJA(2)
3023 sage: J = cartesian_product([J1,J2])
3024 sage: x,y = J.random_elements(2)
3028 The ability to retrieve the original factors is implemented by our
3029 CombinatorialFreeModule Cartesian product superclass::
3031 sage: J1 = HadamardEJA(2, field=QQ)
3032 sage: J2 = JordanSpinEJA(3, field=QQ)
3033 sage: J = cartesian_product([J1,J2])
3034 sage: J.cartesian_factors()
3035 (Euclidean Jordan algebra of dimension 2 over Rational Field,
3036 Euclidean Jordan algebra of dimension 3 over Rational Field)
3038 You can provide more than two factors::
3040 sage: J1 = HadamardEJA(2)
3041 sage: J2 = JordanSpinEJA(3)
3042 sage: J3 = RealSymmetricEJA(3)
3043 sage: cartesian_product([J1,J2,J3])
3044 Euclidean Jordan algebra of dimension 2 over Algebraic Real
3045 Field (+) Euclidean Jordan algebra of dimension 3 over Algebraic
3046 Real Field (+) Euclidean Jordan algebra of dimension 6 over
3047 Algebraic Real Field
3049 Rank is additive on a Cartesian product::
3051 sage: J1 = HadamardEJA(1)
3052 sage: J2 = RealSymmetricEJA(2)
3053 sage: J = cartesian_product([J1,J2])
3054 sage: J1.rank.clear_cache()
3055 sage: J2.rank.clear_cache()
3056 sage: J.rank.clear_cache()
3059 sage: J.rank() == J1.rank() + J2.rank()
3062 The same rank computation works over the rationals, with whatever
3065 sage: J1 = HadamardEJA(1, field=QQ, orthonormalize=False)
3066 sage: J2 = RealSymmetricEJA(2, field=QQ, orthonormalize=False)
3067 sage: J = cartesian_product([J1,J2])
3068 sage: J1.rank.clear_cache()
3069 sage: J2.rank.clear_cache()
3070 sage: J.rank.clear_cache()
3073 sage: J.rank() == J1.rank() + J2.rank()
3076 The product algebra will be associative if and only if all of its
3077 components are associative::
3079 sage: J1 = HadamardEJA(2)
3080 sage: J1.is_associative()
3082 sage: J2 = HadamardEJA(3)
3083 sage: J2.is_associative()
3085 sage: J3 = RealSymmetricEJA(3)
3086 sage: J3.is_associative()
3088 sage: CP1 = cartesian_product([J1,J2])
3089 sage: CP1.is_associative()
3091 sage: CP2 = cartesian_product([J1,J3])
3092 sage: CP2.is_associative()
3095 Cartesian products of Cartesian products work::
3097 sage: J1 = JordanSpinEJA(1)
3098 sage: J2 = JordanSpinEJA(1)
3099 sage: J3 = JordanSpinEJA(1)
3100 sage: J = cartesian_product([J1,cartesian_product([J2,J3])])
3101 sage: J.multiplication_table()
3102 +----++----+----+----+
3103 | * || e0 | e1 | e2 |
3104 +====++====+====+====+
3105 | e0 || e0 | 0 | 0 |
3106 +----++----+----+----+
3107 | e1 || 0 | e1 | 0 |
3108 +----++----+----+----+
3109 | e2 || 0 | 0 | e2 |
3110 +----++----+----+----+
3111 sage: HadamardEJA(3).multiplication_table()
3112 +----++----+----+----+
3113 | * || e0 | e1 | e2 |
3114 +====++====+====+====+
3115 | e0 || e0 | 0 | 0 |
3116 +----++----+----+----+
3117 | e1 || 0 | e1 | 0 |
3118 +----++----+----+----+
3119 | e2 || 0 | 0 | e2 |
3120 +----++----+----+----+
3124 All factors must share the same base field::
3126 sage: J1 = HadamardEJA(2, field=QQ)
3127 sage: J2 = RealSymmetricEJA(2)
3128 sage: CartesianProductEJA((J1,J2))
3129 Traceback (most recent call last):
3131 ValueError: all factors must share the same base field
3133 The cached unit element is the same one that would be computed::
3135 sage: set_random_seed() # long time
3136 sage: J1 = random_eja() # long time
3137 sage: J2 = random_eja() # long time
3138 sage: J = cartesian_product([J1,J2]) # long time
3139 sage: actual = J.one() # long time
3140 sage: J.one.clear_cache() # long time
3141 sage: expected = J.one() # long time
3142 sage: actual == expected # long time
3146 Element
= FiniteDimensionalEJAElement
3149 def __init__(self
, factors
, **kwargs
):
3154 self
._sets
= factors
3156 field
= factors
[0].base_ring()
3157 if not all( J
.base_ring() == field
for J
in factors
):
3158 raise ValueError("all factors must share the same base field")
3160 associative
= all( f
.is_associative() for f
in factors
)
3162 MS
= self
.matrix_space()
3166 for b
in factors
[i
].matrix_basis():
3171 basis
= tuple( MS(b
) for b
in basis
)
3173 # Define jordan/inner products that operate on that matrix_basis.
3174 def jordan_product(x
,y
):
3176 (factors
[i
](x
[i
])*factors
[i
](y
[i
])).to_matrix()
3180 def inner_product(x
, y
):
3182 factors
[i
](x
[i
]).inner_product(factors
[i
](y
[i
]))
3186 # There's no need to check the field since it already came
3187 # from an EJA. Likewise the axioms are guaranteed to be
3188 # satisfied, unless the guy writing this class sucks.
3190 # If you want the basis to be orthonormalized, orthonormalize
3192 FiniteDimensionalEJA
.__init
__(self
,
3197 orthonormalize
=False,
3198 associative
=associative
,
3199 cartesian_product
=True,
3203 ones
= tuple(J
.one().to_matrix() for J
in factors
)
3204 self
.one
.set_cache(self(ones
))
3205 self
.rank
.set_cache(sum(J
.rank() for J
in factors
))
3207 def cartesian_factors(self
):
3208 # Copy/pasted from CombinatorialFreeModule_CartesianProduct.
3211 def cartesian_factor(self
, i
):
3213 Return the ``i``th factor of this algebra.
3215 return self
._sets
[i
]
3218 # Copy/pasted from CombinatorialFreeModule_CartesianProduct.
3219 from sage
.categories
.cartesian_product
import cartesian_product
3220 return cartesian_product
.symbol
.join("%s" % factor
3221 for factor
in self
._sets
)
3223 def matrix_space(self
):
3225 Return the space that our matrix basis lives in as a Cartesian
3230 sage: from mjo.eja.eja_algebra import (HadamardEJA,
3231 ....: RealSymmetricEJA)
3235 sage: J1 = HadamardEJA(1)
3236 sage: J2 = RealSymmetricEJA(2)
3237 sage: J = cartesian_product([J1,J2])
3238 sage: J.matrix_space()
3239 The Cartesian product of (Full MatrixSpace of 1 by 1 dense
3240 matrices over Algebraic Real Field, Full MatrixSpace of 2
3241 by 2 dense matrices over Algebraic Real Field)
3244 from sage
.categories
.cartesian_product
import cartesian_product
3245 return cartesian_product( [J
.matrix_space()
3246 for J
in self
.cartesian_factors()] )
3249 def cartesian_projection(self
, i
):
3253 sage: from mjo.eja.eja_algebra import (random_eja,
3254 ....: JordanSpinEJA,
3256 ....: RealSymmetricEJA,
3257 ....: ComplexHermitianEJA)
3261 The projection morphisms are Euclidean Jordan algebra
3264 sage: J1 = HadamardEJA(2)
3265 sage: J2 = RealSymmetricEJA(2)
3266 sage: J = cartesian_product([J1,J2])
3267 sage: J.cartesian_projection(0)
3268 Linear operator between finite-dimensional Euclidean Jordan
3269 algebras represented by the matrix:
3272 Domain: Euclidean Jordan algebra of dimension 2 over Algebraic
3273 Real Field (+) Euclidean Jordan algebra of dimension 3 over
3274 Algebraic Real Field
3275 Codomain: Euclidean Jordan algebra of dimension 2 over Algebraic
3277 sage: J.cartesian_projection(1)
3278 Linear operator between finite-dimensional Euclidean Jordan
3279 algebras represented by the matrix:
3283 Domain: Euclidean Jordan algebra of dimension 2 over Algebraic
3284 Real Field (+) Euclidean Jordan algebra of dimension 3 over
3285 Algebraic Real Field
3286 Codomain: Euclidean Jordan algebra of dimension 3 over Algebraic
3289 The projections work the way you'd expect on the vector
3290 representation of an element::
3292 sage: J1 = JordanSpinEJA(2)
3293 sage: J2 = ComplexHermitianEJA(2)
3294 sage: J = cartesian_product([J1,J2])
3295 sage: pi_left = J.cartesian_projection(0)
3296 sage: pi_right = J.cartesian_projection(1)
3297 sage: pi_left(J.one()).to_vector()
3299 sage: pi_right(J.one()).to_vector()
3301 sage: J.one().to_vector()
3306 The answer never changes::
3308 sage: set_random_seed()
3309 sage: J1 = random_eja()
3310 sage: J2 = random_eja()
3311 sage: J = cartesian_product([J1,J2])
3312 sage: P0 = J.cartesian_projection(0)
3313 sage: P1 = J.cartesian_projection(0)
3318 offset
= sum( self
.cartesian_factor(k
).dimension()
3320 Ji
= self
.cartesian_factor(i
)
3321 Pi
= self
._module
_morphism
(lambda j
: Ji
.monomial(j
- offset
),
3324 return FiniteDimensionalEJAOperator(self
,Ji
,Pi
.matrix())
3327 def cartesian_embedding(self
, i
):
3331 sage: from mjo.eja.eja_algebra import (random_eja,
3332 ....: JordanSpinEJA,
3334 ....: RealSymmetricEJA)
3338 The embedding morphisms are Euclidean Jordan algebra
3341 sage: J1 = HadamardEJA(2)
3342 sage: J2 = RealSymmetricEJA(2)
3343 sage: J = cartesian_product([J1,J2])
3344 sage: J.cartesian_embedding(0)
3345 Linear operator between finite-dimensional Euclidean Jordan
3346 algebras represented by the matrix:
3352 Domain: Euclidean Jordan algebra of dimension 2 over
3353 Algebraic Real Field
3354 Codomain: Euclidean Jordan algebra of dimension 2 over
3355 Algebraic Real Field (+) Euclidean Jordan algebra of
3356 dimension 3 over Algebraic Real Field
3357 sage: J.cartesian_embedding(1)
3358 Linear operator between finite-dimensional Euclidean Jordan
3359 algebras represented by the matrix:
3365 Domain: Euclidean Jordan algebra of dimension 3 over
3366 Algebraic Real Field
3367 Codomain: Euclidean Jordan algebra of dimension 2 over
3368 Algebraic Real Field (+) Euclidean Jordan algebra of
3369 dimension 3 over Algebraic Real Field
3371 The embeddings work the way you'd expect on the vector
3372 representation of an element::
3374 sage: J1 = JordanSpinEJA(3)
3375 sage: J2 = RealSymmetricEJA(2)
3376 sage: J = cartesian_product([J1,J2])
3377 sage: iota_left = J.cartesian_embedding(0)
3378 sage: iota_right = J.cartesian_embedding(1)
3379 sage: iota_left(J1.zero()) == J.zero()
3381 sage: iota_right(J2.zero()) == J.zero()
3383 sage: J1.one().to_vector()
3385 sage: iota_left(J1.one()).to_vector()
3387 sage: J2.one().to_vector()
3389 sage: iota_right(J2.one()).to_vector()
3391 sage: J.one().to_vector()
3396 The answer never changes::
3398 sage: set_random_seed()
3399 sage: J1 = random_eja()
3400 sage: J2 = random_eja()
3401 sage: J = cartesian_product([J1,J2])
3402 sage: E0 = J.cartesian_embedding(0)
3403 sage: E1 = J.cartesian_embedding(0)
3407 Composing a projection with the corresponding inclusion should
3408 produce the identity map, and mismatching them should produce
3411 sage: set_random_seed()
3412 sage: J1 = random_eja()
3413 sage: J2 = random_eja()
3414 sage: J = cartesian_product([J1,J2])
3415 sage: iota_left = J.cartesian_embedding(0)
3416 sage: iota_right = J.cartesian_embedding(1)
3417 sage: pi_left = J.cartesian_projection(0)
3418 sage: pi_right = J.cartesian_projection(1)
3419 sage: pi_left*iota_left == J1.one().operator()
3421 sage: pi_right*iota_right == J2.one().operator()
3423 sage: (pi_left*iota_right).is_zero()
3425 sage: (pi_right*iota_left).is_zero()
3429 offset
= sum( self
.cartesian_factor(k
).dimension()
3431 Ji
= self
.cartesian_factor(i
)
3432 Ei
= Ji
._module
_morphism
(lambda j
: self
.monomial(j
+ offset
),
3434 return FiniteDimensionalEJAOperator(Ji
,self
,Ei
.matrix())
3438 FiniteDimensionalEJA
.CartesianProduct
= CartesianProductEJA
3440 class RationalBasisCartesianProductEJA(CartesianProductEJA
,
3443 A separate class for products of algebras for which we know a
3448 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
3449 ....: RealSymmetricEJA)
3453 This gives us fast characteristic polynomial computations in
3454 product algebras, too::
3457 sage: J1 = JordanSpinEJA(2)
3458 sage: J2 = RealSymmetricEJA(3)
3459 sage: J = cartesian_product([J1,J2])
3460 sage: J.characteristic_polynomial_of().degree()
3466 def __init__(self
, algebras
, **kwargs
):
3467 CartesianProductEJA
.__init
__(self
, algebras
, **kwargs
)
3469 self
._rational
_algebra
= None
3470 if self
.vector_space().base_field() is not QQ
:
3471 self
._rational
_algebra
= cartesian_product([
3472 r
._rational
_algebra
for r
in algebras
3476 RationalBasisEJA
.CartesianProduct
= RationalBasisCartesianProductEJA
3478 def random_eja(*args
, **kwargs
):
3479 J1
= ConcreteEJA
.random_instance(*args
, **kwargs
)
3481 # This might make Cartesian products appear roughly as often as
3482 # any other ConcreteEJA.
3483 if ZZ
.random_element(len(ConcreteEJA
.__subclasses
__()) + 1) == 0:
3484 # Use random_eja() again so we can get more than two factors.
3485 J2
= random_eja(*args
, **kwargs
)
3486 J
= cartesian_product([J1
,J2
])