2 Representations and constructions for Euclidean Jordan algebras.
4 A Euclidean Jordan algebra is a Jordan algebra that has some
7 1. It is finite-dimensional.
8 2. Its scalar field is the real numbers.
9 3a. An inner product is defined on it, and...
10 3b. That inner product is compatible with the Jordan product
11 in the sense that `<x*y,z> = <y,x*z>` for all elements
12 `x,y,z` in the algebra.
14 Every Euclidean Jordan algebra is formally-real: for any two elements
15 `x` and `y` in the algebra, `x^{2} + y^{2} = 0` implies that `x = y =
16 0`. Conversely, every finite-dimensional formally-real Jordan algebra
17 can be made into a Euclidean Jordan algebra with an appropriate choice
20 Formally-real Jordan algebras were originally studied as a framework
21 for quantum mechanics. Today, Euclidean Jordan algebras are crucial in
22 symmetric cone optimization, since every symmetric cone arises as the
23 cone of squares in some Euclidean Jordan algebra.
25 It is known that every Euclidean Jordan algebra decomposes into an
26 orthogonal direct sum (essentially, a Cartesian product) of simple
27 algebras, and that moreover, up to Jordan-algebra isomorphism, there
28 are only five families of simple algebras. We provide constructions
29 for these simple algebras:
31 * :class:`BilinearFormEJA`
32 * :class:`RealSymmetricEJA`
33 * :class:`ComplexHermitianEJA`
34 * :class:`QuaternionHermitianEJA`
36 Missing from this list is the algebra of three-by-three octononion
37 Hermitian matrices, as there is (as of yet) no implementation of the
38 octonions in SageMath. In addition to these, we provide two other
39 example constructions,
41 * :class:`HadamardEJA`
44 The Jordan spin algebra is a bilinear form algebra where the bilinear
45 form is the identity. The Hadamard EJA is simply a Cartesian product
46 of one-dimensional spin algebras. And last but not least, the trivial
47 EJA is exactly what you think. Cartesian products of these are also
48 supported using the usual ``cartesian_product()`` function; as a
49 result, we support (up to isomorphism) all Euclidean Jordan algebras
50 that don't involve octonions.
54 sage: from mjo.eja.eja_algebra import random_eja
59 Euclidean Jordan algebra of dimension...
62 from itertools
import repeat
64 from sage
.algebras
.quatalg
.quaternion_algebra
import QuaternionAlgebra
65 from sage
.categories
.magmatic_algebras
import MagmaticAlgebras
66 from sage
.categories
.sets_cat
import cartesian_product
67 from sage
.combinat
.free_module
import (CombinatorialFreeModule
,
68 CombinatorialFreeModule_CartesianProduct
)
69 from sage
.matrix
.constructor
import matrix
70 from sage
.matrix
.matrix_space
import MatrixSpace
71 from sage
.misc
.cachefunc
import cached_method
72 from sage
.misc
.table
import table
73 from sage
.modules
.free_module
import FreeModule
, VectorSpace
74 from sage
.rings
.all
import (ZZ
, QQ
, AA
, QQbar
, RR
, RLF
, CLF
,
77 from mjo
.eja
.eja_element
import FiniteDimensionalEJAElement
78 from mjo
.eja
.eja_operator
import FiniteDimensionalEJAOperator
79 from mjo
.eja
.eja_utils
import _all2list
, _mat2vec
81 class FiniteDimensionalEJA(CombinatorialFreeModule
):
83 A finite-dimensional Euclidean Jordan algebra.
87 - ``basis`` -- a tuple; a tuple of basis elements in "matrix
88 form," which must be the same form as the arguments to
89 ``jordan_product`` and ``inner_product``. In reality, "matrix
90 form" can be either vectors, matrices, or a Cartesian product
91 (ordered tuple) of vectors or matrices. All of these would
92 ideally be vector spaces in sage with no special-casing
93 needed; but in reality we turn vectors into column-matrices
94 and Cartesian products `(a,b)` into column matrices
95 `(a,b)^{T}` after converting `a` and `b` themselves.
97 - ``jordan_product`` -- a function; afunction of two ``basis``
98 elements (in matrix form) that returns their jordan product,
99 also in matrix form; this will be applied to ``basis`` to
100 compute a multiplication table for the algebra.
102 - ``inner_product`` -- a function; a function of two ``basis``
103 elements (in matrix form) that returns their inner
104 product. This will be applied to ``basis`` to compute an
105 inner-product table (basically a matrix) for this algebra.
107 - ``field`` -- a subfield of the reals (default: ``AA``); the scalar
108 field for the algebra.
110 - ``orthonormalize`` -- boolean (default: ``True``); whether or
111 not to orthonormalize the basis. Doing so is expensive and
112 generally rules out using the rationals as your ``field``, but
113 is required for spectral decompositions.
117 sage: from mjo.eja.eja_algebra import random_eja
121 We should compute that an element subalgebra is associative even
122 if we circumvent the element method::
124 sage: set_random_seed()
125 sage: J = random_eja(field=QQ,orthonormalize=False)
126 sage: x = J.random_element()
127 sage: A = x.subalgebra_generated_by(orthonormalize=False)
128 sage: basis = tuple(b.superalgebra_element() for b in A.basis())
129 sage: J.subalgebra(basis, orthonormalize=False).is_associative()
133 Element
= FiniteDimensionalEJAElement
142 cartesian_product
=False,
147 # Keep track of whether or not the matrix basis consists of
148 # tuples, since we need special cases for them damned near
149 # everywhere. This is INDEPENDENT of whether or not the
150 # algebra is a cartesian product, since a subalgebra of a
151 # cartesian product will have a basis of tuples, but will not
152 # in general itself be a cartesian product algebra.
153 self
._matrix
_basis
_is
_cartesian
= False
156 if hasattr(basis
[0], 'cartesian_factors'):
157 self
._matrix
_basis
_is
_cartesian
= True
160 if not field
.is_subring(RR
):
161 # Note: this does return true for the real algebraic
162 # field, the rationals, and any quadratic field where
163 # we've specified a real embedding.
164 raise ValueError("scalar field is not real")
166 # If the basis given to us wasn't over the field that it's
167 # supposed to be over, fix that. Or, you know, crash.
168 if not cartesian_product
:
169 # The field for a cartesian product algebra comes from one
170 # of its factors and is the same for all factors, so
171 # there's no need to "reapply" it on product algebras.
172 if self
._matrix
_basis
_is
_cartesian
:
173 # OK since if n == 0, the basis does not consist of tuples.
174 P
= basis
[0].parent()
175 basis
= tuple( P(tuple(b_i
.change_ring(field
) for b_i
in b
))
178 basis
= tuple( b
.change_ring(field
) for b
in basis
)
182 # Check commutativity of the Jordan and inner-products.
183 # This has to be done before we build the multiplication
184 # and inner-product tables/matrices, because we take
185 # advantage of symmetry in the process.
186 if not all( jordan_product(bi
,bj
) == jordan_product(bj
,bi
)
189 raise ValueError("Jordan product is not commutative")
191 if not all( inner_product(bi
,bj
) == inner_product(bj
,bi
)
194 raise ValueError("inner-product is not commutative")
197 category
= MagmaticAlgebras(field
).FiniteDimensional()
198 category
= category
.WithBasis().Unital().Commutative()
200 if associative
is None:
201 # We should figure it out. As with check_axioms, we have to do
202 # this without the help of the _jordan_product_is_associative()
203 # method because we need to know the category before we
204 # initialize the algebra.
205 associative
= all( jordan_product(jordan_product(bi
,bj
),bk
)
207 jordan_product(bi
,jordan_product(bj
,bk
))
213 # Element subalgebras can take advantage of this.
214 category
= category
.Associative()
215 if cartesian_product
:
216 category
= category
.CartesianProducts()
218 # Call the superclass constructor so that we can use its from_vector()
219 # method to build our multiplication table.
220 CombinatorialFreeModule
.__init
__(self
,
227 # Now comes all of the hard work. We'll be constructing an
228 # ambient vector space V that our (vectorized) basis lives in,
229 # as well as a subspace W of V spanned by those (vectorized)
230 # basis elements. The W-coordinates are the coefficients that
231 # we see in things like x = 1*e1 + 2*e2.
236 degree
= len(_all2list(basis
[0]))
238 # Build an ambient space that fits our matrix basis when
239 # written out as "long vectors."
240 V
= VectorSpace(field
, degree
)
242 # The matrix that will hole the orthonormal -> unorthonormal
243 # coordinate transformation.
244 self
._deortho
_matrix
= None
247 # Save a copy of the un-orthonormalized basis for later.
248 # Convert it to ambient V (vector) coordinates while we're
249 # at it, because we'd have to do it later anyway.
250 deortho_vector_basis
= tuple( V(_all2list(b
)) for b
in basis
)
252 from mjo
.eja
.eja_utils
import gram_schmidt
253 basis
= tuple(gram_schmidt(basis
, inner_product
))
255 # Save the (possibly orthonormalized) matrix basis for
257 self
._matrix
_basis
= basis
259 # Now create the vector space for the algebra, which will have
260 # its own set of non-ambient coordinates (in terms of the
262 vector_basis
= tuple( V(_all2list(b
)) for b
in basis
)
263 W
= V
.span_of_basis( vector_basis
, check
=check_axioms
)
266 # Now "W" is the vector space of our algebra coordinates. The
267 # variables "X1", "X2",... refer to the entries of vectors in
268 # W. Thus to convert back and forth between the orthonormal
269 # coordinates and the given ones, we need to stick the original
271 U
= V
.span_of_basis( deortho_vector_basis
, check
=check_axioms
)
272 self
._deortho
_matrix
= matrix( U
.coordinate_vector(q
)
273 for q
in vector_basis
)
276 # Now we actually compute the multiplication and inner-product
277 # tables/matrices using the possibly-orthonormalized basis.
278 self
._inner
_product
_matrix
= matrix
.identity(field
, n
)
279 self
._multiplication
_table
= [ [0 for j
in range(i
+1)]
282 # Note: the Jordan and inner-products are defined in terms
283 # of the ambient basis. It's important that their arguments
284 # are in ambient coordinates as well.
287 # ortho basis w.r.t. ambient coords
291 # The jordan product returns a matrixy answer, so we
292 # have to convert it to the algebra coordinates.
293 elt
= jordan_product(q_i
, q_j
)
294 elt
= W
.coordinate_vector(V(_all2list(elt
)))
295 self
._multiplication
_table
[i
][j
] = self
.from_vector(elt
)
297 if not orthonormalize
:
298 # If we're orthonormalizing the basis with respect
299 # to an inner-product, then the inner-product
300 # matrix with respect to the resulting basis is
301 # just going to be the identity.
302 ip
= inner_product(q_i
, q_j
)
303 self
._inner
_product
_matrix
[i
,j
] = ip
304 self
._inner
_product
_matrix
[j
,i
] = ip
306 self
._inner
_product
_matrix
._cache
= {'hermitian': True}
307 self
._inner
_product
_matrix
.set_immutable()
310 if not self
._is
_jordanian
():
311 raise ValueError("Jordan identity does not hold")
312 if not self
._inner
_product
_is
_associative
():
313 raise ValueError("inner product is not associative")
316 def _coerce_map_from_base_ring(self
):
318 Disable the map from the base ring into the algebra.
320 Performing a nonsense conversion like this automatically
321 is counterpedagogical. The fallback is to try the usual
322 element constructor, which should also fail.
326 sage: from mjo.eja.eja_algebra import random_eja
330 sage: set_random_seed()
331 sage: J = random_eja()
333 Traceback (most recent call last):
335 ValueError: not an element of this algebra
341 def product_on_basis(self
, i
, j
):
343 Returns the Jordan product of the `i` and `j`th basis elements.
345 This completely defines the Jordan product on the algebra, and
346 is used direclty by our superclass machinery to implement
351 sage: from mjo.eja.eja_algebra import random_eja
355 sage: set_random_seed()
356 sage: J = random_eja()
357 sage: n = J.dimension()
360 sage: ei_ej = J.zero()*J.zero()
362 ....: i = ZZ.random_element(n)
363 ....: j = ZZ.random_element(n)
364 ....: ei = J.gens()[i]
365 ....: ej = J.gens()[j]
366 ....: ei_ej = J.product_on_basis(i,j)
371 # We only stored the lower-triangular portion of the
372 # multiplication table.
374 return self
._multiplication
_table
[i
][j
]
376 return self
._multiplication
_table
[j
][i
]
378 def inner_product(self
, x
, y
):
380 The inner product associated with this Euclidean Jordan algebra.
382 Defaults to the trace inner product, but can be overridden by
383 subclasses if they are sure that the necessary properties are
388 sage: from mjo.eja.eja_algebra import (random_eja,
390 ....: BilinearFormEJA)
394 Our inner product is "associative," which means the following for
395 a symmetric bilinear form::
397 sage: set_random_seed()
398 sage: J = random_eja()
399 sage: x,y,z = J.random_elements(3)
400 sage: (x*y).inner_product(z) == y.inner_product(x*z)
405 Ensure that this is the usual inner product for the algebras
408 sage: set_random_seed()
409 sage: J = HadamardEJA.random_instance()
410 sage: x,y = J.random_elements(2)
411 sage: actual = x.inner_product(y)
412 sage: expected = x.to_vector().inner_product(y.to_vector())
413 sage: actual == expected
416 Ensure that this is one-half of the trace inner-product in a
417 BilinearFormEJA that isn't just the reals (when ``n`` isn't
418 one). This is in Faraut and Koranyi, and also my "On the
421 sage: set_random_seed()
422 sage: J = BilinearFormEJA.random_instance()
423 sage: n = J.dimension()
424 sage: x = J.random_element()
425 sage: y = J.random_element()
426 sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
430 B
= self
._inner
_product
_matrix
431 return (B
*x
.to_vector()).inner_product(y
.to_vector())
434 def is_associative(self
):
436 Return whether or not this algebra's Jordan product is associative.
440 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
444 sage: J = ComplexHermitianEJA(3, field=QQ, orthonormalize=False)
445 sage: J.is_associative()
447 sage: x = sum(J.gens())
448 sage: A = x.subalgebra_generated_by(orthonormalize=False)
449 sage: A.is_associative()
453 return "Associative" in self
.category().axioms()
455 def _is_commutative(self
):
457 Whether or not this algebra's multiplication table is commutative.
459 This method should of course always return ``True``, unless
460 this algebra was constructed with ``check_axioms=False`` and
461 passed an invalid multiplication table.
463 return all( x
*y
== y
*x
for x
in self
.gens() for y
in self
.gens() )
465 def _is_jordanian(self
):
467 Whether or not this algebra's multiplication table respects the
468 Jordan identity `(x^{2})(xy) = x(x^{2}y)`.
470 We only check one arrangement of `x` and `y`, so for a
471 ``True`` result to be truly true, you should also check
472 :meth:`_is_commutative`. This method should of course always
473 return ``True``, unless this algebra was constructed with
474 ``check_axioms=False`` and passed an invalid multiplication table.
476 return all( (self
.gens()[i
]**2)*(self
.gens()[i
]*self
.gens()[j
])
478 (self
.gens()[i
])*((self
.gens()[i
]**2)*self
.gens()[j
])
479 for i
in range(self
.dimension())
480 for j
in range(self
.dimension()) )
482 def _jordan_product_is_associative(self
):
484 Return whether or not this algebra's Jordan product is
485 associative; that is, whether or not `x*(y*z) = (x*y)*z`
488 This method should agree with :meth:`is_associative` unless
489 you lied about the value of the ``associative`` parameter
490 when you constructed the algebra.
494 sage: from mjo.eja.eja_algebra import (random_eja,
495 ....: RealSymmetricEJA,
496 ....: ComplexHermitianEJA,
497 ....: QuaternionHermitianEJA)
501 sage: J = RealSymmetricEJA(4, orthonormalize=False)
502 sage: J._jordan_product_is_associative()
504 sage: x = sum(J.gens())
505 sage: A = x.subalgebra_generated_by()
506 sage: A._jordan_product_is_associative()
511 sage: J = ComplexHermitianEJA(2,field=QQ,orthonormalize=False)
512 sage: J._jordan_product_is_associative()
514 sage: x = sum(J.gens())
515 sage: A = x.subalgebra_generated_by(orthonormalize=False)
516 sage: A._jordan_product_is_associative()
521 sage: J = QuaternionHermitianEJA(2)
522 sage: J._jordan_product_is_associative()
524 sage: x = sum(J.gens())
525 sage: A = x.subalgebra_generated_by()
526 sage: A._jordan_product_is_associative()
531 The values we've presupplied to the constructors agree with
534 sage: set_random_seed()
535 sage: J = random_eja()
536 sage: J.is_associative() == J._jordan_product_is_associative()
542 # Used to check whether or not something is zero.
545 # I don't know of any examples that make this magnitude
546 # necessary because I don't know how to make an
547 # associative algebra when the element subalgebra
548 # construction is unreliable (as it is over RDF; we can't
549 # find the degree of an element because we can't compute
550 # the rank of a matrix). But even multiplication of floats
551 # is non-associative, so *some* epsilon is needed... let's
552 # just take the one from _inner_product_is_associative?
555 for i
in range(self
.dimension()):
556 for j
in range(self
.dimension()):
557 for k
in range(self
.dimension()):
561 diff
= (x
*y
)*z
- x
*(y
*z
)
563 if diff
.norm() > epsilon
:
568 def _inner_product_is_associative(self
):
570 Return whether or not this algebra's inner product `B` is
571 associative; that is, whether or not `B(xy,z) = B(x,yz)`.
573 This method should of course always return ``True``, unless
574 this algebra was constructed with ``check_axioms=False`` and
575 passed an invalid Jordan or inner-product.
579 # Used to check whether or not something is zero.
582 # This choice is sufficient to allow the construction of
583 # QuaternionHermitianEJA(2, field=RDF) with check_axioms=True.
586 for i
in range(self
.dimension()):
587 for j
in range(self
.dimension()):
588 for k
in range(self
.dimension()):
592 diff
= (x
*y
).inner_product(z
) - x
.inner_product(y
*z
)
594 if diff
.abs() > epsilon
:
599 def _element_constructor_(self
, elt
):
601 Construct an element of this algebra from its vector or matrix
604 This gets called only after the parent element _call_ method
605 fails to find a coercion for the argument.
609 sage: from mjo.eja.eja_algebra import (random_eja,
612 ....: RealSymmetricEJA)
616 The identity in `S^n` is converted to the identity in the EJA::
618 sage: J = RealSymmetricEJA(3)
619 sage: I = matrix.identity(QQ,3)
620 sage: J(I) == J.one()
623 This skew-symmetric matrix can't be represented in the EJA::
625 sage: J = RealSymmetricEJA(3)
626 sage: A = matrix(QQ,3, lambda i,j: i-j)
628 Traceback (most recent call last):
630 ValueError: not an element of this algebra
632 Tuples work as well, provided that the matrix basis for the
633 algebra consists of them::
635 sage: J1 = HadamardEJA(3)
636 sage: J2 = RealSymmetricEJA(2)
637 sage: J = cartesian_product([J1,J2])
638 sage: J( (J1.matrix_basis()[1], J2.matrix_basis()[2]) )
643 Ensure that we can convert any element back and forth
644 faithfully between its matrix and algebra representations::
646 sage: set_random_seed()
647 sage: J = random_eja()
648 sage: x = J.random_element()
649 sage: J(x.to_matrix()) == x
652 We cannot coerce elements between algebras just because their
653 matrix representations are compatible::
655 sage: J1 = HadamardEJA(3)
656 sage: J2 = JordanSpinEJA(3)
658 Traceback (most recent call last):
660 ValueError: not an element of this algebra
662 Traceback (most recent call last):
664 ValueError: not an element of this algebra
666 msg
= "not an element of this algebra"
667 if elt
in self
.base_ring():
668 # Ensure that no base ring -> algebra coercion is performed
669 # by this method. There's some stupidity in sage that would
670 # otherwise propagate to this method; for example, sage thinks
671 # that the integer 3 belongs to the space of 2-by-2 matrices.
672 raise ValueError(msg
)
675 # Try to convert a vector into a column-matrix...
677 except (AttributeError, TypeError):
678 # and ignore failure, because we weren't really expecting
679 # a vector as an argument anyway.
682 if elt
not in self
.matrix_space():
683 raise ValueError(msg
)
685 # Thanks for nothing! Matrix spaces aren't vector spaces in
686 # Sage, so we have to figure out its matrix-basis coordinates
687 # ourselves. We use the basis space's ring instead of the
688 # element's ring because the basis space might be an algebraic
689 # closure whereas the base ring of the 3-by-3 identity matrix
690 # could be QQ instead of QQbar.
692 # And, we also have to handle Cartesian product bases (when
693 # the matrix basis consists of tuples) here. The "good news"
694 # is that we're already converting everything to long vectors,
695 # and that strategy works for tuples as well.
697 # We pass check=False because the matrix basis is "guaranteed"
698 # to be linearly independent... right? Ha ha.
700 V
= VectorSpace(self
.base_ring(), len(elt
))
701 W
= V
.span_of_basis( (V(_all2list(s
)) for s
in self
.matrix_basis()),
705 coords
= W
.coordinate_vector(V(elt
))
706 except ArithmeticError: # vector is not in free module
707 raise ValueError(msg
)
709 return self
.from_vector(coords
)
713 Return a string representation of ``self``.
717 sage: from mjo.eja.eja_algebra import JordanSpinEJA
721 Ensure that it says what we think it says::
723 sage: JordanSpinEJA(2, field=AA)
724 Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
725 sage: JordanSpinEJA(3, field=RDF)
726 Euclidean Jordan algebra of dimension 3 over Real Double Field
729 fmt
= "Euclidean Jordan algebra of dimension {} over {}"
730 return fmt
.format(self
.dimension(), self
.base_ring())
734 def characteristic_polynomial_of(self
):
736 Return the algebra's "characteristic polynomial of" function,
737 which is itself a multivariate polynomial that, when evaluated
738 at the coordinates of some algebra element, returns that
739 element's characteristic polynomial.
741 The resulting polynomial has `n+1` variables, where `n` is the
742 dimension of this algebra. The first `n` variables correspond to
743 the coordinates of an algebra element: when evaluated at the
744 coordinates of an algebra element with respect to a certain
745 basis, the result is a univariate polynomial (in the one
746 remaining variable ``t``), namely the characteristic polynomial
751 sage: from mjo.eja.eja_algebra import JordanSpinEJA, TrivialEJA
755 The characteristic polynomial in the spin algebra is given in
756 Alizadeh, Example 11.11::
758 sage: J = JordanSpinEJA(3)
759 sage: p = J.characteristic_polynomial_of(); p
760 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
761 sage: xvec = J.one().to_vector()
765 By definition, the characteristic polynomial is a monic
766 degree-zero polynomial in a rank-zero algebra. Note that
767 Cayley-Hamilton is indeed satisfied since the polynomial
768 ``1`` evaluates to the identity element of the algebra on
771 sage: J = TrivialEJA()
772 sage: J.characteristic_polynomial_of()
779 # The list of coefficient polynomials a_0, a_1, a_2, ..., a_(r-1).
780 a
= self
._charpoly
_coefficients
()
782 # We go to a bit of trouble here to reorder the
783 # indeterminates, so that it's easier to evaluate the
784 # characteristic polynomial at x's coordinates and get back
785 # something in terms of t, which is what we want.
786 S
= PolynomialRing(self
.base_ring(),'t')
790 S
= PolynomialRing(S
, R
.variable_names())
793 return (t
**r
+ sum( a
[k
]*(t
**k
) for k
in range(r
) ))
795 def coordinate_polynomial_ring(self
):
797 The multivariate polynomial ring in which this algebra's
798 :meth:`characteristic_polynomial_of` lives.
802 sage: from mjo.eja.eja_algebra import (HadamardEJA,
803 ....: RealSymmetricEJA)
807 sage: J = HadamardEJA(2)
808 sage: J.coordinate_polynomial_ring()
809 Multivariate Polynomial Ring in X1, X2...
810 sage: J = RealSymmetricEJA(3,field=QQ,orthonormalize=False)
811 sage: J.coordinate_polynomial_ring()
812 Multivariate Polynomial Ring in X1, X2, X3, X4, X5, X6...
815 var_names
= tuple( "X%d" % z
for z
in range(1, self
.dimension()+1) )
816 return PolynomialRing(self
.base_ring(), var_names
)
818 def inner_product(self
, x
, y
):
820 The inner product associated with this Euclidean Jordan algebra.
822 Defaults to the trace inner product, but can be overridden by
823 subclasses if they are sure that the necessary properties are
828 sage: from mjo.eja.eja_algebra import (random_eja,
830 ....: BilinearFormEJA)
834 Our inner product is "associative," which means the following for
835 a symmetric bilinear form::
837 sage: set_random_seed()
838 sage: J = random_eja()
839 sage: x,y,z = J.random_elements(3)
840 sage: (x*y).inner_product(z) == y.inner_product(x*z)
845 Ensure that this is the usual inner product for the algebras
848 sage: set_random_seed()
849 sage: J = HadamardEJA.random_instance()
850 sage: x,y = J.random_elements(2)
851 sage: actual = x.inner_product(y)
852 sage: expected = x.to_vector().inner_product(y.to_vector())
853 sage: actual == expected
856 Ensure that this is one-half of the trace inner-product in a
857 BilinearFormEJA that isn't just the reals (when ``n`` isn't
858 one). This is in Faraut and Koranyi, and also my "On the
861 sage: set_random_seed()
862 sage: J = BilinearFormEJA.random_instance()
863 sage: n = J.dimension()
864 sage: x = J.random_element()
865 sage: y = J.random_element()
866 sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
869 B
= self
._inner
_product
_matrix
870 return (B
*x
.to_vector()).inner_product(y
.to_vector())
873 def is_trivial(self
):
875 Return whether or not this algebra is trivial.
877 A trivial algebra contains only the zero element.
881 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
886 sage: J = ComplexHermitianEJA(3)
892 sage: J = TrivialEJA()
897 return self
.dimension() == 0
900 def multiplication_table(self
):
902 Return a visual representation of this algebra's multiplication
903 table (on basis elements).
907 sage: from mjo.eja.eja_algebra import JordanSpinEJA
911 sage: J = JordanSpinEJA(4)
912 sage: J.multiplication_table()
913 +----++----+----+----+----+
914 | * || e0 | e1 | e2 | e3 |
915 +====++====+====+====+====+
916 | e0 || e0 | e1 | e2 | e3 |
917 +----++----+----+----+----+
918 | e1 || e1 | e0 | 0 | 0 |
919 +----++----+----+----+----+
920 | e2 || e2 | 0 | e0 | 0 |
921 +----++----+----+----+----+
922 | e3 || e3 | 0 | 0 | e0 |
923 +----++----+----+----+----+
927 # Prepend the header row.
928 M
= [["*"] + list(self
.gens())]
930 # And to each subsequent row, prepend an entry that belongs to
931 # the left-side "header column."
932 M
+= [ [self
.gens()[i
]] + [ self
.gens()[i
]*self
.gens()[j
]
936 return table(M
, header_row
=True, header_column
=True, frame
=True)
939 def matrix_basis(self
):
941 Return an (often more natural) representation of this algebras
942 basis as an ordered tuple of matrices.
944 Every finite-dimensional Euclidean Jordan Algebra is a, up to
945 Jordan isomorphism, a direct sum of five simple
946 algebras---four of which comprise Hermitian matrices. And the
947 last type of algebra can of course be thought of as `n`-by-`1`
948 column matrices (ambiguusly called column vectors) to avoid
949 special cases. As a result, matrices (and column vectors) are
950 a natural representation format for Euclidean Jordan algebra
953 But, when we construct an algebra from a basis of matrices,
954 those matrix representations are lost in favor of coordinate
955 vectors *with respect to* that basis. We could eventually
956 convert back if we tried hard enough, but having the original
957 representations handy is valuable enough that we simply store
958 them and return them from this method.
960 Why implement this for non-matrix algebras? Avoiding special
961 cases for the :class:`BilinearFormEJA` pays with simplicity in
962 its own right. But mainly, we would like to be able to assume
963 that elements of a :class:`CartesianProductEJA` can be displayed
964 nicely, without having to have special classes for direct sums
965 one of whose components was a matrix algebra.
969 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
970 ....: RealSymmetricEJA)
974 sage: J = RealSymmetricEJA(2)
976 Finite family {0: e0, 1: e1, 2: e2}
977 sage: J.matrix_basis()
979 [1 0] [ 0 0.7071067811865475?] [0 0]
980 [0 0], [0.7071067811865475? 0], [0 1]
985 sage: J = JordanSpinEJA(2)
987 Finite family {0: e0, 1: e1}
988 sage: J.matrix_basis()
994 return self
._matrix
_basis
997 def matrix_space(self
):
999 Return the matrix space in which this algebra's elements live, if
1000 we think of them as matrices (including column vectors of the
1003 "By default" this will be an `n`-by-`1` column-matrix space,
1004 except when the algebra is trivial. There it's `n`-by-`n`
1005 (where `n` is zero), to ensure that two elements of the matrix
1006 space (empty matrices) can be multiplied. For algebras of
1007 matrices, this returns the space in which their
1008 real embeddings live.
1012 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
1013 ....: JordanSpinEJA,
1014 ....: QuaternionHermitianEJA,
1019 By default, the matrix representation is just a column-matrix
1020 equivalent to the vector representation::
1022 sage: J = JordanSpinEJA(3)
1023 sage: J.matrix_space()
1024 Full MatrixSpace of 3 by 1 dense matrices over Algebraic
1027 The matrix representation in the trivial algebra is
1028 zero-by-zero instead of the usual `n`-by-one::
1030 sage: J = TrivialEJA()
1031 sage: J.matrix_space()
1032 Full MatrixSpace of 0 by 0 dense matrices over Algebraic
1035 The matrix space for complex/quaternion Hermitian matrix EJA
1036 is the space in which their real-embeddings live, not the
1037 original complex/quaternion matrix space::
1039 sage: J = ComplexHermitianEJA(2,field=QQ,orthonormalize=False)
1040 sage: J.matrix_space()
1041 Full MatrixSpace of 4 by 4 dense matrices over Rational Field
1042 sage: J = QuaternionHermitianEJA(1,field=QQ,orthonormalize=False)
1043 sage: J.matrix_space()
1044 Full MatrixSpace of 4 by 4 dense matrices over Rational Field
1047 if self
.is_trivial():
1048 return MatrixSpace(self
.base_ring(), 0)
1050 return self
.matrix_basis()[0].parent()
1056 Return the unit element of this algebra.
1060 sage: from mjo.eja.eja_algebra import (HadamardEJA,
1065 We can compute unit element in the Hadamard EJA::
1067 sage: J = HadamardEJA(5)
1069 e0 + e1 + e2 + e3 + e4
1071 The unit element in the Hadamard EJA is inherited in the
1072 subalgebras generated by its elements::
1074 sage: J = HadamardEJA(5)
1076 e0 + e1 + e2 + e3 + e4
1077 sage: x = sum(J.gens())
1078 sage: A = x.subalgebra_generated_by(orthonormalize=False)
1081 sage: A.one().superalgebra_element()
1082 e0 + e1 + e2 + e3 + e4
1086 The identity element acts like the identity, regardless of
1087 whether or not we orthonormalize::
1089 sage: set_random_seed()
1090 sage: J = random_eja()
1091 sage: x = J.random_element()
1092 sage: J.one()*x == x and x*J.one() == x
1094 sage: A = x.subalgebra_generated_by()
1095 sage: y = A.random_element()
1096 sage: A.one()*y == y and y*A.one() == y
1101 sage: set_random_seed()
1102 sage: J = random_eja(field=QQ, orthonormalize=False)
1103 sage: x = J.random_element()
1104 sage: J.one()*x == x and x*J.one() == x
1106 sage: A = x.subalgebra_generated_by(orthonormalize=False)
1107 sage: y = A.random_element()
1108 sage: A.one()*y == y and y*A.one() == y
1111 The matrix of the unit element's operator is the identity,
1112 regardless of the base field and whether or not we
1115 sage: set_random_seed()
1116 sage: J = random_eja()
1117 sage: actual = J.one().operator().matrix()
1118 sage: expected = matrix.identity(J.base_ring(), J.dimension())
1119 sage: actual == expected
1121 sage: x = J.random_element()
1122 sage: A = x.subalgebra_generated_by()
1123 sage: actual = A.one().operator().matrix()
1124 sage: expected = matrix.identity(A.base_ring(), A.dimension())
1125 sage: actual == expected
1130 sage: set_random_seed()
1131 sage: J = random_eja(field=QQ, orthonormalize=False)
1132 sage: actual = J.one().operator().matrix()
1133 sage: expected = matrix.identity(J.base_ring(), J.dimension())
1134 sage: actual == expected
1136 sage: x = J.random_element()
1137 sage: A = x.subalgebra_generated_by(orthonormalize=False)
1138 sage: actual = A.one().operator().matrix()
1139 sage: expected = matrix.identity(A.base_ring(), A.dimension())
1140 sage: actual == expected
1143 Ensure that the cached unit element (often precomputed by
1144 hand) agrees with the computed one::
1146 sage: set_random_seed()
1147 sage: J = random_eja()
1148 sage: cached = J.one()
1149 sage: J.one.clear_cache()
1150 sage: J.one() == cached
1155 sage: set_random_seed()
1156 sage: J = random_eja(field=QQ, orthonormalize=False)
1157 sage: cached = J.one()
1158 sage: J.one.clear_cache()
1159 sage: J.one() == cached
1163 # We can brute-force compute the matrices of the operators
1164 # that correspond to the basis elements of this algebra.
1165 # If some linear combination of those basis elements is the
1166 # algebra identity, then the same linear combination of
1167 # their matrices has to be the identity matrix.
1169 # Of course, matrices aren't vectors in sage, so we have to
1170 # appeal to the "long vectors" isometry.
1171 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
1173 # Now we use basic linear algebra to find the coefficients,
1174 # of the matrices-as-vectors-linear-combination, which should
1175 # work for the original algebra basis too.
1176 A
= matrix(self
.base_ring(), oper_vecs
)
1178 # We used the isometry on the left-hand side already, but we
1179 # still need to do it for the right-hand side. Recall that we
1180 # wanted something that summed to the identity matrix.
1181 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
1183 # Now if there's an identity element in the algebra, this
1184 # should work. We solve on the left to avoid having to
1185 # transpose the matrix "A".
1186 return self
.from_vector(A
.solve_left(b
))
1189 def peirce_decomposition(self
, c
):
1191 The Peirce decomposition of this algebra relative to the
1194 In the future, this can be extended to a complete system of
1195 orthogonal idempotents.
1199 - ``c`` -- an idempotent of this algebra.
1203 A triple (J0, J5, J1) containing two subalgebras and one subspace
1206 - ``J0`` -- the algebra on the eigenspace of ``c.operator()``
1207 corresponding to the eigenvalue zero.
1209 - ``J5`` -- the eigenspace (NOT a subalgebra) of ``c.operator()``
1210 corresponding to the eigenvalue one-half.
1212 - ``J1`` -- the algebra on the eigenspace of ``c.operator()``
1213 corresponding to the eigenvalue one.
1215 These are the only possible eigenspaces for that operator, and this
1216 algebra is a direct sum of them. The spaces ``J0`` and ``J1`` are
1217 orthogonal, and are subalgebras of this algebra with the appropriate
1222 sage: from mjo.eja.eja_algebra import random_eja, RealSymmetricEJA
1226 The canonical example comes from the symmetric matrices, which
1227 decompose into diagonal and off-diagonal parts::
1229 sage: J = RealSymmetricEJA(3)
1230 sage: C = matrix(QQ, [ [1,0,0],
1234 sage: J0,J5,J1 = J.peirce_decomposition(c)
1236 Euclidean Jordan algebra of dimension 1...
1238 Vector space of degree 6 and dimension 2...
1240 Euclidean Jordan algebra of dimension 3...
1241 sage: J0.one().to_matrix()
1245 sage: orig_df = AA.options.display_format
1246 sage: AA.options.display_format = 'radical'
1247 sage: J.from_vector(J5.basis()[0]).to_matrix()
1251 sage: J.from_vector(J5.basis()[1]).to_matrix()
1255 sage: AA.options.display_format = orig_df
1256 sage: J1.one().to_matrix()
1263 Every algebra decomposes trivially with respect to its identity
1266 sage: set_random_seed()
1267 sage: J = random_eja()
1268 sage: J0,J5,J1 = J.peirce_decomposition(J.one())
1269 sage: J0.dimension() == 0 and J5.dimension() == 0
1271 sage: J1.superalgebra() == J and J1.dimension() == J.dimension()
1274 The decomposition is into eigenspaces, and its components are
1275 therefore necessarily orthogonal. Moreover, the identity
1276 elements in the two subalgebras are the projections onto their
1277 respective subspaces of the superalgebra's identity element::
1279 sage: set_random_seed()
1280 sage: J = random_eja()
1281 sage: x = J.random_element()
1282 sage: if not J.is_trivial():
1283 ....: while x.is_nilpotent():
1284 ....: x = J.random_element()
1285 sage: c = x.subalgebra_idempotent()
1286 sage: J0,J5,J1 = J.peirce_decomposition(c)
1288 sage: for (w,y,z) in zip(J0.basis(), J5.basis(), J1.basis()):
1289 ....: w = w.superalgebra_element()
1290 ....: y = J.from_vector(y)
1291 ....: z = z.superalgebra_element()
1292 ....: ipsum += w.inner_product(y).abs()
1293 ....: ipsum += w.inner_product(z).abs()
1294 ....: ipsum += y.inner_product(z).abs()
1297 sage: J1(c) == J1.one()
1299 sage: J0(J.one() - c) == J0.one()
1303 if not c
.is_idempotent():
1304 raise ValueError("element is not idempotent: %s" % c
)
1306 # Default these to what they should be if they turn out to be
1307 # trivial, because eigenspaces_left() won't return eigenvalues
1308 # corresponding to trivial spaces (e.g. it returns only the
1309 # eigenspace corresponding to lambda=1 if you take the
1310 # decomposition relative to the identity element).
1311 trivial
= self
.subalgebra(())
1312 J0
= trivial
# eigenvalue zero
1313 J5
= VectorSpace(self
.base_ring(), 0) # eigenvalue one-half
1314 J1
= trivial
# eigenvalue one
1316 for (eigval
, eigspace
) in c
.operator().matrix().right_eigenspaces():
1317 if eigval
== ~
(self
.base_ring()(2)):
1320 gens
= tuple( self
.from_vector(b
) for b
in eigspace
.basis() )
1321 subalg
= self
.subalgebra(gens
, check_axioms
=False)
1327 raise ValueError("unexpected eigenvalue: %s" % eigval
)
1332 def random_element(self
, thorough
=False):
1334 Return a random element of this algebra.
1336 Our algebra superclass method only returns a linear
1337 combination of at most two basis elements. We instead
1338 want the vector space "random element" method that
1339 returns a more diverse selection.
1343 - ``thorough`` -- (boolean; default False) whether or not we
1344 should generate irrational coefficients for the random
1345 element when our base ring is irrational; this slows the
1346 algebra operations to a crawl, but any truly random method
1350 # For a general base ring... maybe we can trust this to do the
1351 # right thing? Unlikely, but.
1352 V
= self
.vector_space()
1353 v
= V
.random_element()
1355 if self
.base_ring() is AA
:
1356 # The "random element" method of the algebraic reals is
1357 # stupid at the moment, and only returns integers between
1358 # -2 and 2, inclusive:
1360 # https://trac.sagemath.org/ticket/30875
1362 # Instead, we implement our own "random vector" method,
1363 # and then coerce that into the algebra. We use the vector
1364 # space degree here instead of the dimension because a
1365 # subalgebra could (for example) be spanned by only two
1366 # vectors, each with five coordinates. We need to
1367 # generate all five coordinates.
1369 v
*= QQbar
.random_element().real()
1371 v
*= QQ
.random_element()
1373 return self
.from_vector(V
.coordinate_vector(v
))
1375 def random_elements(self
, count
, thorough
=False):
1377 Return ``count`` random elements as a tuple.
1381 - ``thorough`` -- (boolean; default False) whether or not we
1382 should generate irrational coefficients for the random
1383 elements when our base ring is irrational; this slows the
1384 algebra operations to a crawl, but any truly random method
1389 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1393 sage: J = JordanSpinEJA(3)
1394 sage: x,y,z = J.random_elements(3)
1395 sage: all( [ x in J, y in J, z in J ])
1397 sage: len( J.random_elements(10) ) == 10
1401 return tuple( self
.random_element(thorough
)
1402 for idx
in range(count
) )
1406 def _charpoly_coefficients(self
):
1408 The `r` polynomial coefficients of the "characteristic polynomial
1413 sage: from mjo.eja.eja_algebra import random_eja
1417 The theory shows that these are all homogeneous polynomials of
1420 sage: set_random_seed()
1421 sage: J = random_eja()
1422 sage: all(p.is_homogeneous() for p in J._charpoly_coefficients())
1426 n
= self
.dimension()
1427 R
= self
.coordinate_polynomial_ring()
1429 F
= R
.fraction_field()
1432 # From a result in my book, these are the entries of the
1433 # basis representation of L_x.
1434 return sum( vars[k
]*self
.gens()[k
].operator().matrix()[i
,j
]
1437 L_x
= matrix(F
, n
, n
, L_x_i_j
)
1440 if self
.rank
.is_in_cache():
1442 # There's no need to pad the system with redundant
1443 # columns if we *know* they'll be redundant.
1446 # Compute an extra power in case the rank is equal to
1447 # the dimension (otherwise, we would stop at x^(r-1)).
1448 x_powers
= [ (L_x
**k
)*self
.one().to_vector()
1449 for k
in range(n
+1) ]
1450 A
= matrix
.column(F
, x_powers
[:n
])
1451 AE
= A
.extended_echelon_form()
1458 # The theory says that only the first "r" coefficients are
1459 # nonzero, and they actually live in the original polynomial
1460 # ring and not the fraction field. We negate them because in
1461 # the actual characteristic polynomial, they get moved to the
1462 # other side where x^r lives. We don't bother to trim A_rref
1463 # down to a square matrix and solve the resulting system,
1464 # because the upper-left r-by-r portion of A_rref is
1465 # guaranteed to be the identity matrix, so e.g.
1467 # A_rref.solve_right(Y)
1469 # would just be returning Y.
1470 return (-E
*b
)[:r
].change_ring(R
)
1475 Return the rank of this EJA.
1477 This is a cached method because we know the rank a priori for
1478 all of the algebras we can construct. Thus we can avoid the
1479 expensive ``_charpoly_coefficients()`` call unless we truly
1480 need to compute the whole characteristic polynomial.
1484 sage: from mjo.eja.eja_algebra import (HadamardEJA,
1485 ....: JordanSpinEJA,
1486 ....: RealSymmetricEJA,
1487 ....: ComplexHermitianEJA,
1488 ....: QuaternionHermitianEJA,
1493 The rank of the Jordan spin algebra is always two::
1495 sage: JordanSpinEJA(2).rank()
1497 sage: JordanSpinEJA(3).rank()
1499 sage: JordanSpinEJA(4).rank()
1502 The rank of the `n`-by-`n` Hermitian real, complex, or
1503 quaternion matrices is `n`::
1505 sage: RealSymmetricEJA(4).rank()
1507 sage: ComplexHermitianEJA(3).rank()
1509 sage: QuaternionHermitianEJA(2).rank()
1514 Ensure that every EJA that we know how to construct has a
1515 positive integer rank, unless the algebra is trivial in
1516 which case its rank will be zero::
1518 sage: set_random_seed()
1519 sage: J = random_eja()
1523 sage: r > 0 or (r == 0 and J.is_trivial())
1526 Ensure that computing the rank actually works, since the ranks
1527 of all simple algebras are known and will be cached by default::
1529 sage: set_random_seed() # long time
1530 sage: J = random_eja() # long time
1531 sage: cached = J.rank() # long time
1532 sage: J.rank.clear_cache() # long time
1533 sage: J.rank() == cached # long time
1537 return len(self
._charpoly
_coefficients
())
1540 def subalgebra(self
, basis
, **kwargs
):
1542 Create a subalgebra of this algebra from the given basis.
1544 from mjo
.eja
.eja_subalgebra
import FiniteDimensionalEJASubalgebra
1545 return FiniteDimensionalEJASubalgebra(self
, basis
, **kwargs
)
1548 def vector_space(self
):
1550 Return the vector space that underlies this algebra.
1554 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1558 sage: J = RealSymmetricEJA(2)
1559 sage: J.vector_space()
1560 Vector space of dimension 3 over...
1563 return self
.zero().to_vector().parent().ambient_vector_space()
1567 class RationalBasisEJA(FiniteDimensionalEJA
):
1569 New class for algebras whose supplied basis elements have all rational entries.
1573 sage: from mjo.eja.eja_algebra import BilinearFormEJA
1577 The supplied basis is orthonormalized by default::
1579 sage: B = matrix(QQ, [[1, 0, 0], [0, 25, -32], [0, -32, 41]])
1580 sage: J = BilinearFormEJA(B)
1581 sage: J.matrix_basis()
1598 # Abuse the check_field parameter to check that the entries of
1599 # out basis (in ambient coordinates) are in the field QQ.
1600 if not all( all(b_i
in QQ
for b_i
in b
.list()) for b
in basis
):
1601 raise TypeError("basis not rational")
1603 super().__init
__(basis
,
1607 check_field
=check_field
,
1610 self
._rational
_algebra
= None
1612 # There's no point in constructing the extra algebra if this
1613 # one is already rational.
1615 # Note: the same Jordan and inner-products work here,
1616 # because they are necessarily defined with respect to
1617 # ambient coordinates and not any particular basis.
1618 self
._rational
_algebra
= FiniteDimensionalEJA(
1623 associative
=self
.is_associative(),
1624 orthonormalize
=False,
1629 def _charpoly_coefficients(self
):
1633 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
1634 ....: JordanSpinEJA)
1638 The base ring of the resulting polynomial coefficients is what
1639 it should be, and not the rationals (unless the algebra was
1640 already over the rationals)::
1642 sage: J = JordanSpinEJA(3)
1643 sage: J._charpoly_coefficients()
1644 (X1^2 - X2^2 - X3^2, -2*X1)
1645 sage: a0 = J._charpoly_coefficients()[0]
1647 Algebraic Real Field
1648 sage: a0.base_ring()
1649 Algebraic Real Field
1652 if self
._rational
_algebra
is None:
1653 # There's no need to construct *another* algebra over the
1654 # rationals if this one is already over the
1655 # rationals. Likewise, if we never orthonormalized our
1656 # basis, we might as well just use the given one.
1657 return super()._charpoly
_coefficients
()
1659 # Do the computation over the rationals. The answer will be
1660 # the same, because all we've done is a change of basis.
1661 # Then, change back from QQ to our real base ring
1662 a
= ( a_i
.change_ring(self
.base_ring())
1663 for a_i
in self
._rational
_algebra
._charpoly
_coefficients
() )
1665 if self
._deortho
_matrix
is None:
1666 # This can happen if our base ring was, say, AA and we
1667 # chose not to (or didn't need to) orthonormalize. It's
1668 # still faster to do the computations over QQ even if
1669 # the numbers in the boxes stay the same.
1672 # Otherwise, convert the coordinate variables back to the
1673 # deorthonormalized ones.
1674 R
= self
.coordinate_polynomial_ring()
1675 from sage
.modules
.free_module_element
import vector
1676 X
= vector(R
, R
.gens())
1677 BX
= self
._deortho
_matrix
*X
1679 subs_dict
= { X[i]: BX[i] for i in range(len(X)) }
1680 return tuple( a_i
.subs(subs_dict
) for a_i
in a
)
1682 class ConcreteEJA(RationalBasisEJA
):
1684 A class for the Euclidean Jordan algebras that we know by name.
1686 These are the Jordan algebras whose basis, multiplication table,
1687 rank, and so on are known a priori. More to the point, they are
1688 the Euclidean Jordan algebras for which we are able to conjure up
1689 a "random instance."
1693 sage: from mjo.eja.eja_algebra import ConcreteEJA
1697 Our basis is normalized with respect to the algebra's inner
1698 product, unless we specify otherwise::
1700 sage: set_random_seed()
1701 sage: J = ConcreteEJA.random_instance()
1702 sage: all( b.norm() == 1 for b in J.gens() )
1705 Since our basis is orthonormal with respect to the algebra's inner
1706 product, and since we know that this algebra is an EJA, any
1707 left-multiplication operator's matrix will be symmetric because
1708 natural->EJA basis representation is an isometry and within the
1709 EJA the operator is self-adjoint by the Jordan axiom::
1711 sage: set_random_seed()
1712 sage: J = ConcreteEJA.random_instance()
1713 sage: x = J.random_element()
1714 sage: x.operator().is_self_adjoint()
1719 def _max_random_instance_size():
1721 Return an integer "size" that is an upper bound on the size of
1722 this algebra when it is used in a random test
1723 case. Unfortunately, the term "size" is ambiguous -- when
1724 dealing with `R^n` under either the Hadamard or Jordan spin
1725 product, the "size" refers to the dimension `n`. When dealing
1726 with a matrix algebra (real symmetric or complex/quaternion
1727 Hermitian), it refers to the size of the matrix, which is far
1728 less than the dimension of the underlying vector space.
1730 This method must be implemented in each subclass.
1732 raise NotImplementedError
1735 def random_instance(cls
, *args
, **kwargs
):
1737 Return a random instance of this type of algebra.
1739 This method should be implemented in each subclass.
1741 from sage
.misc
.prandom
import choice
1742 eja_class
= choice(cls
.__subclasses
__())
1744 # These all bubble up to the RationalBasisEJA superclass
1745 # constructor, so any (kw)args valid there are also valid
1747 return eja_class
.random_instance(*args
, **kwargs
)
1752 def dimension_over_reals():
1754 The dimension of this matrix's base ring over the reals.
1756 The reals are dimension one over themselves, obviously; that's
1757 just `\mathbb{R}^{1}`. Likewise, the complex numbers `a + bi`
1758 have dimension two. Finally, the quaternions have dimension
1759 four over the reals.
1761 This is used to determine the size of the matrix returned from
1762 :meth:`real_embed`, among other things.
1764 raise NotImplementedError
1767 def real_embed(cls
,M
):
1769 Embed the matrix ``M`` into a space of real matrices.
1771 The matrix ``M`` can have entries in any field at the moment:
1772 the real numbers, complex numbers, or quaternions. And although
1773 they are not a field, we can probably support octonions at some
1774 point, too. This function returns a real matrix that "acts like"
1775 the original with respect to matrix multiplication; i.e.
1777 real_embed(M*N) = real_embed(M)*real_embed(N)
1780 if M
.ncols() != M
.nrows():
1781 raise ValueError("the matrix 'M' must be square")
1786 def real_unembed(cls
,M
):
1788 The inverse of :meth:`real_embed`.
1790 if M
.ncols() != M
.nrows():
1791 raise ValueError("the matrix 'M' must be square")
1792 if not ZZ(M
.nrows()).mod(cls
.dimension_over_reals()).is_zero():
1793 raise ValueError("the matrix 'M' must be a real embedding")
1797 def jordan_product(X
,Y
):
1798 return (X
*Y
+ Y
*X
)/2
1801 def trace_inner_product(cls
,X
,Y
):
1803 Compute the trace inner-product of two real-embeddings.
1807 sage: from mjo.eja.eja_algebra import (RealSymmetricEJA,
1808 ....: ComplexHermitianEJA,
1809 ....: QuaternionHermitianEJA)
1813 This gives the same answer as it would if we computed the trace
1814 from the unembedded (original) matrices::
1816 sage: set_random_seed()
1817 sage: J = RealSymmetricEJA.random_instance()
1818 sage: x,y = J.random_elements(2)
1819 sage: Xe = x.to_matrix()
1820 sage: Ye = y.to_matrix()
1821 sage: X = J.real_unembed(Xe)
1822 sage: Y = J.real_unembed(Ye)
1823 sage: expected = (X*Y).trace()
1824 sage: actual = J.trace_inner_product(Xe,Ye)
1825 sage: actual == expected
1830 sage: set_random_seed()
1831 sage: J = ComplexHermitianEJA.random_instance()
1832 sage: x,y = J.random_elements(2)
1833 sage: Xe = x.to_matrix()
1834 sage: Ye = y.to_matrix()
1835 sage: X = J.real_unembed(Xe)
1836 sage: Y = J.real_unembed(Ye)
1837 sage: expected = (X*Y).trace().real()
1838 sage: actual = J.trace_inner_product(Xe,Ye)
1839 sage: actual == expected
1844 sage: set_random_seed()
1845 sage: J = QuaternionHermitianEJA.random_instance()
1846 sage: x,y = J.random_elements(2)
1847 sage: Xe = x.to_matrix()
1848 sage: Ye = y.to_matrix()
1849 sage: X = J.real_unembed(Xe)
1850 sage: Y = J.real_unembed(Ye)
1851 sage: expected = (X*Y).trace().coefficient_tuple()[0]
1852 sage: actual = J.trace_inner_product(Xe,Ye)
1853 sage: actual == expected
1857 Xu
= cls
.real_unembed(X
)
1858 Yu
= cls
.real_unembed(Y
)
1859 tr
= (Xu
*Yu
).trace()
1862 # Works in QQ, AA, RDF, et cetera.
1864 except AttributeError:
1865 # A quaternion doesn't have a real() method, but does
1866 # have coefficient_tuple() method that returns the
1867 # coefficients of 1, i, j, and k -- in that order.
1868 return tr
.coefficient_tuple()[0]
1871 class RealMatrixEJA(MatrixEJA
):
1873 def dimension_over_reals():
1877 class RealSymmetricEJA(ConcreteEJA
, RealMatrixEJA
):
1879 The rank-n simple EJA consisting of real symmetric n-by-n
1880 matrices, the usual symmetric Jordan product, and the trace inner
1881 product. It has dimension `(n^2 + n)/2` over the reals.
1885 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1889 sage: J = RealSymmetricEJA(2)
1890 sage: e0, e1, e2 = J.gens()
1898 In theory, our "field" can be any subfield of the reals::
1900 sage: RealSymmetricEJA(2, field=RDF, check_axioms=True)
1901 Euclidean Jordan algebra of dimension 3 over Real Double Field
1902 sage: RealSymmetricEJA(2, field=RR, check_axioms=True)
1903 Euclidean Jordan algebra of dimension 3 over Real Field with
1904 53 bits of precision
1908 The dimension of this algebra is `(n^2 + n) / 2`::
1910 sage: set_random_seed()
1911 sage: n_max = RealSymmetricEJA._max_random_instance_size()
1912 sage: n = ZZ.random_element(1, n_max)
1913 sage: J = RealSymmetricEJA(n)
1914 sage: J.dimension() == (n^2 + n)/2
1917 The Jordan multiplication is what we think it is::
1919 sage: set_random_seed()
1920 sage: J = RealSymmetricEJA.random_instance()
1921 sage: x,y = J.random_elements(2)
1922 sage: actual = (x*y).to_matrix()
1923 sage: X = x.to_matrix()
1924 sage: Y = y.to_matrix()
1925 sage: expected = (X*Y + Y*X)/2
1926 sage: actual == expected
1928 sage: J(expected) == x*y
1931 We can change the generator prefix::
1933 sage: RealSymmetricEJA(3, prefix='q').gens()
1934 (q0, q1, q2, q3, q4, q5)
1936 We can construct the (trivial) algebra of rank zero::
1938 sage: RealSymmetricEJA(0)
1939 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1943 def _denormalized_basis(cls
, n
):
1945 Return a basis for the space of real symmetric n-by-n matrices.
1949 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1953 sage: set_random_seed()
1954 sage: n = ZZ.random_element(1,5)
1955 sage: B = RealSymmetricEJA._denormalized_basis(n)
1956 sage: all( M.is_symmetric() for M in B)
1960 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1964 for j
in range(i
+1):
1965 Eij
= matrix(ZZ
, n
, lambda k
,l
: k
==i
and l
==j
)
1969 Sij
= Eij
+ Eij
.transpose()
1975 def _max_random_instance_size():
1976 return 4 # Dimension 10
1979 def random_instance(cls
, **kwargs
):
1981 Return a random instance of this type of algebra.
1983 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1984 return cls(n
, **kwargs
)
1986 def __init__(self
, n
, **kwargs
):
1987 # We know this is a valid EJA, but will double-check
1988 # if the user passes check_axioms=True.
1989 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
1995 super().__init
__(self
._denormalized
_basis
(n
),
1996 self
.jordan_product
,
1997 self
.trace_inner_product
,
1998 associative
=associative
,
2001 # TODO: this could be factored out somehow, but is left here
2002 # because the MatrixEJA is not presently a subclass of the
2003 # FDEJA class that defines rank() and one().
2004 self
.rank
.set_cache(n
)
2005 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
2006 self
.one
.set_cache(self(idV
))
2010 class ComplexMatrixEJA(MatrixEJA
):
2011 # A manual dictionary-cache for the complex_extension() method,
2012 # since apparently @classmethods can't also be @cached_methods.
2013 _complex_extension
= {}
2016 def complex_extension(cls
,field
):
2018 The complex field that we embed/unembed, as an extension
2019 of the given ``field``.
2021 if field
in cls
._complex
_extension
:
2022 return cls
._complex
_extension
[field
]
2024 # Sage doesn't know how to adjoin the complex "i" (the root of
2025 # x^2 + 1) to a field in a general way. Here, we just enumerate
2026 # all of the cases that I have cared to support so far.
2028 # Sage doesn't know how to embed AA into QQbar, i.e. how
2029 # to adjoin sqrt(-1) to AA.
2031 elif not field
.is_exact():
2033 F
= field
.complex_field()
2035 # Works for QQ and... maybe some other fields.
2036 R
= PolynomialRing(field
, 'z')
2038 F
= field
.extension(z
**2 + 1, 'I', embedding
=CLF(-1).sqrt())
2040 cls
._complex
_extension
[field
] = F
2044 def dimension_over_reals():
2048 def real_embed(cls
,M
):
2050 Embed the n-by-n complex matrix ``M`` into the space of real
2051 matrices of size 2n-by-2n via the map the sends each entry `z = a +
2052 bi` to the block matrix ``[[a,b],[-b,a]]``.
2056 sage: from mjo.eja.eja_algebra import ComplexMatrixEJA
2060 sage: F = QuadraticField(-1, 'I')
2061 sage: x1 = F(4 - 2*i)
2062 sage: x2 = F(1 + 2*i)
2065 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
2066 sage: ComplexMatrixEJA.real_embed(M)
2075 Embedding is a homomorphism (isomorphism, in fact)::
2077 sage: set_random_seed()
2078 sage: n = ZZ.random_element(3)
2079 sage: F = QuadraticField(-1, 'I')
2080 sage: X = random_matrix(F, n)
2081 sage: Y = random_matrix(F, n)
2082 sage: Xe = ComplexMatrixEJA.real_embed(X)
2083 sage: Ye = ComplexMatrixEJA.real_embed(Y)
2084 sage: XYe = ComplexMatrixEJA.real_embed(X*Y)
2089 super().real_embed(M
)
2092 # We don't need any adjoined elements...
2093 field
= M
.base_ring().base_ring()
2099 blocks
.append(matrix(field
, 2, [ [ a
, b
],
2102 return matrix
.block(field
, n
, blocks
)
2106 def real_unembed(cls
,M
):
2108 The inverse of _embed_complex_matrix().
2112 sage: from mjo.eja.eja_algebra import ComplexMatrixEJA
2116 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
2117 ....: [-2, 1, -4, 3],
2118 ....: [ 9, 10, 11, 12],
2119 ....: [-10, 9, -12, 11] ])
2120 sage: ComplexMatrixEJA.real_unembed(A)
2122 [ 10*I + 9 12*I + 11]
2126 Unembedding is the inverse of embedding::
2128 sage: set_random_seed()
2129 sage: F = QuadraticField(-1, 'I')
2130 sage: M = random_matrix(F, 3)
2131 sage: Me = ComplexMatrixEJA.real_embed(M)
2132 sage: ComplexMatrixEJA.real_unembed(Me) == M
2136 super().real_unembed(M
)
2138 d
= cls
.dimension_over_reals()
2139 F
= cls
.complex_extension(M
.base_ring())
2142 # Go top-left to bottom-right (reading order), converting every
2143 # 2-by-2 block we see to a single complex element.
2145 for k
in range(n
/d
):
2146 for j
in range(n
/d
):
2147 submat
= M
[d
*k
:d
*k
+d
,d
*j
:d
*j
+d
]
2148 if submat
[0,0] != submat
[1,1]:
2149 raise ValueError('bad on-diagonal submatrix')
2150 if submat
[0,1] != -submat
[1,0]:
2151 raise ValueError('bad off-diagonal submatrix')
2152 z
= submat
[0,0] + submat
[0,1]*i
2155 return matrix(F
, n
/d
, elements
)
2158 class ComplexHermitianEJA(ConcreteEJA
, ComplexMatrixEJA
):
2160 The rank-n simple EJA consisting of complex Hermitian n-by-n
2161 matrices over the real numbers, the usual symmetric Jordan product,
2162 and the real-part-of-trace inner product. It has dimension `n^2` over
2167 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
2171 In theory, our "field" can be any subfield of the reals::
2173 sage: ComplexHermitianEJA(2, field=RDF, check_axioms=True)
2174 Euclidean Jordan algebra of dimension 4 over Real Double Field
2175 sage: ComplexHermitianEJA(2, field=RR, check_axioms=True)
2176 Euclidean Jordan algebra of dimension 4 over Real Field with
2177 53 bits of precision
2181 The dimension of this algebra is `n^2`::
2183 sage: set_random_seed()
2184 sage: n_max = ComplexHermitianEJA._max_random_instance_size()
2185 sage: n = ZZ.random_element(1, n_max)
2186 sage: J = ComplexHermitianEJA(n)
2187 sage: J.dimension() == n^2
2190 The Jordan multiplication is what we think it is::
2192 sage: set_random_seed()
2193 sage: J = ComplexHermitianEJA.random_instance()
2194 sage: x,y = J.random_elements(2)
2195 sage: actual = (x*y).to_matrix()
2196 sage: X = x.to_matrix()
2197 sage: Y = y.to_matrix()
2198 sage: expected = (X*Y + Y*X)/2
2199 sage: actual == expected
2201 sage: J(expected) == x*y
2204 We can change the generator prefix::
2206 sage: ComplexHermitianEJA(2, prefix='z').gens()
2209 We can construct the (trivial) algebra of rank zero::
2211 sage: ComplexHermitianEJA(0)
2212 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2217 def _denormalized_basis(cls
, n
):
2219 Returns a basis for the space of complex Hermitian n-by-n matrices.
2221 Why do we embed these? Basically, because all of numerical linear
2222 algebra assumes that you're working with vectors consisting of `n`
2223 entries from a field and scalars from the same field. There's no way
2224 to tell SageMath that (for example) the vectors contain complex
2225 numbers, while the scalar field is real.
2229 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
2233 sage: set_random_seed()
2234 sage: n = ZZ.random_element(1,5)
2235 sage: B = ComplexHermitianEJA._denormalized_basis(n)
2236 sage: all( M.is_symmetric() for M in B)
2241 R
= PolynomialRing(field
, 'z')
2243 F
= field
.extension(z
**2 + 1, 'I')
2246 # This is like the symmetric case, but we need to be careful:
2248 # * We want conjugate-symmetry, not just symmetry.
2249 # * The diagonal will (as a result) be real.
2252 Eij
= matrix
.zero(F
,n
)
2254 for j
in range(i
+1):
2258 Sij
= cls
.real_embed(Eij
)
2261 # The second one has a minus because it's conjugated.
2262 Eij
[j
,i
] = 1 # Eij = Eij + Eij.transpose()
2263 Sij_real
= cls
.real_embed(Eij
)
2265 # Eij = I*Eij - I*Eij.transpose()
2268 Sij_imag
= cls
.real_embed(Eij
)
2274 # Since we embedded these, we can drop back to the "field" that we
2275 # started with instead of the complex extension "F".
2276 return tuple( s
.change_ring(field
) for s
in S
)
2279 def __init__(self
, n
, **kwargs
):
2280 # We know this is a valid EJA, but will double-check
2281 # if the user passes check_axioms=True.
2282 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2288 super().__init
__(self
._denormalized
_basis
(n
),
2289 self
.jordan_product
,
2290 self
.trace_inner_product
,
2291 associative
=associative
,
2293 # TODO: this could be factored out somehow, but is left here
2294 # because the MatrixEJA is not presently a subclass of the
2295 # FDEJA class that defines rank() and one().
2296 self
.rank
.set_cache(n
)
2297 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
2298 self
.one
.set_cache(self(idV
))
2301 def _max_random_instance_size():
2302 return 3 # Dimension 9
2305 def random_instance(cls
, **kwargs
):
2307 Return a random instance of this type of algebra.
2309 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2310 return cls(n
, **kwargs
)
2312 class QuaternionMatrixEJA(MatrixEJA
):
2314 # A manual dictionary-cache for the quaternion_extension() method,
2315 # since apparently @classmethods can't also be @cached_methods.
2316 _quaternion_extension
= {}
2319 def quaternion_extension(cls
,field
):
2321 The quaternion field that we embed/unembed, as an extension
2322 of the given ``field``.
2324 if field
in cls
._quaternion
_extension
:
2325 return cls
._quaternion
_extension
[field
]
2327 Q
= QuaternionAlgebra(field
,-1,-1)
2329 cls
._quaternion
_extension
[field
] = Q
2333 def dimension_over_reals():
2337 def real_embed(cls
,M
):
2339 Embed the n-by-n quaternion matrix ``M`` into the space of real
2340 matrices of size 4n-by-4n by first sending each quaternion entry `z
2341 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
2342 c+di],[-c + di, a-bi]]`, and then embedding those into a real
2347 sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA
2351 sage: Q = QuaternionAlgebra(QQ,-1,-1)
2352 sage: i,j,k = Q.gens()
2353 sage: x = 1 + 2*i + 3*j + 4*k
2354 sage: M = matrix(Q, 1, [[x]])
2355 sage: QuaternionMatrixEJA.real_embed(M)
2361 Embedding is a homomorphism (isomorphism, in fact)::
2363 sage: set_random_seed()
2364 sage: n = ZZ.random_element(2)
2365 sage: Q = QuaternionAlgebra(QQ,-1,-1)
2366 sage: X = random_matrix(Q, n)
2367 sage: Y = random_matrix(Q, n)
2368 sage: Xe = QuaternionMatrixEJA.real_embed(X)
2369 sage: Ye = QuaternionMatrixEJA.real_embed(Y)
2370 sage: XYe = QuaternionMatrixEJA.real_embed(X*Y)
2375 super().real_embed(M
)
2376 quaternions
= M
.base_ring()
2379 F
= QuadraticField(-1, 'I')
2384 t
= z
.coefficient_tuple()
2389 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
2390 [-c
+ d
*i
, a
- b
*i
]])
2391 realM
= ComplexMatrixEJA
.real_embed(cplxM
)
2392 blocks
.append(realM
)
2394 # We should have real entries by now, so use the realest field
2395 # we've got for the return value.
2396 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
2401 def real_unembed(cls
,M
):
2403 The inverse of _embed_quaternion_matrix().
2407 sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA
2411 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
2412 ....: [-2, 1, -4, 3],
2413 ....: [-3, 4, 1, -2],
2414 ....: [-4, -3, 2, 1]])
2415 sage: QuaternionMatrixEJA.real_unembed(M)
2416 [1 + 2*i + 3*j + 4*k]
2420 Unembedding is the inverse of embedding::
2422 sage: set_random_seed()
2423 sage: Q = QuaternionAlgebra(QQ, -1, -1)
2424 sage: M = random_matrix(Q, 3)
2425 sage: Me = QuaternionMatrixEJA.real_embed(M)
2426 sage: QuaternionMatrixEJA.real_unembed(Me) == M
2430 super().real_unembed(M
)
2432 d
= cls
.dimension_over_reals()
2434 # Use the base ring of the matrix to ensure that its entries can be
2435 # multiplied by elements of the quaternion algebra.
2436 Q
= cls
.quaternion_extension(M
.base_ring())
2439 # Go top-left to bottom-right (reading order), converting every
2440 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
2443 for l
in range(n
/d
):
2444 for m
in range(n
/d
):
2445 submat
= ComplexMatrixEJA
.real_unembed(
2446 M
[d
*l
:d
*l
+d
,d
*m
:d
*m
+d
] )
2447 if submat
[0,0] != submat
[1,1].conjugate():
2448 raise ValueError('bad on-diagonal submatrix')
2449 if submat
[0,1] != -submat
[1,0].conjugate():
2450 raise ValueError('bad off-diagonal submatrix')
2451 z
= submat
[0,0].real()
2452 z
+= submat
[0,0].imag()*i
2453 z
+= submat
[0,1].real()*j
2454 z
+= submat
[0,1].imag()*k
2457 return matrix(Q
, n
/d
, elements
)
2460 class QuaternionHermitianEJA(ConcreteEJA
, QuaternionMatrixEJA
):
2462 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
2463 matrices, the usual symmetric Jordan product, and the
2464 real-part-of-trace inner product. It has dimension `2n^2 - n` over
2469 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2473 In theory, our "field" can be any subfield of the reals::
2475 sage: QuaternionHermitianEJA(2, field=RDF, check_axioms=True)
2476 Euclidean Jordan algebra of dimension 6 over Real Double Field
2477 sage: QuaternionHermitianEJA(2, field=RR, check_axioms=True)
2478 Euclidean Jordan algebra of dimension 6 over Real Field with
2479 53 bits of precision
2483 The dimension of this algebra is `2*n^2 - n`::
2485 sage: set_random_seed()
2486 sage: n_max = QuaternionHermitianEJA._max_random_instance_size()
2487 sage: n = ZZ.random_element(1, n_max)
2488 sage: J = QuaternionHermitianEJA(n)
2489 sage: J.dimension() == 2*(n^2) - n
2492 The Jordan multiplication is what we think it is::
2494 sage: set_random_seed()
2495 sage: J = QuaternionHermitianEJA.random_instance()
2496 sage: x,y = J.random_elements(2)
2497 sage: actual = (x*y).to_matrix()
2498 sage: X = x.to_matrix()
2499 sage: Y = y.to_matrix()
2500 sage: expected = (X*Y + Y*X)/2
2501 sage: actual == expected
2503 sage: J(expected) == x*y
2506 We can change the generator prefix::
2508 sage: QuaternionHermitianEJA(2, prefix='a').gens()
2509 (a0, a1, a2, a3, a4, a5)
2511 We can construct the (trivial) algebra of rank zero::
2513 sage: QuaternionHermitianEJA(0)
2514 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2518 def _denormalized_basis(cls
, n
):
2520 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
2522 Why do we embed these? Basically, because all of numerical
2523 linear algebra assumes that you're working with vectors consisting
2524 of `n` entries from a field and scalars from the same field. There's
2525 no way to tell SageMath that (for example) the vectors contain
2526 complex numbers, while the scalar field is real.
2530 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2534 sage: set_random_seed()
2535 sage: n = ZZ.random_element(1,5)
2536 sage: B = QuaternionHermitianEJA._denormalized_basis(n)
2537 sage: all( M.is_symmetric() for M in B )
2542 Q
= QuaternionAlgebra(QQ
,-1,-1)
2545 # This is like the symmetric case, but we need to be careful:
2547 # * We want conjugate-symmetry, not just symmetry.
2548 # * The diagonal will (as a result) be real.
2551 Eij
= matrix
.zero(Q
,n
)
2553 for j
in range(i
+1):
2557 Sij
= cls
.real_embed(Eij
)
2560 # The second, third, and fourth ones have a minus
2561 # because they're conjugated.
2562 # Eij = Eij + Eij.transpose()
2564 Sij_real
= cls
.real_embed(Eij
)
2566 # Eij = I*(Eij - Eij.transpose())
2569 Sij_I
= cls
.real_embed(Eij
)
2571 # Eij = J*(Eij - Eij.transpose())
2574 Sij_J
= cls
.real_embed(Eij
)
2576 # Eij = K*(Eij - Eij.transpose())
2579 Sij_K
= cls
.real_embed(Eij
)
2585 # Since we embedded these, we can drop back to the "field" that we
2586 # started with instead of the quaternion algebra "Q".
2587 return tuple( s
.change_ring(field
) for s
in S
)
2590 def __init__(self
, n
, **kwargs
):
2591 # We know this is a valid EJA, but will double-check
2592 # if the user passes check_axioms=True.
2593 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2599 super().__init
__(self
._denormalized
_basis
(n
),
2600 self
.jordan_product
,
2601 self
.trace_inner_product
,
2602 associative
=associative
,
2605 # TODO: this could be factored out somehow, but is left here
2606 # because the MatrixEJA is not presently a subclass of the
2607 # FDEJA class that defines rank() and one().
2608 self
.rank
.set_cache(n
)
2609 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
2610 self
.one
.set_cache(self(idV
))
2614 def _max_random_instance_size():
2616 The maximum rank of a random QuaternionHermitianEJA.
2618 return 2 # Dimension 6
2621 def random_instance(cls
, **kwargs
):
2623 Return a random instance of this type of algebra.
2625 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2626 return cls(n
, **kwargs
)
2629 class HadamardEJA(ConcreteEJA
):
2631 Return the Euclidean Jordan Algebra corresponding to the set
2632 `R^n` under the Hadamard product.
2634 Note: this is nothing more than the Cartesian product of ``n``
2635 copies of the spin algebra. Once Cartesian product algebras
2636 are implemented, this can go.
2640 sage: from mjo.eja.eja_algebra import HadamardEJA
2644 This multiplication table can be verified by hand::
2646 sage: J = HadamardEJA(3)
2647 sage: e0,e1,e2 = J.gens()
2663 We can change the generator prefix::
2665 sage: HadamardEJA(3, prefix='r').gens()
2669 def __init__(self
, n
, **kwargs
):
2671 jordan_product
= lambda x
,y
: x
2672 inner_product
= lambda x
,y
: x
2674 def jordan_product(x
,y
):
2676 return P( xi
*yi
for (xi
,yi
) in zip(x
,y
) )
2678 def inner_product(x
,y
):
2681 # New defaults for keyword arguments. Don't orthonormalize
2682 # because our basis is already orthonormal with respect to our
2683 # inner-product. Don't check the axioms, because we know this
2684 # is a valid EJA... but do double-check if the user passes
2685 # check_axioms=True. Note: we DON'T override the "check_field"
2686 # default here, because the user can pass in a field!
2687 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2688 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2690 column_basis
= tuple( b
.column() for b
in FreeModule(ZZ
, n
).basis() )
2691 super().__init
__(column_basis
,
2696 self
.rank
.set_cache(n
)
2699 self
.one
.set_cache( self
.zero() )
2701 self
.one
.set_cache( sum(self
.gens()) )
2704 def _max_random_instance_size():
2706 The maximum dimension of a random HadamardEJA.
2711 def random_instance(cls
, **kwargs
):
2713 Return a random instance of this type of algebra.
2715 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2716 return cls(n
, **kwargs
)
2719 class BilinearFormEJA(ConcreteEJA
):
2721 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2722 with the half-trace inner product and jordan product ``x*y =
2723 (<Bx,y>,y_bar>, x0*y_bar + y0*x_bar)`` where `B = 1 \times B22` is
2724 a symmetric positive-definite "bilinear form" matrix. Its
2725 dimension is the size of `B`, and it has rank two in dimensions
2726 larger than two. It reduces to the ``JordanSpinEJA`` when `B` is
2727 the identity matrix of order ``n``.
2729 We insist that the one-by-one upper-left identity block of `B` be
2730 passed in as well so that we can be passed a matrix of size zero
2731 to construct a trivial algebra.
2735 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
2736 ....: JordanSpinEJA)
2740 When no bilinear form is specified, the identity matrix is used,
2741 and the resulting algebra is the Jordan spin algebra::
2743 sage: B = matrix.identity(AA,3)
2744 sage: J0 = BilinearFormEJA(B)
2745 sage: J1 = JordanSpinEJA(3)
2746 sage: J0.multiplication_table() == J0.multiplication_table()
2749 An error is raised if the matrix `B` does not correspond to a
2750 positive-definite bilinear form::
2752 sage: B = matrix.random(QQ,2,3)
2753 sage: J = BilinearFormEJA(B)
2754 Traceback (most recent call last):
2756 ValueError: bilinear form is not positive-definite
2757 sage: B = matrix.zero(QQ,3)
2758 sage: J = BilinearFormEJA(B)
2759 Traceback (most recent call last):
2761 ValueError: bilinear form is not positive-definite
2765 We can create a zero-dimensional algebra::
2767 sage: B = matrix.identity(AA,0)
2768 sage: J = BilinearFormEJA(B)
2772 We can check the multiplication condition given in the Jordan, von
2773 Neumann, and Wigner paper (and also discussed on my "On the
2774 symmetry..." paper). Note that this relies heavily on the standard
2775 choice of basis, as does anything utilizing the bilinear form
2776 matrix. We opt not to orthonormalize the basis, because if we
2777 did, we would have to normalize the `s_{i}` in a similar manner::
2779 sage: set_random_seed()
2780 sage: n = ZZ.random_element(5)
2781 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2782 sage: B11 = matrix.identity(QQ,1)
2783 sage: B22 = M.transpose()*M
2784 sage: B = block_matrix(2,2,[ [B11,0 ],
2786 sage: J = BilinearFormEJA(B, orthonormalize=False)
2787 sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
2788 sage: V = J.vector_space()
2789 sage: sis = [ J( V([0] + (M.inverse()*ei).list()).column() )
2790 ....: for ei in eis ]
2791 sage: actual = [ sis[i]*sis[j]
2792 ....: for i in range(n-1)
2793 ....: for j in range(n-1) ]
2794 sage: expected = [ J.one() if i == j else J.zero()
2795 ....: for i in range(n-1)
2796 ....: for j in range(n-1) ]
2797 sage: actual == expected
2801 def __init__(self
, B
, **kwargs
):
2802 # The matrix "B" is supplied by the user in most cases,
2803 # so it makes sense to check whether or not its positive-
2804 # definite unless we are specifically asked not to...
2805 if ("check_axioms" not in kwargs
) or kwargs
["check_axioms"]:
2806 if not B
.is_positive_definite():
2807 raise ValueError("bilinear form is not positive-definite")
2809 # However, all of the other data for this EJA is computed
2810 # by us in manner that guarantees the axioms are
2811 # satisfied. So, again, unless we are specifically asked to
2812 # verify things, we'll skip the rest of the checks.
2813 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2815 def inner_product(x
,y
):
2816 return (y
.T
*B
*x
)[0,0]
2818 def jordan_product(x
,y
):
2824 z0
= inner_product(y
,x
)
2825 zbar
= y0
*xbar
+ x0
*ybar
2826 return P([z0
] + zbar
.list())
2829 column_basis
= tuple( b
.column() for b
in FreeModule(ZZ
, n
).basis() )
2831 # TODO: I haven't actually checked this, but it seems legit.
2836 super().__init
__(column_basis
,
2839 associative
=associative
,
2842 # The rank of this algebra is two, unless we're in a
2843 # one-dimensional ambient space (because the rank is bounded
2844 # by the ambient dimension).
2845 self
.rank
.set_cache(min(n
,2))
2848 self
.one
.set_cache( self
.zero() )
2850 self
.one
.set_cache( self
.monomial(0) )
2853 def _max_random_instance_size():
2855 The maximum dimension of a random BilinearFormEJA.
2860 def random_instance(cls
, **kwargs
):
2862 Return a random instance of this algebra.
2864 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2866 B
= matrix
.identity(ZZ
, n
)
2867 return cls(B
, **kwargs
)
2869 B11
= matrix
.identity(ZZ
, 1)
2870 M
= matrix
.random(ZZ
, n
-1)
2871 I
= matrix
.identity(ZZ
, n
-1)
2873 while alpha
.is_zero():
2874 alpha
= ZZ
.random_element().abs()
2875 B22
= M
.transpose()*M
+ alpha
*I
2877 from sage
.matrix
.special
import block_matrix
2878 B
= block_matrix(2,2, [ [B11
, ZZ(0) ],
2881 return cls(B
, **kwargs
)
2884 class JordanSpinEJA(BilinearFormEJA
):
2886 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2887 with the usual inner product and jordan product ``x*y =
2888 (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2893 sage: from mjo.eja.eja_algebra import JordanSpinEJA
2897 This multiplication table can be verified by hand::
2899 sage: J = JordanSpinEJA(4)
2900 sage: e0,e1,e2,e3 = J.gens()
2916 We can change the generator prefix::
2918 sage: JordanSpinEJA(2, prefix='B').gens()
2923 Ensure that we have the usual inner product on `R^n`::
2925 sage: set_random_seed()
2926 sage: J = JordanSpinEJA.random_instance()
2927 sage: x,y = J.random_elements(2)
2928 sage: actual = x.inner_product(y)
2929 sage: expected = x.to_vector().inner_product(y.to_vector())
2930 sage: actual == expected
2934 def __init__(self
, n
, **kwargs
):
2935 # This is a special case of the BilinearFormEJA with the
2936 # identity matrix as its bilinear form.
2937 B
= matrix
.identity(ZZ
, n
)
2939 # Don't orthonormalize because our basis is already
2940 # orthonormal with respect to our inner-product.
2941 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2943 # But also don't pass check_field=False here, because the user
2944 # can pass in a field!
2945 super().__init
__(B
, **kwargs
)
2948 def _max_random_instance_size():
2950 The maximum dimension of a random JordanSpinEJA.
2955 def random_instance(cls
, **kwargs
):
2957 Return a random instance of this type of algebra.
2959 Needed here to override the implementation for ``BilinearFormEJA``.
2961 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2962 return cls(n
, **kwargs
)
2965 class TrivialEJA(ConcreteEJA
):
2967 The trivial Euclidean Jordan algebra consisting of only a zero element.
2971 sage: from mjo.eja.eja_algebra import TrivialEJA
2975 sage: J = TrivialEJA()
2982 sage: 7*J.one()*12*J.one()
2984 sage: J.one().inner_product(J.one())
2986 sage: J.one().norm()
2988 sage: J.one().subalgebra_generated_by()
2989 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2994 def __init__(self
, **kwargs
):
2995 jordan_product
= lambda x
,y
: x
2996 inner_product
= lambda x
,y
: 0
2999 # New defaults for keyword arguments
3000 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
3001 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
3003 super().__init
__(basis
,
3009 # The rank is zero using my definition, namely the dimension of the
3010 # largest subalgebra generated by any element.
3011 self
.rank
.set_cache(0)
3012 self
.one
.set_cache( self
.zero() )
3015 def random_instance(cls
, **kwargs
):
3016 # We don't take a "size" argument so the superclass method is
3017 # inappropriate for us.
3018 return cls(**kwargs
)
3021 class CartesianProductEJA(CombinatorialFreeModule_CartesianProduct
,
3022 FiniteDimensionalEJA
):
3024 The external (orthogonal) direct sum of two or more Euclidean
3025 Jordan algebras. Every Euclidean Jordan algebra decomposes into an
3026 orthogonal direct sum of simple Euclidean Jordan algebras which is
3027 then isometric to a Cartesian product, so no generality is lost by
3028 providing only this construction.
3032 sage: from mjo.eja.eja_algebra import (random_eja,
3033 ....: CartesianProductEJA,
3035 ....: JordanSpinEJA,
3036 ....: RealSymmetricEJA)
3040 The Jordan product is inherited from our factors and implemented by
3041 our CombinatorialFreeModule Cartesian product superclass::
3043 sage: set_random_seed()
3044 sage: J1 = HadamardEJA(2)
3045 sage: J2 = RealSymmetricEJA(2)
3046 sage: J = cartesian_product([J1,J2])
3047 sage: x,y = J.random_elements(2)
3051 The ability to retrieve the original factors is implemented by our
3052 CombinatorialFreeModule Cartesian product superclass::
3054 sage: J1 = HadamardEJA(2, field=QQ)
3055 sage: J2 = JordanSpinEJA(3, field=QQ)
3056 sage: J = cartesian_product([J1,J2])
3057 sage: J.cartesian_factors()
3058 (Euclidean Jordan algebra of dimension 2 over Rational Field,
3059 Euclidean Jordan algebra of dimension 3 over Rational Field)
3061 You can provide more than two factors::
3063 sage: J1 = HadamardEJA(2)
3064 sage: J2 = JordanSpinEJA(3)
3065 sage: J3 = RealSymmetricEJA(3)
3066 sage: cartesian_product([J1,J2,J3])
3067 Euclidean Jordan algebra of dimension 2 over Algebraic Real
3068 Field (+) Euclidean Jordan algebra of dimension 3 over Algebraic
3069 Real Field (+) Euclidean Jordan algebra of dimension 6 over
3070 Algebraic Real Field
3072 Rank is additive on a Cartesian product::
3074 sage: J1 = HadamardEJA(1)
3075 sage: J2 = RealSymmetricEJA(2)
3076 sage: J = cartesian_product([J1,J2])
3077 sage: J1.rank.clear_cache()
3078 sage: J2.rank.clear_cache()
3079 sage: J.rank.clear_cache()
3082 sage: J.rank() == J1.rank() + J2.rank()
3085 The same rank computation works over the rationals, with whatever
3088 sage: J1 = HadamardEJA(1, field=QQ, orthonormalize=False)
3089 sage: J2 = RealSymmetricEJA(2, field=QQ, orthonormalize=False)
3090 sage: J = cartesian_product([J1,J2])
3091 sage: J1.rank.clear_cache()
3092 sage: J2.rank.clear_cache()
3093 sage: J.rank.clear_cache()
3096 sage: J.rank() == J1.rank() + J2.rank()
3099 The product algebra will be associative if and only if all of its
3100 components are associative::
3102 sage: J1 = HadamardEJA(2)
3103 sage: J1.is_associative()
3105 sage: J2 = HadamardEJA(3)
3106 sage: J2.is_associative()
3108 sage: J3 = RealSymmetricEJA(3)
3109 sage: J3.is_associative()
3111 sage: CP1 = cartesian_product([J1,J2])
3112 sage: CP1.is_associative()
3114 sage: CP2 = cartesian_product([J1,J3])
3115 sage: CP2.is_associative()
3120 All factors must share the same base field::
3122 sage: J1 = HadamardEJA(2, field=QQ)
3123 sage: J2 = RealSymmetricEJA(2)
3124 sage: CartesianProductEJA((J1,J2))
3125 Traceback (most recent call last):
3127 ValueError: all factors must share the same base field
3129 The cached unit element is the same one that would be computed::
3131 sage: set_random_seed() # long time
3132 sage: J1 = random_eja() # long time
3133 sage: J2 = random_eja() # long time
3134 sage: J = cartesian_product([J1,J2]) # long time
3135 sage: actual = J.one() # long time
3136 sage: J.one.clear_cache() # long time
3137 sage: expected = J.one() # long time
3138 sage: actual == expected # long time
3142 Element
= FiniteDimensionalEJAElement
3145 def __init__(self
, algebras
, **kwargs
):
3146 CombinatorialFreeModule_CartesianProduct
.__init
__(self
,
3149 field
= algebras
[0].base_ring()
3150 if not all( J
.base_ring() == field
for J
in algebras
):
3151 raise ValueError("all factors must share the same base field")
3153 associative
= all( m
.is_associative() for m
in algebras
)
3155 # The definition of matrix_space() and self.basis() relies
3156 # only on the stuff in the CFM_CartesianProduct class, which
3157 # we've already initialized.
3158 Js
= self
.cartesian_factors()
3160 MS
= self
.matrix_space()
3162 MS(tuple( self
.cartesian_projection(i
)(b
).to_matrix()
3163 for i
in range(m
) ))
3164 for b
in self
.basis()
3167 # Define jordan/inner products that operate on that matrix_basis.
3168 def jordan_product(x
,y
):
3170 (Js
[i
](x
[i
])*Js
[i
](y
[i
])).to_matrix() for i
in range(m
)
3173 def inner_product(x
, y
):
3175 Js
[i
](x
[i
]).inner_product(Js
[i
](y
[i
])) for i
in range(m
)
3178 # There's no need to check the field since it already came
3179 # from an EJA. Likewise the axioms are guaranteed to be
3180 # satisfied, unless the guy writing this class sucks.
3182 # If you want the basis to be orthonormalized, orthonormalize
3184 FiniteDimensionalEJA
.__init
__(self
,
3189 orthonormalize
=False,
3190 associative
=associative
,
3191 cartesian_product
=True,
3195 ones
= tuple(J
.one() for J
in algebras
)
3196 self
.one
.set_cache(self
._cartesian
_product
_of
_elements
(ones
))
3197 self
.rank
.set_cache(sum(J
.rank() for J
in algebras
))
3199 def _monomial_to_generator(self
, mon
):
3201 Convert a monomial index into a generator index.
3205 sage: from mjo.eja.eja_algebra import random_eja
3209 sage: J1 = random_eja(field=QQ, orthonormalize=False)
3210 sage: J2 = random_eja(field=QQ, orthonormalize=False)
3211 sage: J = cartesian_product([J1,J2])
3212 sage: all( J.monomial(m)
3214 ....: J.gens()[J._monomial_to_generator(m)]
3215 ....: for m in J.basis().keys() )
3219 # The superclass method indexes into a matrix, so we have to
3220 # turn the tuples i and j into integers. This is easy enough
3221 # given that the first coordinate of i and j corresponds to
3222 # the factor, and the second coordinate corresponds to the
3223 # index of the generator within that factor.
3226 except TypeError: # 'int' object is not subscriptable
3228 idx_in_factor
= self
._monomial
_to
_generator
(mon
[1])
3230 offset
= sum( f
.dimension()
3231 for f
in self
.cartesian_factors()[:factor
] )
3232 return offset
+ idx_in_factor
3234 def product_on_basis(self
, i
, j
):
3236 Return the product of the monomials indexed by ``i`` and ``j``.
3238 This overrides the superclass method because here, both ``i``
3239 and ``j`` will be ordered pairs.
3243 sage: from mjo.eja.eja_algebra import (HadamardEJA,
3244 ....: JordanSpinEJA,
3245 ....: QuaternionHermitianEJA,
3246 ....: RealSymmetricEJA,)
3250 sage: J1 = JordanSpinEJA(2, field=QQ)
3251 sage: J2 = RealSymmetricEJA(2, field=QQ, orthonormalize=False)
3252 sage: J3 = HadamardEJA(1, field=QQ)
3253 sage: K1 = cartesian_product([J1,J2])
3254 sage: K2 = cartesian_product([K1,J3])
3255 sage: list(K2.basis())
3256 [e(0, (0, 0)), e(0, (0, 1)), e(0, (1, 0)), e(0, (1, 1)),
3257 e(0, (1, 2)), e(1, 0)]
3259 sage: (g[0] + 2*g[3]) * (g[1] - 4*g[2])
3260 e(0, (0, 1)) - 4*e(0, (1, 1))
3264 sage: J1 = RealSymmetricEJA(1,field=QQ)
3265 sage: J2 = QuaternionHermitianEJA(1,field=QQ)
3266 sage: J = cartesian_product([J1,J2])
3267 sage: x = sum(J.gens())
3274 l
= self
._monomial
_to
_generator
(i
)
3275 m
= self
._monomial
_to
_generator
(j
)
3276 return FiniteDimensionalEJA
.product_on_basis(self
, l
, m
)
3278 def matrix_space(self
):
3280 Return the space that our matrix basis lives in as a Cartesian
3285 sage: from mjo.eja.eja_algebra import (HadamardEJA,
3286 ....: RealSymmetricEJA)
3290 sage: J1 = HadamardEJA(1)
3291 sage: J2 = RealSymmetricEJA(2)
3292 sage: J = cartesian_product([J1,J2])
3293 sage: J.matrix_space()
3294 The Cartesian product of (Full MatrixSpace of 1 by 1 dense
3295 matrices over Algebraic Real Field, Full MatrixSpace of 2
3296 by 2 dense matrices over Algebraic Real Field)
3299 from sage
.categories
.cartesian_product
import cartesian_product
3300 return cartesian_product( [J
.matrix_space()
3301 for J
in self
.cartesian_factors()] )
3304 def cartesian_projection(self
, i
):
3308 sage: from mjo.eja.eja_algebra import (random_eja,
3309 ....: JordanSpinEJA,
3311 ....: RealSymmetricEJA,
3312 ....: ComplexHermitianEJA)
3316 The projection morphisms are Euclidean Jordan algebra
3319 sage: J1 = HadamardEJA(2)
3320 sage: J2 = RealSymmetricEJA(2)
3321 sage: J = cartesian_product([J1,J2])
3322 sage: J.cartesian_projection(0)
3323 Linear operator between finite-dimensional Euclidean Jordan
3324 algebras represented by the matrix:
3327 Domain: Euclidean Jordan algebra of dimension 2 over Algebraic
3328 Real Field (+) Euclidean Jordan algebra of dimension 3 over
3329 Algebraic Real Field
3330 Codomain: Euclidean Jordan algebra of dimension 2 over Algebraic
3332 sage: J.cartesian_projection(1)
3333 Linear operator between finite-dimensional Euclidean Jordan
3334 algebras represented by the matrix:
3338 Domain: Euclidean Jordan algebra of dimension 2 over Algebraic
3339 Real Field (+) Euclidean Jordan algebra of dimension 3 over
3340 Algebraic Real Field
3341 Codomain: Euclidean Jordan algebra of dimension 3 over Algebraic
3344 The projections work the way you'd expect on the vector
3345 representation of an element::
3347 sage: J1 = JordanSpinEJA(2)
3348 sage: J2 = ComplexHermitianEJA(2)
3349 sage: J = cartesian_product([J1,J2])
3350 sage: pi_left = J.cartesian_projection(0)
3351 sage: pi_right = J.cartesian_projection(1)
3352 sage: pi_left(J.one()).to_vector()
3354 sage: pi_right(J.one()).to_vector()
3356 sage: J.one().to_vector()
3361 The answer never changes::
3363 sage: set_random_seed()
3364 sage: J1 = random_eja()
3365 sage: J2 = random_eja()
3366 sage: J = cartesian_product([J1,J2])
3367 sage: P0 = J.cartesian_projection(0)
3368 sage: P1 = J.cartesian_projection(0)
3373 Ji
= self
.cartesian_factors()[i
]
3374 # Requires the fix on Trac 31421/31422 to work!
3375 Pi
= super().cartesian_projection(i
)
3376 return FiniteDimensionalEJAOperator(self
,Ji
,Pi
.matrix())
3379 def cartesian_embedding(self
, i
):
3383 sage: from mjo.eja.eja_algebra import (random_eja,
3384 ....: JordanSpinEJA,
3386 ....: RealSymmetricEJA)
3390 The embedding morphisms are Euclidean Jordan algebra
3393 sage: J1 = HadamardEJA(2)
3394 sage: J2 = RealSymmetricEJA(2)
3395 sage: J = cartesian_product([J1,J2])
3396 sage: J.cartesian_embedding(0)
3397 Linear operator between finite-dimensional Euclidean Jordan
3398 algebras represented by the matrix:
3404 Domain: Euclidean Jordan algebra of dimension 2 over
3405 Algebraic Real Field
3406 Codomain: Euclidean Jordan algebra of dimension 2 over
3407 Algebraic Real Field (+) Euclidean Jordan algebra of
3408 dimension 3 over Algebraic Real Field
3409 sage: J.cartesian_embedding(1)
3410 Linear operator between finite-dimensional Euclidean Jordan
3411 algebras represented by the matrix:
3417 Domain: Euclidean Jordan algebra of dimension 3 over
3418 Algebraic Real Field
3419 Codomain: Euclidean Jordan algebra of dimension 2 over
3420 Algebraic Real Field (+) Euclidean Jordan algebra of
3421 dimension 3 over Algebraic Real Field
3423 The embeddings work the way you'd expect on the vector
3424 representation of an element::
3426 sage: J1 = JordanSpinEJA(3)
3427 sage: J2 = RealSymmetricEJA(2)
3428 sage: J = cartesian_product([J1,J2])
3429 sage: iota_left = J.cartesian_embedding(0)
3430 sage: iota_right = J.cartesian_embedding(1)
3431 sage: iota_left(J1.zero()) == J.zero()
3433 sage: iota_right(J2.zero()) == J.zero()
3435 sage: J1.one().to_vector()
3437 sage: iota_left(J1.one()).to_vector()
3439 sage: J2.one().to_vector()
3441 sage: iota_right(J2.one()).to_vector()
3443 sage: J.one().to_vector()
3448 The answer never changes::
3450 sage: set_random_seed()
3451 sage: J1 = random_eja()
3452 sage: J2 = random_eja()
3453 sage: J = cartesian_product([J1,J2])
3454 sage: E0 = J.cartesian_embedding(0)
3455 sage: E1 = J.cartesian_embedding(0)
3459 Composing a projection with the corresponding inclusion should
3460 produce the identity map, and mismatching them should produce
3463 sage: set_random_seed()
3464 sage: J1 = random_eja()
3465 sage: J2 = random_eja()
3466 sage: J = cartesian_product([J1,J2])
3467 sage: iota_left = J.cartesian_embedding(0)
3468 sage: iota_right = J.cartesian_embedding(1)
3469 sage: pi_left = J.cartesian_projection(0)
3470 sage: pi_right = J.cartesian_projection(1)
3471 sage: pi_left*iota_left == J1.one().operator()
3473 sage: pi_right*iota_right == J2.one().operator()
3475 sage: (pi_left*iota_right).is_zero()
3477 sage: (pi_right*iota_left).is_zero()
3481 Ji
= self
.cartesian_factors()[i
]
3482 # Requires the fix on Trac 31421/31422 to work!
3483 Ei
= super().cartesian_embedding(i
)
3484 return FiniteDimensionalEJAOperator(Ji
,self
,Ei
.matrix())
3488 FiniteDimensionalEJA
.CartesianProduct
= CartesianProductEJA
3490 class RationalBasisCartesianProductEJA(CartesianProductEJA
,
3493 A separate class for products of algebras for which we know a
3498 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
3499 ....: RealSymmetricEJA)
3503 This gives us fast characteristic polynomial computations in
3504 product algebras, too::
3507 sage: J1 = JordanSpinEJA(2)
3508 sage: J2 = RealSymmetricEJA(3)
3509 sage: J = cartesian_product([J1,J2])
3510 sage: J.characteristic_polynomial_of().degree()
3516 def __init__(self
, algebras
, **kwargs
):
3517 CartesianProductEJA
.__init
__(self
, algebras
, **kwargs
)
3519 self
._rational
_algebra
= None
3520 if self
.vector_space().base_field() is not QQ
:
3521 self
._rational
_algebra
= cartesian_product([
3522 r
._rational
_algebra
for r
in algebras
3526 RationalBasisEJA
.CartesianProduct
= RationalBasisCartesianProductEJA
3528 random_eja
= ConcreteEJA
.random_instance
3530 # def random_eja(*args, **kwargs):
3531 # J1 = ConcreteEJA.random_instance(*args, **kwargs)
3533 # # This might make Cartesian products appear roughly as often as
3534 # # any other ConcreteEJA.
3535 # if ZZ.random_element(len(ConcreteEJA.__subclasses__()) + 1) == 0:
3536 # # Use random_eja() again so we can get more than two factors.
3537 # J2 = random_eja(*args, **kwargs)
3538 # J = cartesian_product([J1,J2])