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`
35 * :class:`OctonionHermitianEJA`
37 In addition to these, we provide two other example constructions,
39 * :class:`JordanSpinEJA`
40 * :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. The Albert EJA is simply a special
47 case of the :class:`OctonionHermitianEJA` where the matrices are
48 three-by-three and the resulting space has dimension 27. And
49 last/least, the trivial EJA is exactly what you think it is; it could
50 also be obtained by constructing a dimension-zero instance of any of
51 the other algebras. Cartesian products of these are also supported
52 using the usual ``cartesian_product()`` function; as a result, we
53 support (up to isomorphism) all Euclidean Jordan algebras.
55 At a minimum, the following are required to construct a Euclidean
58 * A basis of matrices, column vectors, or MatrixAlgebra elements
59 * A Jordan product defined on the basis
60 * Its inner product defined on the basis
62 The real numbers form a Euclidean Jordan algebra when both the Jordan
63 and inner products are the usual multiplication. We use this as our
64 example, and demonstrate a few ways to construct an EJA.
66 First, we can use one-by-one SageMath matrices with algebraic real
67 entries to represent real numbers. We define the Jordan and inner
68 products to be essentially real-number multiplication, with the only
69 difference being that the Jordan product again returns a one-by-one
70 matrix, whereas the inner product must return a scalar. Our basis for
71 the one-by-one matrices is of course the set consisting of a single
72 matrix with its sole entry non-zero::
74 sage: from mjo.eja.eja_algebra import FiniteDimensionalEJA
75 sage: jp = lambda X,Y: X*Y
76 sage: ip = lambda X,Y: X[0,0]*Y[0,0]
77 sage: b1 = matrix(AA, [[1]])
78 sage: J1 = FiniteDimensionalEJA((b1,), jp, ip)
80 Euclidean Jordan algebra of dimension 1 over Algebraic Real Field
82 In fact, any positive scalar multiple of that inner-product would work::
84 sage: ip2 = lambda X,Y: 16*ip(X,Y)
85 sage: J2 = FiniteDimensionalEJA((b1,), jp, ip2)
87 Euclidean Jordan algebra of dimension 1 over Algebraic Real Field
89 But beware that your basis will be orthonormalized _with respect to the
90 given inner-product_ unless you pass ``orthonormalize=False`` to the
91 constructor. For example::
93 sage: J3 = FiniteDimensionalEJA((b1,), jp, ip2, orthonormalize=False)
95 Euclidean Jordan algebra of dimension 1 over Algebraic Real Field
97 To see the difference, you can take the first and only basis element
98 of the resulting algebra, and ask for it to be converted back into
101 sage: J1.basis()[0].to_matrix()
103 sage: J2.basis()[0].to_matrix()
105 sage: J3.basis()[0].to_matrix()
108 Since square roots are used in that process, the default scalar field
109 that we use is the field of algebraic real numbers, ``AA``. You can
110 also Use rational numbers, but only if you either pass
111 ``orthonormalize=False`` or know that orthonormalizing your basis
112 won't stray beyond the rational numbers. The example above would
113 have worked only because ``sqrt(16) == 4`` is rational.
115 Another option for your basis is to use elemebts of a
116 :class:`MatrixAlgebra`::
118 sage: from mjo.matrix_algebra import MatrixAlgebra
119 sage: A = MatrixAlgebra(1,AA,AA)
120 sage: J4 = FiniteDimensionalEJA(A.gens(), jp, ip)
122 Euclidean Jordan algebra of dimension 1 over Algebraic Real Field
123 sage: J4.basis()[0].to_matrix()
128 An easier way to view the entire EJA basis in its original (but
129 perhaps orthonormalized) matrix form is to use the ``matrix_basis``
132 sage: J4.matrix_basis()
137 In particular, a :class:`MatrixAlgebra` is needed to work around the
138 fact that matrices in SageMath must have entries in the same
139 (commutative and associative) ring as its scalars. There are many
140 Euclidean Jordan algebras whose elements are matrices that violate
141 those assumptions. The complex, quaternion, and octonion Hermitian
142 matrices all have entries in a ring (the complex numbers, quaternions,
143 or octonions...) that differs from the algebra's scalar ring (the real
144 numbers). Quaternions are also non-commutative; the octonions are
145 neither commutative nor associative.
149 sage: from mjo.eja.eja_algebra import random_eja
154 Euclidean Jordan algebra of dimension...
157 from sage
.algebras
.quatalg
.quaternion_algebra
import QuaternionAlgebra
158 from sage
.categories
.magmatic_algebras
import MagmaticAlgebras
159 from sage
.categories
.sets_cat
import cartesian_product
160 from sage
.combinat
.free_module
import CombinatorialFreeModule
161 from sage
.matrix
.constructor
import matrix
162 from sage
.matrix
.matrix_space
import MatrixSpace
163 from sage
.misc
.cachefunc
import cached_method
164 from sage
.misc
.table
import table
165 from sage
.modules
.free_module
import FreeModule
, VectorSpace
166 from sage
.rings
.all
import (ZZ
, QQ
, AA
, QQbar
, RR
, RLF
, CLF
,
169 from mjo
.eja
.eja_element
import (CartesianProductEJAElement
,
170 FiniteDimensionalEJAElement
)
171 from mjo
.eja
.eja_operator
import FiniteDimensionalEJAOperator
172 from mjo
.eja
.eja_utils
import _all2list
174 def EuclideanJordanAlgebras(field
):
176 The category of Euclidean Jordan algebras over ``field``, which
177 must be a subfield of the real numbers. For now this is just a
178 convenient wrapper around all of the other category axioms that
181 category
= MagmaticAlgebras(field
).FiniteDimensional()
182 category
= category
.WithBasis().Unital().Commutative()
185 class FiniteDimensionalEJA(CombinatorialFreeModule
):
187 A finite-dimensional Euclidean Jordan algebra.
191 - ``basis`` -- a tuple; a tuple of basis elements in "matrix
192 form," which must be the same form as the arguments to
193 ``jordan_product`` and ``inner_product``. In reality, "matrix
194 form" can be either vectors, matrices, or a Cartesian product
195 (ordered tuple) of vectors or matrices. All of these would
196 ideally be vector spaces in sage with no special-casing
197 needed; but in reality we turn vectors into column-matrices
198 and Cartesian products `(a,b)` into column matrices
199 `(a,b)^{T}` after converting `a` and `b` themselves.
201 - ``jordan_product`` -- a function; afunction of two ``basis``
202 elements (in matrix form) that returns their jordan product,
203 also in matrix form; this will be applied to ``basis`` to
204 compute a multiplication table for the algebra.
206 - ``inner_product`` -- a function; a function of two ``basis``
207 elements (in matrix form) that returns their inner
208 product. This will be applied to ``basis`` to compute an
209 inner-product table (basically a matrix) for this algebra.
211 - ``matrix_space`` -- the space that your matrix basis lives in,
212 or ``None`` (the default). So long as your basis does not have
213 length zero you can omit this. But in trivial algebras, it is
216 - ``field`` -- a subfield of the reals (default: ``AA``); the scalar
217 field for the algebra.
219 - ``orthonormalize`` -- boolean (default: ``True``); whether or
220 not to orthonormalize the basis. Doing so is expensive and
221 generally rules out using the rationals as your ``field``, but
222 is required for spectral decompositions.
226 sage: from mjo.eja.eja_algebra import random_eja
230 We should compute that an element subalgebra is associative even
231 if we circumvent the element method::
233 sage: set_random_seed()
234 sage: J = random_eja(field=QQ,orthonormalize=False)
235 sage: x = J.random_element()
236 sage: A = x.subalgebra_generated_by(orthonormalize=False)
237 sage: basis = tuple(b.superalgebra_element() for b in A.basis())
238 sage: J.subalgebra(basis, orthonormalize=False).is_associative()
241 Element
= FiniteDimensionalEJAElement
244 def _check_input_field(field
):
245 if not field
.is_subring(RR
):
246 # Note: this does return true for the real algebraic
247 # field, the rationals, and any quadratic field where
248 # we've specified a real embedding.
249 raise ValueError("scalar field is not real")
252 def _check_input_axioms(basis
, jordan_product
, inner_product
):
253 if not all( jordan_product(bi
,bj
) == jordan_product(bj
,bi
)
256 raise ValueError("Jordan product is not commutative")
258 if not all( inner_product(bi
,bj
) == inner_product(bj
,bi
)
261 raise ValueError("inner-product is not commutative")
278 self
._check
_input
_field
(field
)
281 # Check commutativity of the Jordan and inner-products.
282 # This has to be done before we build the multiplication
283 # and inner-product tables/matrices, because we take
284 # advantage of symmetry in the process.
285 self
._check
_input
_axioms
(basis
, jordan_product
, inner_product
)
288 # All zero- and one-dimensional algebras are just the real
289 # numbers with (some positive multiples of) the usual
290 # multiplication as its Jordan and inner-product.
292 if associative
is None:
293 # We should figure it out. As with check_axioms, we have to do
294 # this without the help of the _jordan_product_is_associative()
295 # method because we need to know the category before we
296 # initialize the algebra.
297 associative
= all( jordan_product(jordan_product(bi
,bj
),bk
)
299 jordan_product(bi
,jordan_product(bj
,bk
))
304 category
= EuclideanJordanAlgebras(field
)
307 # Element subalgebras can take advantage of this.
308 category
= category
.Associative()
310 # Call the superclass constructor so that we can use its from_vector()
311 # method to build our multiplication table.
312 CombinatorialFreeModule
.__init
__(self
,
319 # Now comes all of the hard work. We'll be constructing an
320 # ambient vector space V that our (vectorized) basis lives in,
321 # as well as a subspace W of V spanned by those (vectorized)
322 # basis elements. The W-coordinates are the coefficients that
323 # we see in things like x = 1*b1 + 2*b2.
327 degree
= len(_all2list(basis
[0]))
329 # Build an ambient space that fits our matrix basis when
330 # written out as "long vectors."
331 V
= VectorSpace(field
, degree
)
333 # The matrix that will hold the orthonormal -> unorthonormal
334 # coordinate transformation. Default to an identity matrix of
335 # the appropriate size to avoid special cases for None
337 self
._deortho
_matrix
= matrix
.identity(field
,n
)
340 # Save a copy of the un-orthonormalized basis for later.
341 # Convert it to ambient V (vector) coordinates while we're
342 # at it, because we'd have to do it later anyway.
343 deortho_vector_basis
= tuple( V(_all2list(b
)) for b
in basis
)
345 from mjo
.eja
.eja_utils
import gram_schmidt
346 basis
= tuple(gram_schmidt(basis
, inner_product
))
348 # Save the (possibly orthonormalized) matrix basis for
349 # later, as well as the space that its elements live in.
350 # In most cases we can deduce the matrix space, but when
351 # n == 0 (that is, there are no basis elements) we cannot.
352 self
._matrix
_basis
= basis
353 if matrix_space
is None:
354 self
._matrix
_space
= self
._matrix
_basis
[0].parent()
356 self
._matrix
_space
= matrix_space
358 # Now create the vector space for the algebra, which will have
359 # its own set of non-ambient coordinates (in terms of the
361 vector_basis
= tuple( V(_all2list(b
)) for b
in basis
)
363 # Save the span of our matrix basis (when written out as long
364 # vectors) because otherwise we'll have to reconstruct it
365 # every time we want to coerce a matrix into the algebra.
366 self
._matrix
_span
= V
.span_of_basis( vector_basis
, check
=check_axioms
)
369 # Now "self._matrix_span" is the vector space of our
370 # algebra coordinates. The variables "X0", "X1",... refer
371 # to the entries of vectors in self._matrix_span. Thus to
372 # convert back and forth between the orthonormal
373 # coordinates and the given ones, we need to stick the
374 # original basis in self._matrix_span.
375 U
= V
.span_of_basis( deortho_vector_basis
, check
=check_axioms
)
376 self
._deortho
_matrix
= matrix
.column( U
.coordinate_vector(q
)
377 for q
in vector_basis
)
380 # Now we actually compute the multiplication and inner-product
381 # tables/matrices using the possibly-orthonormalized basis.
382 self
._inner
_product
_matrix
= matrix
.identity(field
, n
)
384 self
._multiplication
_table
= [ [zed
for j
in range(i
+1)]
387 # Note: the Jordan and inner-products are defined in terms
388 # of the ambient basis. It's important that their arguments
389 # are in ambient coordinates as well.
392 # ortho basis w.r.t. ambient coords
396 # The jordan product returns a matrixy answer, so we
397 # have to convert it to the algebra coordinates.
398 elt
= jordan_product(q_i
, q_j
)
399 elt
= self
._matrix
_span
.coordinate_vector(V(_all2list(elt
)))
400 self
._multiplication
_table
[i
][j
] = self
.from_vector(elt
)
402 if not orthonormalize
:
403 # If we're orthonormalizing the basis with respect
404 # to an inner-product, then the inner-product
405 # matrix with respect to the resulting basis is
406 # just going to be the identity.
407 ip
= inner_product(q_i
, q_j
)
408 self
._inner
_product
_matrix
[i
,j
] = ip
409 self
._inner
_product
_matrix
[j
,i
] = ip
411 self
._inner
_product
_matrix
._cache
= {'hermitian': True}
412 self
._inner
_product
_matrix
.set_immutable()
415 if not self
._is
_jordanian
():
416 raise ValueError("Jordan identity does not hold")
417 if not self
._inner
_product
_is
_associative
():
418 raise ValueError("inner product is not associative")
421 def _coerce_map_from_base_ring(self
):
423 Disable the map from the base ring into the algebra.
425 Performing a nonsense conversion like this automatically
426 is counterpedagogical. The fallback is to try the usual
427 element constructor, which should also fail.
431 sage: from mjo.eja.eja_algebra import random_eja
435 sage: set_random_seed()
436 sage: J = random_eja()
438 Traceback (most recent call last):
440 ValueError: not an element of this algebra
446 def product_on_basis(self
, i
, j
):
448 Returns the Jordan product of the `i` and `j`th basis elements.
450 This completely defines the Jordan product on the algebra, and
451 is used direclty by our superclass machinery to implement
456 sage: from mjo.eja.eja_algebra import random_eja
460 sage: set_random_seed()
461 sage: J = random_eja()
462 sage: n = J.dimension()
465 sage: bi_bj = J.zero()*J.zero()
467 ....: i = ZZ.random_element(n)
468 ....: j = ZZ.random_element(n)
469 ....: bi = J.monomial(i)
470 ....: bj = J.monomial(j)
471 ....: bi_bj = J.product_on_basis(i,j)
476 # We only stored the lower-triangular portion of the
477 # multiplication table.
479 return self
._multiplication
_table
[i
][j
]
481 return self
._multiplication
_table
[j
][i
]
483 def inner_product(self
, x
, y
):
485 The inner product associated with this Euclidean Jordan algebra.
487 Defaults to the trace inner product, but can be overridden by
488 subclasses if they are sure that the necessary properties are
493 sage: from mjo.eja.eja_algebra import (random_eja,
495 ....: BilinearFormEJA)
499 Our inner product is "associative," which means the following for
500 a symmetric bilinear form::
502 sage: set_random_seed()
503 sage: J = random_eja()
504 sage: x,y,z = J.random_elements(3)
505 sage: (x*y).inner_product(z) == y.inner_product(x*z)
510 Ensure that this is the usual inner product for the algebras
513 sage: set_random_seed()
514 sage: J = HadamardEJA.random_instance()
515 sage: x,y = J.random_elements(2)
516 sage: actual = x.inner_product(y)
517 sage: expected = x.to_vector().inner_product(y.to_vector())
518 sage: actual == expected
521 Ensure that this is one-half of the trace inner-product in a
522 BilinearFormEJA that isn't just the reals (when ``n`` isn't
523 one). This is in Faraut and Koranyi, and also my "On the
526 sage: set_random_seed()
527 sage: J = BilinearFormEJA.random_instance()
528 sage: n = J.dimension()
529 sage: x = J.random_element()
530 sage: y = J.random_element()
531 sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
535 B
= self
._inner
_product
_matrix
536 return (B
*x
.to_vector()).inner_product(y
.to_vector())
539 def is_associative(self
):
541 Return whether or not this algebra's Jordan product is associative.
545 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
549 sage: J = ComplexHermitianEJA(3, field=QQ, orthonormalize=False)
550 sage: J.is_associative()
552 sage: x = sum(J.gens())
553 sage: A = x.subalgebra_generated_by(orthonormalize=False)
554 sage: A.is_associative()
558 return "Associative" in self
.category().axioms()
560 def _is_commutative(self
):
562 Whether or not this algebra's multiplication table is commutative.
564 This method should of course always return ``True``, unless
565 this algebra was constructed with ``check_axioms=False`` and
566 passed an invalid multiplication table.
568 return all( x
*y
== y
*x
for x
in self
.gens() for y
in self
.gens() )
570 def _is_jordanian(self
):
572 Whether or not this algebra's multiplication table respects the
573 Jordan identity `(x^{2})(xy) = x(x^{2}y)`.
575 We only check one arrangement of `x` and `y`, so for a
576 ``True`` result to be truly true, you should also check
577 :meth:`_is_commutative`. This method should of course always
578 return ``True``, unless this algebra was constructed with
579 ``check_axioms=False`` and passed an invalid multiplication table.
581 return all( (self
.monomial(i
)**2)*(self
.monomial(i
)*self
.monomial(j
))
583 (self
.monomial(i
))*((self
.monomial(i
)**2)*self
.monomial(j
))
584 for i
in range(self
.dimension())
585 for j
in range(self
.dimension()) )
587 def _jordan_product_is_associative(self
):
589 Return whether or not this algebra's Jordan product is
590 associative; that is, whether or not `x*(y*z) = (x*y)*z`
593 This method should agree with :meth:`is_associative` unless
594 you lied about the value of the ``associative`` parameter
595 when you constructed the algebra.
599 sage: from mjo.eja.eja_algebra import (random_eja,
600 ....: RealSymmetricEJA,
601 ....: ComplexHermitianEJA,
602 ....: QuaternionHermitianEJA)
606 sage: J = RealSymmetricEJA(4, orthonormalize=False)
607 sage: J._jordan_product_is_associative()
609 sage: x = sum(J.gens())
610 sage: A = x.subalgebra_generated_by()
611 sage: A._jordan_product_is_associative()
616 sage: J = ComplexHermitianEJA(2,field=QQ,orthonormalize=False)
617 sage: J._jordan_product_is_associative()
619 sage: x = sum(J.gens())
620 sage: A = x.subalgebra_generated_by(orthonormalize=False)
621 sage: A._jordan_product_is_associative()
626 sage: J = QuaternionHermitianEJA(2)
627 sage: J._jordan_product_is_associative()
629 sage: x = sum(J.gens())
630 sage: A = x.subalgebra_generated_by()
631 sage: A._jordan_product_is_associative()
636 The values we've presupplied to the constructors agree with
639 sage: set_random_seed()
640 sage: J = random_eja()
641 sage: J.is_associative() == J._jordan_product_is_associative()
647 # Used to check whether or not something is zero.
650 # I don't know of any examples that make this magnitude
651 # necessary because I don't know how to make an
652 # associative algebra when the element subalgebra
653 # construction is unreliable (as it is over RDF; we can't
654 # find the degree of an element because we can't compute
655 # the rank of a matrix). But even multiplication of floats
656 # is non-associative, so *some* epsilon is needed... let's
657 # just take the one from _inner_product_is_associative?
660 for i
in range(self
.dimension()):
661 for j
in range(self
.dimension()):
662 for k
in range(self
.dimension()):
666 diff
= (x
*y
)*z
- x
*(y
*z
)
668 if diff
.norm() > epsilon
:
673 def _inner_product_is_associative(self
):
675 Return whether or not this algebra's inner product `B` is
676 associative; that is, whether or not `B(xy,z) = B(x,yz)`.
678 This method should of course always return ``True``, unless
679 this algebra was constructed with ``check_axioms=False`` and
680 passed an invalid Jordan or inner-product.
684 # Used to check whether or not something is zero.
687 # This choice is sufficient to allow the construction of
688 # QuaternionHermitianEJA(2, field=RDF) with check_axioms=True.
691 for i
in range(self
.dimension()):
692 for j
in range(self
.dimension()):
693 for k
in range(self
.dimension()):
697 diff
= (x
*y
).inner_product(z
) - x
.inner_product(y
*z
)
699 if diff
.abs() > epsilon
:
704 def _element_constructor_(self
, elt
):
706 Construct an element of this algebra or a subalgebra from its
707 EJA element, vector, or matrix representation.
709 This gets called only after the parent element _call_ method
710 fails to find a coercion for the argument.
714 sage: from mjo.eja.eja_algebra import (random_eja,
717 ....: RealSymmetricEJA)
721 The identity in `S^n` is converted to the identity in the EJA::
723 sage: J = RealSymmetricEJA(3)
724 sage: I = matrix.identity(QQ,3)
725 sage: J(I) == J.one()
728 This skew-symmetric matrix can't be represented in the EJA::
730 sage: J = RealSymmetricEJA(3)
731 sage: A = matrix(QQ,3, lambda i,j: i-j)
733 Traceback (most recent call last):
735 ValueError: not an element of this algebra
737 Tuples work as well, provided that the matrix basis for the
738 algebra consists of them::
740 sage: J1 = HadamardEJA(3)
741 sage: J2 = RealSymmetricEJA(2)
742 sage: J = cartesian_product([J1,J2])
743 sage: J( (J1.matrix_basis()[1], J2.matrix_basis()[2]) )
746 Subalgebra elements are embedded into the superalgebra::
748 sage: J = JordanSpinEJA(3)
751 sage: x = sum(J.gens())
752 sage: A = x.subalgebra_generated_by()
758 Ensure that we can convert any element back and forth
759 faithfully between its matrix and algebra representations::
761 sage: set_random_seed()
762 sage: J = random_eja()
763 sage: x = J.random_element()
764 sage: J(x.to_matrix()) == x
767 We cannot coerce elements between algebras just because their
768 matrix representations are compatible::
770 sage: J1 = HadamardEJA(3)
771 sage: J2 = JordanSpinEJA(3)
773 Traceback (most recent call last):
775 ValueError: not an element of this algebra
777 Traceback (most recent call last):
779 ValueError: not an element of this algebra
782 msg
= "not an element of this algebra"
783 if elt
in self
.base_ring():
784 # Ensure that no base ring -> algebra coercion is performed
785 # by this method. There's some stupidity in sage that would
786 # otherwise propagate to this method; for example, sage thinks
787 # that the integer 3 belongs to the space of 2-by-2 matrices.
788 raise ValueError(msg
)
790 if hasattr(elt
, 'superalgebra_element'):
791 # Handle subalgebra elements
792 if elt
.parent().superalgebra() == self
:
793 return elt
.superalgebra_element()
795 if hasattr(elt
, 'sparse_vector'):
796 # Convert a vector into a column-matrix. We check for
797 # "sparse_vector" and not "column" because matrices also
798 # have a "column" method.
801 if elt
not in self
.matrix_space():
802 raise ValueError(msg
)
804 # Thanks for nothing! Matrix spaces aren't vector spaces in
805 # Sage, so we have to figure out its matrix-basis coordinates
806 # ourselves. We use the basis space's ring instead of the
807 # element's ring because the basis space might be an algebraic
808 # closure whereas the base ring of the 3-by-3 identity matrix
809 # could be QQ instead of QQbar.
811 # And, we also have to handle Cartesian product bases (when
812 # the matrix basis consists of tuples) here. The "good news"
813 # is that we're already converting everything to long vectors,
814 # and that strategy works for tuples as well.
816 elt
= self
._matrix
_span
.ambient_vector_space()(_all2list(elt
))
819 coords
= self
._matrix
_span
.coordinate_vector(elt
)
820 except ArithmeticError: # vector is not in free module
821 raise ValueError(msg
)
823 return self
.from_vector(coords
)
827 Return a string representation of ``self``.
831 sage: from mjo.eja.eja_algebra import JordanSpinEJA
835 Ensure that it says what we think it says::
837 sage: JordanSpinEJA(2, field=AA)
838 Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
839 sage: JordanSpinEJA(3, field=RDF)
840 Euclidean Jordan algebra of dimension 3 over Real Double Field
843 fmt
= "Euclidean Jordan algebra of dimension {} over {}"
844 return fmt
.format(self
.dimension(), self
.base_ring())
848 def characteristic_polynomial_of(self
):
850 Return the algebra's "characteristic polynomial of" function,
851 which is itself a multivariate polynomial that, when evaluated
852 at the coordinates of some algebra element, returns that
853 element's characteristic polynomial.
855 The resulting polynomial has `n+1` variables, where `n` is the
856 dimension of this algebra. The first `n` variables correspond to
857 the coordinates of an algebra element: when evaluated at the
858 coordinates of an algebra element with respect to a certain
859 basis, the result is a univariate polynomial (in the one
860 remaining variable ``t``), namely the characteristic polynomial
865 sage: from mjo.eja.eja_algebra import JordanSpinEJA, TrivialEJA
869 The characteristic polynomial in the spin algebra is given in
870 Alizadeh, Example 11.11::
872 sage: J = JordanSpinEJA(3)
873 sage: p = J.characteristic_polynomial_of(); p
874 X0^2 - X1^2 - X2^2 + (-2*t)*X0 + t^2
875 sage: xvec = J.one().to_vector()
879 By definition, the characteristic polynomial is a monic
880 degree-zero polynomial in a rank-zero algebra. Note that
881 Cayley-Hamilton is indeed satisfied since the polynomial
882 ``1`` evaluates to the identity element of the algebra on
885 sage: J = TrivialEJA()
886 sage: J.characteristic_polynomial_of()
893 # The list of coefficient polynomials a_0, a_1, a_2, ..., a_(r-1).
894 a
= self
._charpoly
_coefficients
()
896 # We go to a bit of trouble here to reorder the
897 # indeterminates, so that it's easier to evaluate the
898 # characteristic polynomial at x's coordinates and get back
899 # something in terms of t, which is what we want.
900 S
= PolynomialRing(self
.base_ring(),'t')
904 S
= PolynomialRing(S
, R
.variable_names())
907 return (t
**r
+ sum( a
[k
]*(t
**k
) for k
in range(r
) ))
909 def coordinate_polynomial_ring(self
):
911 The multivariate polynomial ring in which this algebra's
912 :meth:`characteristic_polynomial_of` lives.
916 sage: from mjo.eja.eja_algebra import (HadamardEJA,
917 ....: RealSymmetricEJA)
921 sage: J = HadamardEJA(2)
922 sage: J.coordinate_polynomial_ring()
923 Multivariate Polynomial Ring in X0, X1...
924 sage: J = RealSymmetricEJA(3,field=QQ,orthonormalize=False)
925 sage: J.coordinate_polynomial_ring()
926 Multivariate Polynomial Ring in X0, X1, X2, X3, X4, X5...
929 var_names
= tuple( "X%d" % z
for z
in range(self
.dimension()) )
930 return PolynomialRing(self
.base_ring(), var_names
)
932 def inner_product(self
, x
, y
):
934 The inner product associated with this Euclidean Jordan algebra.
936 Defaults to the trace inner product, but can be overridden by
937 subclasses if they are sure that the necessary properties are
942 sage: from mjo.eja.eja_algebra import (random_eja,
944 ....: BilinearFormEJA)
948 Our inner product is "associative," which means the following for
949 a symmetric bilinear form::
951 sage: set_random_seed()
952 sage: J = random_eja()
953 sage: x,y,z = J.random_elements(3)
954 sage: (x*y).inner_product(z) == y.inner_product(x*z)
959 Ensure that this is the usual inner product for the algebras
962 sage: set_random_seed()
963 sage: J = HadamardEJA.random_instance()
964 sage: x,y = J.random_elements(2)
965 sage: actual = x.inner_product(y)
966 sage: expected = x.to_vector().inner_product(y.to_vector())
967 sage: actual == expected
970 Ensure that this is one-half of the trace inner-product in a
971 BilinearFormEJA that isn't just the reals (when ``n`` isn't
972 one). This is in Faraut and Koranyi, and also my "On the
975 sage: set_random_seed()
976 sage: J = BilinearFormEJA.random_instance()
977 sage: n = J.dimension()
978 sage: x = J.random_element()
979 sage: y = J.random_element()
980 sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
983 B
= self
._inner
_product
_matrix
984 return (B
*x
.to_vector()).inner_product(y
.to_vector())
987 def is_trivial(self
):
989 Return whether or not this algebra is trivial.
991 A trivial algebra contains only the zero element.
995 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
1000 sage: J = ComplexHermitianEJA(3)
1001 sage: J.is_trivial()
1006 sage: J = TrivialEJA()
1007 sage: J.is_trivial()
1011 return self
.dimension() == 0
1014 def multiplication_table(self
):
1016 Return a visual representation of this algebra's multiplication
1017 table (on basis elements).
1021 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1025 sage: J = JordanSpinEJA(4)
1026 sage: J.multiplication_table()
1027 +----++----+----+----+----+
1028 | * || b0 | b1 | b2 | b3 |
1029 +====++====+====+====+====+
1030 | b0 || b0 | b1 | b2 | b3 |
1031 +----++----+----+----+----+
1032 | b1 || b1 | b0 | 0 | 0 |
1033 +----++----+----+----+----+
1034 | b2 || b2 | 0 | b0 | 0 |
1035 +----++----+----+----+----+
1036 | b3 || b3 | 0 | 0 | b0 |
1037 +----++----+----+----+----+
1040 n
= self
.dimension()
1041 # Prepend the header row.
1042 M
= [["*"] + list(self
.gens())]
1044 # And to each subsequent row, prepend an entry that belongs to
1045 # the left-side "header column."
1046 M
+= [ [self
.monomial(i
)] + [ self
.monomial(i
)*self
.monomial(j
)
1050 return table(M
, header_row
=True, header_column
=True, frame
=True)
1053 def matrix_basis(self
):
1055 Return an (often more natural) representation of this algebras
1056 basis as an ordered tuple of matrices.
1058 Every finite-dimensional Euclidean Jordan Algebra is a, up to
1059 Jordan isomorphism, a direct sum of five simple
1060 algebras---four of which comprise Hermitian matrices. And the
1061 last type of algebra can of course be thought of as `n`-by-`1`
1062 column matrices (ambiguusly called column vectors) to avoid
1063 special cases. As a result, matrices (and column vectors) are
1064 a natural representation format for Euclidean Jordan algebra
1067 But, when we construct an algebra from a basis of matrices,
1068 those matrix representations are lost in favor of coordinate
1069 vectors *with respect to* that basis. We could eventually
1070 convert back if we tried hard enough, but having the original
1071 representations handy is valuable enough that we simply store
1072 them and return them from this method.
1074 Why implement this for non-matrix algebras? Avoiding special
1075 cases for the :class:`BilinearFormEJA` pays with simplicity in
1076 its own right. But mainly, we would like to be able to assume
1077 that elements of a :class:`CartesianProductEJA` can be displayed
1078 nicely, without having to have special classes for direct sums
1079 one of whose components was a matrix algebra.
1083 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1084 ....: RealSymmetricEJA)
1088 sage: J = RealSymmetricEJA(2)
1090 Finite family {0: b0, 1: b1, 2: b2}
1091 sage: J.matrix_basis()
1093 [1 0] [ 0 0.7071067811865475?] [0 0]
1094 [0 0], [0.7071067811865475? 0], [0 1]
1099 sage: J = JordanSpinEJA(2)
1101 Finite family {0: b0, 1: b1}
1102 sage: J.matrix_basis()
1108 return self
._matrix
_basis
1111 def matrix_space(self
):
1113 Return the matrix space in which this algebra's elements live, if
1114 we think of them as matrices (including column vectors of the
1117 "By default" this will be an `n`-by-`1` column-matrix space,
1118 except when the algebra is trivial. There it's `n`-by-`n`
1119 (where `n` is zero), to ensure that two elements of the matrix
1120 space (empty matrices) can be multiplied. For algebras of
1121 matrices, this returns the space in which their
1122 real embeddings live.
1126 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
1127 ....: JordanSpinEJA,
1128 ....: QuaternionHermitianEJA,
1133 By default, the matrix representation is just a column-matrix
1134 equivalent to the vector representation::
1136 sage: J = JordanSpinEJA(3)
1137 sage: J.matrix_space()
1138 Full MatrixSpace of 3 by 1 dense matrices over Algebraic
1141 The matrix representation in the trivial algebra is
1142 zero-by-zero instead of the usual `n`-by-one::
1144 sage: J = TrivialEJA()
1145 sage: J.matrix_space()
1146 Full MatrixSpace of 0 by 0 dense matrices over Algebraic
1149 The matrix space for complex/quaternion Hermitian matrix EJA
1150 is the space in which their real-embeddings live, not the
1151 original complex/quaternion matrix space::
1153 sage: J = ComplexHermitianEJA(2,field=QQ,orthonormalize=False)
1154 sage: J.matrix_space()
1155 Module of 2 by 2 matrices with entries in Algebraic Field over
1156 the scalar ring Rational Field
1157 sage: J = QuaternionHermitianEJA(1,field=QQ,orthonormalize=False)
1158 sage: J.matrix_space()
1159 Module of 1 by 1 matrices with entries in Quaternion
1160 Algebra (-1, -1) with base ring Rational Field over
1161 the scalar ring Rational Field
1164 return self
._matrix
_space
1170 Return the unit element of this algebra.
1174 sage: from mjo.eja.eja_algebra import (HadamardEJA,
1179 We can compute unit element in the Hadamard EJA::
1181 sage: J = HadamardEJA(5)
1183 b0 + b1 + b2 + b3 + b4
1185 The unit element in the Hadamard EJA is inherited in the
1186 subalgebras generated by its elements::
1188 sage: J = HadamardEJA(5)
1190 b0 + b1 + b2 + b3 + b4
1191 sage: x = sum(J.gens())
1192 sage: A = x.subalgebra_generated_by(orthonormalize=False)
1195 sage: A.one().superalgebra_element()
1196 b0 + b1 + b2 + b3 + b4
1200 The identity element acts like the identity, regardless of
1201 whether or not we orthonormalize::
1203 sage: set_random_seed()
1204 sage: J = random_eja()
1205 sage: x = J.random_element()
1206 sage: J.one()*x == x and x*J.one() == x
1208 sage: A = x.subalgebra_generated_by()
1209 sage: y = A.random_element()
1210 sage: A.one()*y == y and y*A.one() == y
1215 sage: set_random_seed()
1216 sage: J = random_eja(field=QQ, orthonormalize=False)
1217 sage: x = J.random_element()
1218 sage: J.one()*x == x and x*J.one() == x
1220 sage: A = x.subalgebra_generated_by(orthonormalize=False)
1221 sage: y = A.random_element()
1222 sage: A.one()*y == y and y*A.one() == y
1225 The matrix of the unit element's operator is the identity,
1226 regardless of the base field and whether or not we
1229 sage: set_random_seed()
1230 sage: J = random_eja()
1231 sage: actual = J.one().operator().matrix()
1232 sage: expected = matrix.identity(J.base_ring(), J.dimension())
1233 sage: actual == expected
1235 sage: x = J.random_element()
1236 sage: A = x.subalgebra_generated_by()
1237 sage: actual = A.one().operator().matrix()
1238 sage: expected = matrix.identity(A.base_ring(), A.dimension())
1239 sage: actual == expected
1244 sage: set_random_seed()
1245 sage: J = random_eja(field=QQ, orthonormalize=False)
1246 sage: actual = J.one().operator().matrix()
1247 sage: expected = matrix.identity(J.base_ring(), J.dimension())
1248 sage: actual == expected
1250 sage: x = J.random_element()
1251 sage: A = x.subalgebra_generated_by(orthonormalize=False)
1252 sage: actual = A.one().operator().matrix()
1253 sage: expected = matrix.identity(A.base_ring(), A.dimension())
1254 sage: actual == expected
1257 Ensure that the cached unit element (often precomputed by
1258 hand) agrees with the computed one::
1260 sage: set_random_seed()
1261 sage: J = random_eja()
1262 sage: cached = J.one()
1263 sage: J.one.clear_cache()
1264 sage: J.one() == cached
1269 sage: set_random_seed()
1270 sage: J = random_eja(field=QQ, orthonormalize=False)
1271 sage: cached = J.one()
1272 sage: J.one.clear_cache()
1273 sage: J.one() == cached
1277 # We can brute-force compute the matrices of the operators
1278 # that correspond to the basis elements of this algebra.
1279 # If some linear combination of those basis elements is the
1280 # algebra identity, then the same linear combination of
1281 # their matrices has to be the identity matrix.
1283 # Of course, matrices aren't vectors in sage, so we have to
1284 # appeal to the "long vectors" isometry.
1286 V
= VectorSpace(self
.base_ring(), self
.dimension()**2)
1287 oper_vecs
= [ V(g
.operator().matrix().list()) for g
in self
.gens() ]
1289 # Now we use basic linear algebra to find the coefficients,
1290 # of the matrices-as-vectors-linear-combination, which should
1291 # work for the original algebra basis too.
1292 A
= matrix(self
.base_ring(), oper_vecs
)
1294 # We used the isometry on the left-hand side already, but we
1295 # still need to do it for the right-hand side. Recall that we
1296 # wanted something that summed to the identity matrix.
1297 b
= V( matrix
.identity(self
.base_ring(), self
.dimension()).list() )
1299 # Now if there's an identity element in the algebra, this
1300 # should work. We solve on the left to avoid having to
1301 # transpose the matrix "A".
1302 return self
.from_vector(A
.solve_left(b
))
1305 def peirce_decomposition(self
, c
):
1307 The Peirce decomposition of this algebra relative to the
1310 In the future, this can be extended to a complete system of
1311 orthogonal idempotents.
1315 - ``c`` -- an idempotent of this algebra.
1319 A triple (J0, J5, J1) containing two subalgebras and one subspace
1322 - ``J0`` -- the algebra on the eigenspace of ``c.operator()``
1323 corresponding to the eigenvalue zero.
1325 - ``J5`` -- the eigenspace (NOT a subalgebra) of ``c.operator()``
1326 corresponding to the eigenvalue one-half.
1328 - ``J1`` -- the algebra on the eigenspace of ``c.operator()``
1329 corresponding to the eigenvalue one.
1331 These are the only possible eigenspaces for that operator, and this
1332 algebra is a direct sum of them. The spaces ``J0`` and ``J1`` are
1333 orthogonal, and are subalgebras of this algebra with the appropriate
1338 sage: from mjo.eja.eja_algebra import random_eja, RealSymmetricEJA
1342 The canonical example comes from the symmetric matrices, which
1343 decompose into diagonal and off-diagonal parts::
1345 sage: J = RealSymmetricEJA(3)
1346 sage: C = matrix(QQ, [ [1,0,0],
1350 sage: J0,J5,J1 = J.peirce_decomposition(c)
1352 Euclidean Jordan algebra of dimension 1...
1354 Vector space of degree 6 and dimension 2...
1356 Euclidean Jordan algebra of dimension 3...
1357 sage: J0.one().to_matrix()
1361 sage: orig_df = AA.options.display_format
1362 sage: AA.options.display_format = 'radical'
1363 sage: J.from_vector(J5.basis()[0]).to_matrix()
1367 sage: J.from_vector(J5.basis()[1]).to_matrix()
1371 sage: AA.options.display_format = orig_df
1372 sage: J1.one().to_matrix()
1379 Every algebra decomposes trivially with respect to its identity
1382 sage: set_random_seed()
1383 sage: J = random_eja()
1384 sage: J0,J5,J1 = J.peirce_decomposition(J.one())
1385 sage: J0.dimension() == 0 and J5.dimension() == 0
1387 sage: J1.superalgebra() == J and J1.dimension() == J.dimension()
1390 The decomposition is into eigenspaces, and its components are
1391 therefore necessarily orthogonal. Moreover, the identity
1392 elements in the two subalgebras are the projections onto their
1393 respective subspaces of the superalgebra's identity element::
1395 sage: set_random_seed()
1396 sage: J = random_eja()
1397 sage: x = J.random_element()
1398 sage: if not J.is_trivial():
1399 ....: while x.is_nilpotent():
1400 ....: x = J.random_element()
1401 sage: c = x.subalgebra_idempotent()
1402 sage: J0,J5,J1 = J.peirce_decomposition(c)
1404 sage: for (w,y,z) in zip(J0.basis(), J5.basis(), J1.basis()):
1405 ....: w = w.superalgebra_element()
1406 ....: y = J.from_vector(y)
1407 ....: z = z.superalgebra_element()
1408 ....: ipsum += w.inner_product(y).abs()
1409 ....: ipsum += w.inner_product(z).abs()
1410 ....: ipsum += y.inner_product(z).abs()
1413 sage: J1(c) == J1.one()
1415 sage: J0(J.one() - c) == J0.one()
1419 if not c
.is_idempotent():
1420 raise ValueError("element is not idempotent: %s" % c
)
1422 # Default these to what they should be if they turn out to be
1423 # trivial, because eigenspaces_left() won't return eigenvalues
1424 # corresponding to trivial spaces (e.g. it returns only the
1425 # eigenspace corresponding to lambda=1 if you take the
1426 # decomposition relative to the identity element).
1427 trivial
= self
.subalgebra((), check_axioms
=False)
1428 J0
= trivial
# eigenvalue zero
1429 J5
= VectorSpace(self
.base_ring(), 0) # eigenvalue one-half
1430 J1
= trivial
# eigenvalue one
1432 for (eigval
, eigspace
) in c
.operator().matrix().right_eigenspaces():
1433 if eigval
== ~
(self
.base_ring()(2)):
1436 gens
= tuple( self
.from_vector(b
) for b
in eigspace
.basis() )
1437 subalg
= self
.subalgebra(gens
, check_axioms
=False)
1443 raise ValueError("unexpected eigenvalue: %s" % eigval
)
1448 def random_element(self
, thorough
=False):
1450 Return a random element of this algebra.
1452 Our algebra superclass method only returns a linear
1453 combination of at most two basis elements. We instead
1454 want the vector space "random element" method that
1455 returns a more diverse selection.
1459 - ``thorough`` -- (boolean; default False) whether or not we
1460 should generate irrational coefficients for the random
1461 element when our base ring is irrational; this slows the
1462 algebra operations to a crawl, but any truly random method
1466 # For a general base ring... maybe we can trust this to do the
1467 # right thing? Unlikely, but.
1468 V
= self
.vector_space()
1469 v
= V
.random_element()
1471 if self
.base_ring() is AA
:
1472 # The "random element" method of the algebraic reals is
1473 # stupid at the moment, and only returns integers between
1474 # -2 and 2, inclusive:
1476 # https://trac.sagemath.org/ticket/30875
1478 # Instead, we implement our own "random vector" method,
1479 # and then coerce that into the algebra. We use the vector
1480 # space degree here instead of the dimension because a
1481 # subalgebra could (for example) be spanned by only two
1482 # vectors, each with five coordinates. We need to
1483 # generate all five coordinates.
1485 v
*= QQbar
.random_element().real()
1487 v
*= QQ
.random_element()
1489 return self
.from_vector(V
.coordinate_vector(v
))
1491 def random_elements(self
, count
, thorough
=False):
1493 Return ``count`` random elements as a tuple.
1497 - ``thorough`` -- (boolean; default False) whether or not we
1498 should generate irrational coefficients for the random
1499 elements when our base ring is irrational; this slows the
1500 algebra operations to a crawl, but any truly random method
1505 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1509 sage: J = JordanSpinEJA(3)
1510 sage: x,y,z = J.random_elements(3)
1511 sage: all( [ x in J, y in J, z in J ])
1513 sage: len( J.random_elements(10) ) == 10
1517 return tuple( self
.random_element(thorough
)
1518 for idx
in range(count
) )
1521 def operator_polynomial_matrix(self
):
1523 Return the matrix of polynomials (over this algebra's
1524 :meth:`coordinate_polynomial_ring`) that, when evaluated at
1525 the basis coordinates of an element `x`, produces the basis
1526 representation of `L_{x}`.
1530 sage: from mjo.eja.eja_algebra import (HadamardEJA,
1531 ....: JordanSpinEJA)
1535 sage: J = HadamardEJA(4)
1536 sage: L_x = J.operator_polynomial_matrix()
1543 sage: d = zip(J.coordinate_polynomial_ring().gens(), x.to_vector())
1544 sage: L_x.subs(dict(d))
1552 sage: J = JordanSpinEJA(4)
1553 sage: L_x = J.operator_polynomial_matrix()
1560 sage: d = zip(J.coordinate_polynomial_ring().gens(), x.to_vector())
1561 sage: L_x.subs(dict(d))
1568 R
= self
.coordinate_polynomial_ring()
1571 # From a result in my book, these are the entries of the
1572 # basis representation of L_x.
1573 return sum( v
*self
.monomial(k
).operator().matrix()[i
,j
]
1574 for (k
,v
) in enumerate(R
.gens()) )
1576 n
= self
.dimension()
1577 return matrix(R
, n
, n
, L_x_i_j
)
1580 def _charpoly_coefficients(self
):
1582 The `r` polynomial coefficients of the "characteristic polynomial
1587 sage: from mjo.eja.eja_algebra import random_eja
1591 The theory shows that these are all homogeneous polynomials of
1594 sage: set_random_seed()
1595 sage: J = random_eja()
1596 sage: all(p.is_homogeneous() for p in J._charpoly_coefficients())
1600 n
= self
.dimension()
1601 R
= self
.coordinate_polynomial_ring()
1602 F
= R
.fraction_field()
1604 L_x
= self
.operator_polynomial_matrix()
1607 if self
.rank
.is_in_cache():
1609 # There's no need to pad the system with redundant
1610 # columns if we *know* they'll be redundant.
1613 # Compute an extra power in case the rank is equal to
1614 # the dimension (otherwise, we would stop at x^(r-1)).
1615 x_powers
= [ (L_x
**k
)*self
.one().to_vector()
1616 for k
in range(n
+1) ]
1617 A
= matrix
.column(F
, x_powers
[:n
])
1618 AE
= A
.extended_echelon_form()
1625 # The theory says that only the first "r" coefficients are
1626 # nonzero, and they actually live in the original polynomial
1627 # ring and not the fraction field. We negate them because in
1628 # the actual characteristic polynomial, they get moved to the
1629 # other side where x^r lives. We don't bother to trim A_rref
1630 # down to a square matrix and solve the resulting system,
1631 # because the upper-left r-by-r portion of A_rref is
1632 # guaranteed to be the identity matrix, so e.g.
1634 # A_rref.solve_right(Y)
1636 # would just be returning Y.
1637 return (-E
*b
)[:r
].change_ring(R
)
1642 Return the rank of this EJA.
1644 This is a cached method because we know the rank a priori for
1645 all of the algebras we can construct. Thus we can avoid the
1646 expensive ``_charpoly_coefficients()`` call unless we truly
1647 need to compute the whole characteristic polynomial.
1651 sage: from mjo.eja.eja_algebra import (HadamardEJA,
1652 ....: JordanSpinEJA,
1653 ....: RealSymmetricEJA,
1654 ....: ComplexHermitianEJA,
1655 ....: QuaternionHermitianEJA,
1660 The rank of the Jordan spin algebra is always two::
1662 sage: JordanSpinEJA(2).rank()
1664 sage: JordanSpinEJA(3).rank()
1666 sage: JordanSpinEJA(4).rank()
1669 The rank of the `n`-by-`n` Hermitian real, complex, or
1670 quaternion matrices is `n`::
1672 sage: RealSymmetricEJA(4).rank()
1674 sage: ComplexHermitianEJA(3).rank()
1676 sage: QuaternionHermitianEJA(2).rank()
1681 Ensure that every EJA that we know how to construct has a
1682 positive integer rank, unless the algebra is trivial in
1683 which case its rank will be zero::
1685 sage: set_random_seed()
1686 sage: J = random_eja()
1690 sage: r > 0 or (r == 0 and J.is_trivial())
1693 Ensure that computing the rank actually works, since the ranks
1694 of all simple algebras are known and will be cached by default::
1696 sage: set_random_seed() # long time
1697 sage: J = random_eja() # long time
1698 sage: cached = J.rank() # long time
1699 sage: J.rank.clear_cache() # long time
1700 sage: J.rank() == cached # long time
1704 return len(self
._charpoly
_coefficients
())
1707 def subalgebra(self
, basis
, **kwargs
):
1709 Create a subalgebra of this algebra from the given basis.
1711 from mjo
.eja
.eja_subalgebra
import FiniteDimensionalEJASubalgebra
1712 return FiniteDimensionalEJASubalgebra(self
, basis
, **kwargs
)
1715 def vector_space(self
):
1717 Return the vector space that underlies this algebra.
1721 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1725 sage: J = RealSymmetricEJA(2)
1726 sage: J.vector_space()
1727 Vector space of dimension 3 over...
1730 return self
.zero().to_vector().parent().ambient_vector_space()
1734 class RationalBasisEJA(FiniteDimensionalEJA
):
1736 Algebras whose supplied basis elements have all rational entries.
1740 sage: from mjo.eja.eja_algebra import BilinearFormEJA
1744 The supplied basis is orthonormalized by default::
1746 sage: B = matrix(QQ, [[1, 0, 0], [0, 25, -32], [0, -32, 41]])
1747 sage: J = BilinearFormEJA(B)
1748 sage: J.matrix_basis()
1765 # Abuse the check_field parameter to check that the entries of
1766 # out basis (in ambient coordinates) are in the field QQ.
1767 # Use _all2list to get the vector coordinates of octonion
1768 # entries and not the octonions themselves (which are not
1770 if not all( all(b_i
in QQ
for b_i
in _all2list(b
))
1772 raise TypeError("basis not rational")
1774 super().__init
__(basis
,
1778 check_field
=check_field
,
1781 self
._rational
_algebra
= None
1783 # There's no point in constructing the extra algebra if this
1784 # one is already rational.
1786 # Note: the same Jordan and inner-products work here,
1787 # because they are necessarily defined with respect to
1788 # ambient coordinates and not any particular basis.
1789 self
._rational
_algebra
= FiniteDimensionalEJA(
1794 matrix_space
=self
.matrix_space(),
1795 associative
=self
.is_associative(),
1796 orthonormalize
=False,
1800 def rational_algebra(self
):
1801 # Using None as a flag here (rather than just assigning "self"
1802 # to self._rational_algebra by default) feels a little bit
1803 # more sane to me in a garbage-collected environment.
1804 if self
._rational
_algebra
is None:
1807 return self
._rational
_algebra
1810 def _charpoly_coefficients(self
):
1814 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
1815 ....: JordanSpinEJA)
1819 The base ring of the resulting polynomial coefficients is what
1820 it should be, and not the rationals (unless the algebra was
1821 already over the rationals)::
1823 sage: J = JordanSpinEJA(3)
1824 sage: J._charpoly_coefficients()
1825 (X0^2 - X1^2 - X2^2, -2*X0)
1826 sage: a0 = J._charpoly_coefficients()[0]
1828 Algebraic Real Field
1829 sage: a0.base_ring()
1830 Algebraic Real Field
1833 if self
.rational_algebra() is self
:
1834 # Bypass the hijinks if they won't benefit us.
1835 return super()._charpoly
_coefficients
()
1837 # Do the computation over the rationals. The answer will be
1838 # the same, because all we've done is a change of basis.
1839 # Then, change back from QQ to our real base ring
1840 a
= ( a_i
.change_ring(self
.base_ring())
1841 for a_i
in self
.rational_algebra()._charpoly
_coefficients
() )
1843 # Otherwise, convert the coordinate variables back to the
1844 # deorthonormalized ones.
1845 R
= self
.coordinate_polynomial_ring()
1846 from sage
.modules
.free_module_element
import vector
1847 X
= vector(R
, R
.gens())
1848 BX
= self
._deortho
_matrix
*X
1850 subs_dict
= { X[i]: BX[i] for i in range(len(X)) }
1851 return tuple( a_i
.subs(subs_dict
) for a_i
in a
)
1853 class ConcreteEJA(FiniteDimensionalEJA
):
1855 A class for the Euclidean Jordan algebras that we know by name.
1857 These are the Jordan algebras whose basis, multiplication table,
1858 rank, and so on are known a priori. More to the point, they are
1859 the Euclidean Jordan algebras for which we are able to conjure up
1860 a "random instance."
1864 sage: from mjo.eja.eja_algebra import ConcreteEJA
1868 Our basis is normalized with respect to the algebra's inner
1869 product, unless we specify otherwise::
1871 sage: set_random_seed()
1872 sage: J = ConcreteEJA.random_instance()
1873 sage: all( b.norm() == 1 for b in J.gens() )
1876 Since our basis is orthonormal with respect to the algebra's inner
1877 product, and since we know that this algebra is an EJA, any
1878 left-multiplication operator's matrix will be symmetric because
1879 natural->EJA basis representation is an isometry and within the
1880 EJA the operator is self-adjoint by the Jordan axiom::
1882 sage: set_random_seed()
1883 sage: J = ConcreteEJA.random_instance()
1884 sage: x = J.random_element()
1885 sage: x.operator().is_self_adjoint()
1890 def _max_random_instance_dimension():
1892 The maximum dimension of any random instance. Ten dimensions seems
1893 to be about the point where everything takes a turn for the
1894 worse. And dimension ten (but not nine) allows the 4-by-4 real
1895 Hermitian matrices, the 2-by-2 quaternion Hermitian matrices,
1896 and the 2-by-2 octonion Hermitian matrices.
1901 def _max_random_instance_size(max_dimension
):
1903 Return an integer "size" that is an upper bound on the size of
1904 this algebra when it is used in a random test case. This size
1905 (which can be passed to the algebra's constructor) is itself
1906 based on the ``max_dimension`` parameter.
1908 This method must be implemented in each subclass.
1910 raise NotImplementedError
1913 def random_instance(cls
, max_dimension
=None, *args
, **kwargs
):
1915 Return a random instance of this type of algebra whose dimension
1916 is less than or equal to the lesser of ``max_dimension`` and
1917 the value returned by ``_max_random_instance_dimension()``. If
1918 the dimension bound is omitted, then only the
1919 ``_max_random_instance_dimension()`` is used as a bound.
1921 This method should be implemented in each subclass.
1925 sage: from mjo.eja.eja_algebra import ConcreteEJA
1929 Both the class bound and the ``max_dimension`` argument are upper
1930 bounds on the dimension of the algebra returned::
1932 sage: from sage.misc.prandom import choice
1933 sage: eja_class = choice(ConcreteEJA.__subclasses__())
1934 sage: class_max_d = eja_class._max_random_instance_dimension()
1935 sage: J = eja_class.random_instance(max_dimension=20,
1937 ....: orthonormalize=False)
1938 sage: J.dimension() <= class_max_d
1940 sage: J = eja_class.random_instance(max_dimension=2,
1942 ....: orthonormalize=False)
1943 sage: J.dimension() <= 2
1947 from sage
.misc
.prandom
import choice
1948 eja_class
= choice(cls
.__subclasses
__())
1950 # These all bubble up to the RationalBasisEJA superclass
1951 # constructor, so any (kw)args valid there are also valid
1953 return eja_class
.random_instance(max_dimension
, *args
, **kwargs
)
1956 class MatrixEJA(FiniteDimensionalEJA
):
1958 def _denormalized_basis(A
):
1960 Returns a basis for the space of complex Hermitian n-by-n matrices.
1962 Why do we embed these? Basically, because all of numerical linear
1963 algebra assumes that you're working with vectors consisting of `n`
1964 entries from a field and scalars from the same field. There's no way
1965 to tell SageMath that (for example) the vectors contain complex
1966 numbers, while the scalar field is real.
1970 sage: from mjo.hurwitz import (ComplexMatrixAlgebra,
1971 ....: QuaternionMatrixAlgebra,
1972 ....: OctonionMatrixAlgebra)
1973 sage: from mjo.eja.eja_algebra import MatrixEJA
1977 sage: set_random_seed()
1978 sage: n = ZZ.random_element(1,5)
1979 sage: A = MatrixSpace(QQ, n)
1980 sage: B = MatrixEJA._denormalized_basis(A)
1981 sage: all( M.is_hermitian() for M in B)
1986 sage: set_random_seed()
1987 sage: n = ZZ.random_element(1,5)
1988 sage: A = ComplexMatrixAlgebra(n, scalars=QQ)
1989 sage: B = MatrixEJA._denormalized_basis(A)
1990 sage: all( M.is_hermitian() for M in B)
1995 sage: set_random_seed()
1996 sage: n = ZZ.random_element(1,5)
1997 sage: A = QuaternionMatrixAlgebra(n, scalars=QQ)
1998 sage: B = MatrixEJA._denormalized_basis(A)
1999 sage: all( M.is_hermitian() for M in B )
2004 sage: set_random_seed()
2005 sage: n = ZZ.random_element(1,5)
2006 sage: A = OctonionMatrixAlgebra(n, scalars=QQ)
2007 sage: B = MatrixEJA._denormalized_basis(A)
2008 sage: all( M.is_hermitian() for M in B )
2012 # These work for real MatrixSpace, whose monomials only have
2013 # two coordinates (because the last one would always be "1").
2014 es
= A
.base_ring().gens()
2015 gen
= lambda A
,m
: A
.monomial(m
[:2])
2017 if hasattr(A
, 'entry_algebra_gens'):
2018 # We've got a MatrixAlgebra, and its monomials will have
2019 # three coordinates.
2020 es
= A
.entry_algebra_gens()
2021 gen
= lambda A
,m
: A
.monomial(m
)
2024 for i
in range(A
.nrows()):
2025 for j
in range(i
+1):
2027 E_ii
= gen(A
, (i
,j
,es
[0]))
2031 E_ij
= gen(A
, (i
,j
,e
))
2032 E_ij
+= E_ij
.conjugate_transpose()
2035 return tuple( basis
)
2038 def jordan_product(X
,Y
):
2039 return (X
*Y
+ Y
*X
)/2
2042 def trace_inner_product(X
,Y
):
2044 A trace inner-product for matrices that aren't embedded in the
2045 reals. It takes MATRICES as arguments, not EJA elements.
2049 sage: from mjo.eja.eja_algebra import (RealSymmetricEJA,
2050 ....: ComplexHermitianEJA,
2051 ....: QuaternionHermitianEJA,
2052 ....: OctonionHermitianEJA)
2056 sage: J = RealSymmetricEJA(2,field=QQ,orthonormalize=False)
2057 sage: I = J.one().to_matrix()
2058 sage: J.trace_inner_product(I, -I)
2063 sage: J = ComplexHermitianEJA(2,field=QQ,orthonormalize=False)
2064 sage: I = J.one().to_matrix()
2065 sage: J.trace_inner_product(I, -I)
2070 sage: J = QuaternionHermitianEJA(2,field=QQ,orthonormalize=False)
2071 sage: I = J.one().to_matrix()
2072 sage: J.trace_inner_product(I, -I)
2077 sage: J = OctonionHermitianEJA(2,field=QQ,orthonormalize=False)
2078 sage: I = J.one().to_matrix()
2079 sage: J.trace_inner_product(I, -I)
2084 if hasattr(tr
, 'coefficient'):
2085 # Works for octonions, and has to come first because they
2086 # also have a "real()" method that doesn't return an
2087 # element of the scalar ring.
2088 return tr
.coefficient(0)
2089 elif hasattr(tr
, 'coefficient_tuple'):
2090 # Works for quaternions.
2091 return tr
.coefficient_tuple()[0]
2093 # Works for real and complex numbers.
2097 def __init__(self
, matrix_space
, **kwargs
):
2098 # We know this is a valid EJA, but will double-check
2099 # if the user passes check_axioms=True.
2100 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2102 super().__init
__(self
._denormalized
_basis
(matrix_space
),
2103 self
.jordan_product
,
2104 self
.trace_inner_product
,
2105 field
=matrix_space
.base_ring(),
2106 matrix_space
=matrix_space
,
2109 self
.rank
.set_cache(matrix_space
.nrows())
2110 self
.one
.set_cache( self(matrix_space
.one()) )
2112 class RealSymmetricEJA(MatrixEJA
, RationalBasisEJA
, ConcreteEJA
):
2114 The rank-n simple EJA consisting of real symmetric n-by-n
2115 matrices, the usual symmetric Jordan product, and the trace inner
2116 product. It has dimension `(n^2 + n)/2` over the reals.
2120 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
2124 sage: J = RealSymmetricEJA(2)
2125 sage: b0, b1, b2 = J.gens()
2133 In theory, our "field" can be any subfield of the reals::
2135 sage: RealSymmetricEJA(2, field=RDF, check_axioms=True)
2136 Euclidean Jordan algebra of dimension 3 over Real Double Field
2137 sage: RealSymmetricEJA(2, field=RR, check_axioms=True)
2138 Euclidean Jordan algebra of dimension 3 over Real Field with
2139 53 bits of precision
2143 The dimension of this algebra is `(n^2 + n) / 2`::
2145 sage: set_random_seed()
2146 sage: d = RealSymmetricEJA._max_random_instance_dimension()
2147 sage: n = RealSymmetricEJA._max_random_instance_size(d)
2148 sage: J = RealSymmetricEJA(n)
2149 sage: J.dimension() == (n^2 + n)/2
2152 The Jordan multiplication is what we think it is::
2154 sage: set_random_seed()
2155 sage: J = RealSymmetricEJA.random_instance()
2156 sage: x,y = J.random_elements(2)
2157 sage: actual = (x*y).to_matrix()
2158 sage: X = x.to_matrix()
2159 sage: Y = y.to_matrix()
2160 sage: expected = (X*Y + Y*X)/2
2161 sage: actual == expected
2163 sage: J(expected) == x*y
2166 We can change the generator prefix::
2168 sage: RealSymmetricEJA(3, prefix='q').gens()
2169 (q0, q1, q2, q3, q4, q5)
2171 We can construct the (trivial) algebra of rank zero::
2173 sage: RealSymmetricEJA(0)
2174 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2178 def _max_random_instance_size(max_dimension
):
2179 # Obtained by solving d = (n^2 + n)/2.
2180 # The ZZ-int-ZZ thing is just "floor."
2181 return ZZ(int(ZZ(8*max_dimension
+ 1).sqrt()/2 - 1/2))
2184 def random_instance(cls
, max_dimension
=None, *args
, **kwargs
):
2186 Return a random instance of this type of algebra.
2188 class_max_d
= cls
._max
_random
_instance
_dimension
()
2189 if (max_dimension
is None or max_dimension
> class_max_d
):
2190 max_dimension
= class_max_d
2191 max_size
= cls
._max
_random
_instance
_size
(max_dimension
)
2192 n
= ZZ
.random_element(max_size
+ 1)
2193 return cls(n
, **kwargs
)
2195 def __init__(self
, n
, field
=AA
, **kwargs
):
2196 A
= MatrixSpace(field
, n
)
2197 super().__init
__(A
, **kwargs
)
2199 from mjo
.eja
.eja_cache
import real_symmetric_eja_coeffs
2200 a
= real_symmetric_eja_coeffs(self
)
2202 self
.rational_algebra()._charpoly
_coefficients
.set_cache(a
)
2206 class ComplexHermitianEJA(MatrixEJA
, RationalBasisEJA
, ConcreteEJA
):
2208 The rank-n simple EJA consisting of complex Hermitian n-by-n
2209 matrices over the real numbers, the usual symmetric Jordan product,
2210 and the real-part-of-trace inner product. It has dimension `n^2` over
2215 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
2219 In theory, our "field" can be any subfield of the reals, but we
2220 can't use inexact real fields at the moment because SageMath
2221 doesn't know how to convert their elements into complex numbers,
2222 or even into algebraic reals::
2225 Traceback (most recent call last):
2227 TypeError: Illegal initializer for algebraic number
2229 Traceback (most recent call last):
2231 TypeError: Illegal initializer for algebraic number
2233 This causes the following error when we try to scale a matrix of
2234 complex numbers by an inexact real number::
2236 sage: ComplexHermitianEJA(2,field=RR)
2237 Traceback (most recent call last):
2239 TypeError: Unable to coerce entries (=(1.00000000000000,
2240 -0.000000000000000)) to coefficients in Algebraic Real Field
2244 The dimension of this algebra is `n^2`::
2246 sage: set_random_seed()
2247 sage: d = ComplexHermitianEJA._max_random_instance_dimension()
2248 sage: n = ComplexHermitianEJA._max_random_instance_size(d)
2249 sage: J = ComplexHermitianEJA(n)
2250 sage: J.dimension() == n^2
2253 The Jordan multiplication is what we think it is::
2255 sage: set_random_seed()
2256 sage: J = ComplexHermitianEJA.random_instance()
2257 sage: x,y = J.random_elements(2)
2258 sage: actual = (x*y).to_matrix()
2259 sage: X = x.to_matrix()
2260 sage: Y = y.to_matrix()
2261 sage: expected = (X*Y + Y*X)/2
2262 sage: actual == expected
2264 sage: J(expected) == x*y
2267 We can change the generator prefix::
2269 sage: ComplexHermitianEJA(2, prefix='z').gens()
2272 We can construct the (trivial) algebra of rank zero::
2274 sage: ComplexHermitianEJA(0)
2275 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2278 def __init__(self
, n
, field
=AA
, **kwargs
):
2279 from mjo
.hurwitz
import ComplexMatrixAlgebra
2280 A
= ComplexMatrixAlgebra(n
, scalars
=field
)
2281 super().__init
__(A
, **kwargs
)
2283 from mjo
.eja
.eja_cache
import complex_hermitian_eja_coeffs
2284 a
= complex_hermitian_eja_coeffs(self
)
2286 self
.rational_algebra()._charpoly
_coefficients
.set_cache(a
)
2289 def _max_random_instance_size(max_dimension
):
2290 # Obtained by solving d = n^2.
2291 # The ZZ-int-ZZ thing is just "floor."
2292 return ZZ(int(ZZ(max_dimension
).sqrt()))
2295 def random_instance(cls
, max_dimension
=None, *args
, **kwargs
):
2297 Return a random instance of this type of algebra.
2299 class_max_d
= cls
._max
_random
_instance
_dimension
()
2300 if (max_dimension
is None or max_dimension
> class_max_d
):
2301 max_dimension
= class_max_d
2302 max_size
= cls
._max
_random
_instance
_size
(max_dimension
)
2303 n
= ZZ
.random_element(max_size
+ 1)
2304 return cls(n
, **kwargs
)
2307 class QuaternionHermitianEJA(MatrixEJA
, RationalBasisEJA
, ConcreteEJA
):
2309 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
2310 matrices, the usual symmetric Jordan product, and the
2311 real-part-of-trace inner product. It has dimension `2n^2 - n` over
2316 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2320 In theory, our "field" can be any subfield of the reals::
2322 sage: QuaternionHermitianEJA(2, field=RDF, check_axioms=True)
2323 Euclidean Jordan algebra of dimension 6 over Real Double Field
2324 sage: QuaternionHermitianEJA(2, field=RR, check_axioms=True)
2325 Euclidean Jordan algebra of dimension 6 over Real Field with
2326 53 bits of precision
2330 The dimension of this algebra is `2*n^2 - n`::
2332 sage: set_random_seed()
2333 sage: d = QuaternionHermitianEJA._max_random_instance_dimension()
2334 sage: n = QuaternionHermitianEJA._max_random_instance_size(d)
2335 sage: J = QuaternionHermitianEJA(n)
2336 sage: J.dimension() == 2*(n^2) - n
2339 The Jordan multiplication is what we think it is::
2341 sage: set_random_seed()
2342 sage: J = QuaternionHermitianEJA.random_instance()
2343 sage: x,y = J.random_elements(2)
2344 sage: actual = (x*y).to_matrix()
2345 sage: X = x.to_matrix()
2346 sage: Y = y.to_matrix()
2347 sage: expected = (X*Y + Y*X)/2
2348 sage: actual == expected
2350 sage: J(expected) == x*y
2353 We can change the generator prefix::
2355 sage: QuaternionHermitianEJA(2, prefix='a').gens()
2356 (a0, a1, a2, a3, a4, a5)
2358 We can construct the (trivial) algebra of rank zero::
2360 sage: QuaternionHermitianEJA(0)
2361 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2364 def __init__(self
, n
, field
=AA
, **kwargs
):
2365 from mjo
.hurwitz
import QuaternionMatrixAlgebra
2366 A
= QuaternionMatrixAlgebra(n
, scalars
=field
)
2367 super().__init
__(A
, **kwargs
)
2369 from mjo
.eja
.eja_cache
import quaternion_hermitian_eja_coeffs
2370 a
= quaternion_hermitian_eja_coeffs(self
)
2372 self
.rational_algebra()._charpoly
_coefficients
.set_cache(a
)
2377 def _max_random_instance_size(max_dimension
):
2379 The maximum rank of a random QuaternionHermitianEJA.
2381 # Obtained by solving d = 2n^2 - n.
2382 # The ZZ-int-ZZ thing is just "floor."
2383 return ZZ(int(ZZ(8*max_dimension
+ 1).sqrt()/4 + 1/4))
2386 def random_instance(cls
, max_dimension
=None, *args
, **kwargs
):
2388 Return a random instance of this type of algebra.
2390 class_max_d
= cls
._max
_random
_instance
_dimension
()
2391 if (max_dimension
is None or max_dimension
> class_max_d
):
2392 max_dimension
= class_max_d
2393 max_size
= cls
._max
_random
_instance
_size
(max_dimension
)
2394 n
= ZZ
.random_element(max_size
+ 1)
2395 return cls(n
, **kwargs
)
2397 class OctonionHermitianEJA(MatrixEJA
, RationalBasisEJA
, ConcreteEJA
):
2401 sage: from mjo.eja.eja_algebra import (FiniteDimensionalEJA,
2402 ....: OctonionHermitianEJA)
2403 sage: from mjo.hurwitz import Octonions, OctonionMatrixAlgebra
2407 The 3-by-3 algebra satisfies the axioms of an EJA::
2409 sage: OctonionHermitianEJA(3, # long time
2410 ....: field=QQ, # long time
2411 ....: orthonormalize=False, # long time
2412 ....: check_axioms=True) # long time
2413 Euclidean Jordan algebra of dimension 27 over Rational Field
2415 After a change-of-basis, the 2-by-2 algebra has the same
2416 multiplication table as the ten-dimensional Jordan spin algebra::
2418 sage: A = OctonionMatrixAlgebra(2,Octonions(QQ),QQ)
2419 sage: b = OctonionHermitianEJA._denormalized_basis(A)
2420 sage: basis = (b[0] + b[9],) + b[1:9] + (b[0] - b[9],)
2421 sage: jp = OctonionHermitianEJA.jordan_product
2422 sage: ip = OctonionHermitianEJA.trace_inner_product
2423 sage: J = FiniteDimensionalEJA(basis,
2427 ....: orthonormalize=False)
2428 sage: J.multiplication_table()
2429 +----++----+----+----+----+----+----+----+----+----+----+
2430 | * || b0 | b1 | b2 | b3 | b4 | b5 | b6 | b7 | b8 | b9 |
2431 +====++====+====+====+====+====+====+====+====+====+====+
2432 | b0 || b0 | b1 | b2 | b3 | b4 | b5 | b6 | b7 | b8 | b9 |
2433 +----++----+----+----+----+----+----+----+----+----+----+
2434 | b1 || b1 | b0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
2435 +----++----+----+----+----+----+----+----+----+----+----+
2436 | b2 || b2 | 0 | b0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
2437 +----++----+----+----+----+----+----+----+----+----+----+
2438 | b3 || b3 | 0 | 0 | b0 | 0 | 0 | 0 | 0 | 0 | 0 |
2439 +----++----+----+----+----+----+----+----+----+----+----+
2440 | b4 || b4 | 0 | 0 | 0 | b0 | 0 | 0 | 0 | 0 | 0 |
2441 +----++----+----+----+----+----+----+----+----+----+----+
2442 | b5 || b5 | 0 | 0 | 0 | 0 | b0 | 0 | 0 | 0 | 0 |
2443 +----++----+----+----+----+----+----+----+----+----+----+
2444 | b6 || b6 | 0 | 0 | 0 | 0 | 0 | b0 | 0 | 0 | 0 |
2445 +----++----+----+----+----+----+----+----+----+----+----+
2446 | b7 || b7 | 0 | 0 | 0 | 0 | 0 | 0 | b0 | 0 | 0 |
2447 +----++----+----+----+----+----+----+----+----+----+----+
2448 | b8 || b8 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | b0 | 0 |
2449 +----++----+----+----+----+----+----+----+----+----+----+
2450 | b9 || b9 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | b0 |
2451 +----++----+----+----+----+----+----+----+----+----+----+
2455 We can actually construct the 27-dimensional Albert algebra,
2456 and we get the right unit element if we recompute it::
2458 sage: J = OctonionHermitianEJA(3, # long time
2459 ....: field=QQ, # long time
2460 ....: orthonormalize=False) # long time
2461 sage: J.one.clear_cache() # long time
2462 sage: J.one() # long time
2464 sage: J.one().to_matrix() # long time
2473 The 2-by-2 algebra is isomorphic to the ten-dimensional Jordan
2474 spin algebra, but just to be sure, we recompute its rank::
2476 sage: J = OctonionHermitianEJA(2, # long time
2477 ....: field=QQ, # long time
2478 ....: orthonormalize=False) # long time
2479 sage: J.rank.clear_cache() # long time
2480 sage: J.rank() # long time
2485 def _max_random_instance_size(max_dimension
):
2487 The maximum rank of a random QuaternionHermitianEJA.
2489 # There's certainly a formula for this, but with only four
2490 # cases to worry about, I'm not that motivated to derive it.
2491 if max_dimension
>= 27:
2493 elif max_dimension
>= 10:
2495 elif max_dimension
>= 1:
2501 def random_instance(cls
, max_dimension
=None, *args
, **kwargs
):
2503 Return a random instance of this type of algebra.
2505 class_max_d
= cls
._max
_random
_instance
_dimension
()
2506 if (max_dimension
is None or max_dimension
> class_max_d
):
2507 max_dimension
= class_max_d
2508 max_size
= cls
._max
_random
_instance
_size
(max_dimension
)
2509 n
= ZZ
.random_element(max_size
+ 1)
2510 return cls(n
, **kwargs
)
2512 def __init__(self
, n
, field
=AA
, **kwargs
):
2514 # Otherwise we don't get an EJA.
2515 raise ValueError("n cannot exceed 3")
2517 # We know this is a valid EJA, but will double-check
2518 # if the user passes check_axioms=True.
2519 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2521 from mjo
.hurwitz
import OctonionMatrixAlgebra
2522 A
= OctonionMatrixAlgebra(n
, scalars
=field
)
2523 super().__init
__(A
, **kwargs
)
2525 from mjo
.eja
.eja_cache
import octonion_hermitian_eja_coeffs
2526 a
= octonion_hermitian_eja_coeffs(self
)
2528 self
.rational_algebra()._charpoly
_coefficients
.set_cache(a
)
2531 class AlbertEJA(OctonionHermitianEJA
):
2533 The Albert algebra is the algebra of three-by-three Hermitian
2534 matrices whose entries are octonions.
2538 sage: from mjo.eja.eja_algebra import AlbertEJA
2542 sage: AlbertEJA(field=QQ, orthonormalize=False)
2543 Euclidean Jordan algebra of dimension 27 over Rational Field
2544 sage: AlbertEJA() # long time
2545 Euclidean Jordan algebra of dimension 27 over Algebraic Real Field
2548 def __init__(self
, *args
, **kwargs
):
2549 super().__init
__(3, *args
, **kwargs
)
2552 class HadamardEJA(RationalBasisEJA
, ConcreteEJA
):
2554 Return the Euclidean Jordan algebra on `R^n` with the Hadamard
2555 (pointwise real-number multiplication) Jordan product and the
2556 usual inner-product.
2558 This is nothing more than the Cartesian product of ``n`` copies of
2559 the one-dimensional Jordan spin algebra, and is the most common
2560 example of a non-simple Euclidean Jordan algebra.
2564 sage: from mjo.eja.eja_algebra import HadamardEJA
2568 This multiplication table can be verified by hand::
2570 sage: J = HadamardEJA(3)
2571 sage: b0,b1,b2 = J.gens()
2587 We can change the generator prefix::
2589 sage: HadamardEJA(3, prefix='r').gens()
2592 def __init__(self
, n
, field
=AA
, **kwargs
):
2593 MS
= MatrixSpace(field
, n
, 1)
2596 jordan_product
= lambda x
,y
: x
2597 inner_product
= lambda x
,y
: x
2599 def jordan_product(x
,y
):
2600 return MS( xi
*yi
for (xi
,yi
) in zip(x
,y
) )
2602 def inner_product(x
,y
):
2605 # New defaults for keyword arguments. Don't orthonormalize
2606 # because our basis is already orthonormal with respect to our
2607 # inner-product. Don't check the axioms, because we know this
2608 # is a valid EJA... but do double-check if the user passes
2609 # check_axioms=True. Note: we DON'T override the "check_field"
2610 # default here, because the user can pass in a field!
2611 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2612 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2614 column_basis
= tuple( MS(b
) for b
in FreeModule(field
, n
).basis() )
2615 super().__init
__(column_basis
,
2622 self
.rank
.set_cache(n
)
2624 self
.one
.set_cache( self
.sum(self
.gens()) )
2627 def _max_random_instance_dimension():
2629 There's no reason to go higher than five here. That's
2630 enough to get the point across.
2635 def _max_random_instance_size(max_dimension
):
2637 The maximum size (=dimension) of a random HadamardEJA.
2639 return max_dimension
2642 def random_instance(cls
, max_dimension
=None, *args
, **kwargs
):
2644 Return a random instance of this type of algebra.
2646 class_max_d
= cls
._max
_random
_instance
_dimension
()
2647 if (max_dimension
is None or max_dimension
> class_max_d
):
2648 max_dimension
= class_max_d
2649 max_size
= cls
._max
_random
_instance
_size
(max_dimension
)
2650 n
= ZZ
.random_element(max_size
+ 1)
2651 return cls(n
, **kwargs
)
2654 class BilinearFormEJA(RationalBasisEJA
, ConcreteEJA
):
2656 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2657 with the half-trace inner product and jordan product ``x*y =
2658 (<Bx,y>,y_bar>, x0*y_bar + y0*x_bar)`` where `B = 1 \times B22` is
2659 a symmetric positive-definite "bilinear form" matrix. Its
2660 dimension is the size of `B`, and it has rank two in dimensions
2661 larger than two. It reduces to the ``JordanSpinEJA`` when `B` is
2662 the identity matrix of order ``n``.
2664 We insist that the one-by-one upper-left identity block of `B` be
2665 passed in as well so that we can be passed a matrix of size zero
2666 to construct a trivial algebra.
2670 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
2671 ....: JordanSpinEJA)
2675 When no bilinear form is specified, the identity matrix is used,
2676 and the resulting algebra is the Jordan spin algebra::
2678 sage: B = matrix.identity(AA,3)
2679 sage: J0 = BilinearFormEJA(B)
2680 sage: J1 = JordanSpinEJA(3)
2681 sage: J0.multiplication_table() == J0.multiplication_table()
2684 An error is raised if the matrix `B` does not correspond to a
2685 positive-definite bilinear form::
2687 sage: B = matrix.random(QQ,2,3)
2688 sage: J = BilinearFormEJA(B)
2689 Traceback (most recent call last):
2691 ValueError: bilinear form is not positive-definite
2692 sage: B = matrix.zero(QQ,3)
2693 sage: J = BilinearFormEJA(B)
2694 Traceback (most recent call last):
2696 ValueError: bilinear form is not positive-definite
2700 We can create a zero-dimensional algebra::
2702 sage: B = matrix.identity(AA,0)
2703 sage: J = BilinearFormEJA(B)
2707 We can check the multiplication condition given in the Jordan, von
2708 Neumann, and Wigner paper (and also discussed on my "On the
2709 symmetry..." paper). Note that this relies heavily on the standard
2710 choice of basis, as does anything utilizing the bilinear form
2711 matrix. We opt not to orthonormalize the basis, because if we
2712 did, we would have to normalize the `s_{i}` in a similar manner::
2714 sage: set_random_seed()
2715 sage: n = ZZ.random_element(5)
2716 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2717 sage: B11 = matrix.identity(QQ,1)
2718 sage: B22 = M.transpose()*M
2719 sage: B = block_matrix(2,2,[ [B11,0 ],
2721 sage: J = BilinearFormEJA(B, orthonormalize=False)
2722 sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
2723 sage: V = J.vector_space()
2724 sage: sis = [ J( V([0] + (M.inverse()*ei).list()).column() )
2725 ....: for ei in eis ]
2726 sage: actual = [ sis[i]*sis[j]
2727 ....: for i in range(n-1)
2728 ....: for j in range(n-1) ]
2729 sage: expected = [ J.one() if i == j else J.zero()
2730 ....: for i in range(n-1)
2731 ....: for j in range(n-1) ]
2732 sage: actual == expected
2736 def __init__(self
, B
, field
=AA
, **kwargs
):
2737 # The matrix "B" is supplied by the user in most cases,
2738 # so it makes sense to check whether or not its positive-
2739 # definite unless we are specifically asked not to...
2740 if ("check_axioms" not in kwargs
) or kwargs
["check_axioms"]:
2741 if not B
.is_positive_definite():
2742 raise ValueError("bilinear form is not positive-definite")
2744 # However, all of the other data for this EJA is computed
2745 # by us in manner that guarantees the axioms are
2746 # satisfied. So, again, unless we are specifically asked to
2747 # verify things, we'll skip the rest of the checks.
2748 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2751 MS
= MatrixSpace(field
, n
, 1)
2753 def inner_product(x
,y
):
2754 return (y
.T
*B
*x
)[0,0]
2756 def jordan_product(x
,y
):
2761 z0
= inner_product(y
,x
)
2762 zbar
= y0
*xbar
+ x0
*ybar
2763 return MS([z0
] + zbar
.list())
2765 column_basis
= tuple( MS(b
) for b
in FreeModule(field
, n
).basis() )
2767 # TODO: I haven't actually checked this, but it seems legit.
2772 super().__init
__(column_basis
,
2777 associative
=associative
,
2780 # The rank of this algebra is two, unless we're in a
2781 # one-dimensional ambient space (because the rank is bounded
2782 # by the ambient dimension).
2783 self
.rank
.set_cache(min(n
,2))
2785 self
.one
.set_cache( self
.zero() )
2787 self
.one
.set_cache( self
.monomial(0) )
2790 def _max_random_instance_dimension():
2792 There's no reason to go higher than five here. That's
2793 enough to get the point across.
2798 def _max_random_instance_size(max_dimension
):
2800 The maximum size (=dimension) of a random BilinearFormEJA.
2802 return max_dimension
2805 def random_instance(cls
, max_dimension
=None, *args
, **kwargs
):
2807 Return a random instance of this algebra.
2809 class_max_d
= cls
._max
_random
_instance
_dimension
()
2810 if (max_dimension
is None or max_dimension
> class_max_d
):
2811 max_dimension
= class_max_d
2812 max_size
= cls
._max
_random
_instance
_size
(max_dimension
)
2813 n
= ZZ
.random_element(max_size
+ 1)
2816 B
= matrix
.identity(ZZ
, n
)
2817 return cls(B
, **kwargs
)
2819 B11
= matrix
.identity(ZZ
, 1)
2820 M
= matrix
.random(ZZ
, n
-1)
2821 I
= matrix
.identity(ZZ
, n
-1)
2823 while alpha
.is_zero():
2824 alpha
= ZZ
.random_element().abs()
2826 B22
= M
.transpose()*M
+ alpha
*I
2828 from sage
.matrix
.special
import block_matrix
2829 B
= block_matrix(2,2, [ [B11
, ZZ(0) ],
2832 return cls(B
, **kwargs
)
2835 class JordanSpinEJA(BilinearFormEJA
):
2837 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2838 with the usual inner product and jordan product ``x*y =
2839 (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2844 sage: from mjo.eja.eja_algebra import JordanSpinEJA
2848 This multiplication table can be verified by hand::
2850 sage: J = JordanSpinEJA(4)
2851 sage: b0,b1,b2,b3 = J.gens()
2867 We can change the generator prefix::
2869 sage: JordanSpinEJA(2, prefix='B').gens()
2874 Ensure that we have the usual inner product on `R^n`::
2876 sage: set_random_seed()
2877 sage: J = JordanSpinEJA.random_instance()
2878 sage: x,y = J.random_elements(2)
2879 sage: actual = x.inner_product(y)
2880 sage: expected = x.to_vector().inner_product(y.to_vector())
2881 sage: actual == expected
2885 def __init__(self
, n
, *args
, **kwargs
):
2886 # This is a special case of the BilinearFormEJA with the
2887 # identity matrix as its bilinear form.
2888 B
= matrix
.identity(ZZ
, n
)
2890 # Don't orthonormalize because our basis is already
2891 # orthonormal with respect to our inner-product.
2892 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2894 # But also don't pass check_field=False here, because the user
2895 # can pass in a field!
2896 super().__init
__(B
, *args
, **kwargs
)
2899 def random_instance(cls
, max_dimension
=None, *args
, **kwargs
):
2901 Return a random instance of this type of algebra.
2903 Needed here to override the implementation for ``BilinearFormEJA``.
2905 class_max_d
= cls
._max
_random
_instance
_dimension
()
2906 if (max_dimension
is None or max_dimension
> class_max_d
):
2907 max_dimension
= class_max_d
2908 max_size
= cls
._max
_random
_instance
_size
(max_dimension
)
2909 n
= ZZ
.random_element(max_size
+ 1)
2910 return cls(n
, **kwargs
)
2913 class TrivialEJA(RationalBasisEJA
, ConcreteEJA
):
2915 The trivial Euclidean Jordan algebra consisting of only a zero element.
2919 sage: from mjo.eja.eja_algebra import TrivialEJA
2923 sage: J = TrivialEJA()
2930 sage: 7*J.one()*12*J.one()
2932 sage: J.one().inner_product(J.one())
2934 sage: J.one().norm()
2936 sage: J.one().subalgebra_generated_by()
2937 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2942 def __init__(self
, field
=AA
, **kwargs
):
2943 jordan_product
= lambda x
,y
: x
2944 inner_product
= lambda x
,y
: field
.zero()
2946 MS
= MatrixSpace(field
,0)
2948 # New defaults for keyword arguments
2949 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2950 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2952 super().__init
__(basis
,
2960 # The rank is zero using my definition, namely the dimension of the
2961 # largest subalgebra generated by any element.
2962 self
.rank
.set_cache(0)
2963 self
.one
.set_cache( self
.zero() )
2966 def random_instance(cls
, max_dimension
=None, *args
, **kwargs
):
2967 # We don't take a "size" argument so the superclass method is
2968 # inappropriate for us. The ``max_dimension`` argument is
2969 # included so that if this method is called generically with a
2970 # ``max_dimension=<whatever>`` argument, we don't try to pass
2971 # it on to the algebra constructor.
2972 return cls(**kwargs
)
2975 class CartesianProductEJA(FiniteDimensionalEJA
):
2977 The external (orthogonal) direct sum of two or more Euclidean
2978 Jordan algebras. Every Euclidean Jordan algebra decomposes into an
2979 orthogonal direct sum of simple Euclidean Jordan algebras which is
2980 then isometric to a Cartesian product, so no generality is lost by
2981 providing only this construction.
2985 sage: from mjo.eja.eja_algebra import (random_eja,
2986 ....: CartesianProductEJA,
2987 ....: ComplexHermitianEJA,
2989 ....: JordanSpinEJA,
2990 ....: RealSymmetricEJA)
2994 The Jordan product is inherited from our factors and implemented by
2995 our CombinatorialFreeModule Cartesian product superclass::
2997 sage: set_random_seed()
2998 sage: J1 = HadamardEJA(2)
2999 sage: J2 = RealSymmetricEJA(2)
3000 sage: J = cartesian_product([J1,J2])
3001 sage: x,y = J.random_elements(2)
3005 The ability to retrieve the original factors is implemented by our
3006 CombinatorialFreeModule Cartesian product superclass::
3008 sage: J1 = HadamardEJA(2, field=QQ)
3009 sage: J2 = JordanSpinEJA(3, field=QQ)
3010 sage: J = cartesian_product([J1,J2])
3011 sage: J.cartesian_factors()
3012 (Euclidean Jordan algebra of dimension 2 over Rational Field,
3013 Euclidean Jordan algebra of dimension 3 over Rational Field)
3015 You can provide more than two factors::
3017 sage: J1 = HadamardEJA(2)
3018 sage: J2 = JordanSpinEJA(3)
3019 sage: J3 = RealSymmetricEJA(3)
3020 sage: cartesian_product([J1,J2,J3])
3021 Euclidean Jordan algebra of dimension 2 over Algebraic Real
3022 Field (+) Euclidean Jordan algebra of dimension 3 over Algebraic
3023 Real Field (+) Euclidean Jordan algebra of dimension 6 over
3024 Algebraic Real Field
3026 Rank is additive on a Cartesian product::
3028 sage: J1 = HadamardEJA(1)
3029 sage: J2 = RealSymmetricEJA(2)
3030 sage: J = cartesian_product([J1,J2])
3031 sage: J1.rank.clear_cache()
3032 sage: J2.rank.clear_cache()
3033 sage: J.rank.clear_cache()
3036 sage: J.rank() == J1.rank() + J2.rank()
3039 The same rank computation works over the rationals, with whatever
3042 sage: J1 = HadamardEJA(1, field=QQ, orthonormalize=False)
3043 sage: J2 = RealSymmetricEJA(2, field=QQ, orthonormalize=False)
3044 sage: J = cartesian_product([J1,J2])
3045 sage: J1.rank.clear_cache()
3046 sage: J2.rank.clear_cache()
3047 sage: J.rank.clear_cache()
3050 sage: J.rank() == J1.rank() + J2.rank()
3053 The product algebra will be associative if and only if all of its
3054 components are associative::
3056 sage: J1 = HadamardEJA(2)
3057 sage: J1.is_associative()
3059 sage: J2 = HadamardEJA(3)
3060 sage: J2.is_associative()
3062 sage: J3 = RealSymmetricEJA(3)
3063 sage: J3.is_associative()
3065 sage: CP1 = cartesian_product([J1,J2])
3066 sage: CP1.is_associative()
3068 sage: CP2 = cartesian_product([J1,J3])
3069 sage: CP2.is_associative()
3072 Cartesian products of Cartesian products work::
3074 sage: J1 = JordanSpinEJA(1)
3075 sage: J2 = JordanSpinEJA(1)
3076 sage: J3 = JordanSpinEJA(1)
3077 sage: J = cartesian_product([J1,cartesian_product([J2,J3])])
3078 sage: J.multiplication_table()
3079 +----++----+----+----+
3080 | * || b0 | b1 | b2 |
3081 +====++====+====+====+
3082 | b0 || b0 | 0 | 0 |
3083 +----++----+----+----+
3084 | b1 || 0 | b1 | 0 |
3085 +----++----+----+----+
3086 | b2 || 0 | 0 | b2 |
3087 +----++----+----+----+
3088 sage: HadamardEJA(3).multiplication_table()
3089 +----++----+----+----+
3090 | * || b0 | b1 | b2 |
3091 +====++====+====+====+
3092 | b0 || b0 | 0 | 0 |
3093 +----++----+----+----+
3094 | b1 || 0 | b1 | 0 |
3095 +----++----+----+----+
3096 | b2 || 0 | 0 | b2 |
3097 +----++----+----+----+
3099 The "matrix space" of a Cartesian product always consists of
3100 ordered pairs (or triples, or...) whose components are the
3101 matrix spaces of its factors::
3103 sage: J1 = HadamardEJA(2)
3104 sage: J2 = ComplexHermitianEJA(2)
3105 sage: J = cartesian_product([J1,J2])
3106 sage: J.matrix_space()
3107 The Cartesian product of (Full MatrixSpace of 2 by 1 dense
3108 matrices over Algebraic Real Field, Module of 2 by 2 matrices
3109 with entries in Algebraic Field over the scalar ring Algebraic
3111 sage: J.one().to_matrix()[0]
3114 sage: J.one().to_matrix()[1]
3123 All factors must share the same base field::
3125 sage: J1 = HadamardEJA(2, field=QQ)
3126 sage: J2 = RealSymmetricEJA(2)
3127 sage: CartesianProductEJA((J1,J2))
3128 Traceback (most recent call last):
3130 ValueError: all factors must share the same base field
3132 The cached unit element is the same one that would be computed::
3134 sage: set_random_seed() # long time
3135 sage: J1 = random_eja() # long time
3136 sage: J2 = random_eja() # long time
3137 sage: J = cartesian_product([J1,J2]) # long time
3138 sage: actual = J.one() # long time
3139 sage: J.one.clear_cache() # long time
3140 sage: expected = J.one() # long time
3141 sage: actual == expected # long time
3144 Element
= CartesianProductEJAElement
3145 def __init__(self
, factors
, **kwargs
):
3150 self
._sets
= factors
3152 field
= factors
[0].base_ring()
3153 if not all( J
.base_ring() == field
for J
in factors
):
3154 raise ValueError("all factors must share the same base field")
3156 # Figure out the category to use.
3157 associative
= all( f
.is_associative() for f
in factors
)
3158 category
= EuclideanJordanAlgebras(field
)
3159 if associative
: category
= category
.Associative()
3160 category
= category
.join([category
, category
.CartesianProducts()])
3162 # Compute my matrix space. We don't simply use the
3163 # ``cartesian_product()`` functor here because it acts
3164 # differently on SageMath MatrixSpaces and our custom
3165 # MatrixAlgebras, which are CombinatorialFreeModules. We
3166 # always want the result to be represented (and indexed) as an
3167 # ordered tuple. This category isn't perfect, but is good
3168 # enough for what we need to do.
3169 MS_cat
= MagmaticAlgebras(field
).FiniteDimensional().WithBasis()
3170 MS_cat
= MS_cat
.Unital().CartesianProducts()
3171 MS_factors
= tuple( J
.matrix_space() for J
in factors
)
3172 from sage
.sets
.cartesian_product
import CartesianProduct
3173 self
._matrix
_space
= CartesianProduct(MS_factors
, MS_cat
)
3175 self
._matrix
_basis
= []
3176 zero
= self
._matrix
_space
.zero()
3178 for b
in factors
[i
].matrix_basis():
3181 self
._matrix
_basis
.append(z
)
3183 self
._matrix
_basis
= tuple( self
._matrix
_space
(b
)
3184 for b
in self
._matrix
_basis
)
3185 n
= len(self
._matrix
_basis
)
3187 # We already have what we need for the super-superclass constructor.
3188 CombinatorialFreeModule
.__init
__(self
,
3195 # Now create the vector space for the algebra, which will have
3196 # its own set of non-ambient coordinates (in terms of the
3198 degree
= sum( f
._matrix
_span
.ambient_vector_space().degree()
3200 V
= VectorSpace(field
, degree
)
3201 vector_basis
= tuple( V(_all2list(b
)) for b
in self
._matrix
_basis
)
3203 # Save the span of our matrix basis (when written out as long
3204 # vectors) because otherwise we'll have to reconstruct it
3205 # every time we want to coerce a matrix into the algebra.
3206 self
._matrix
_span
= V
.span_of_basis( vector_basis
, check
=False)
3208 # Since we don't (re)orthonormalize the basis, the FDEJA
3209 # constructor is going to set self._deortho_matrix to the
3210 # identity matrix. Here we set it to the correct value using
3211 # the deortho matrices from our factors.
3212 self
._deortho
_matrix
= matrix
.block_diagonal(
3213 [J
._deortho
_matrix
for J
in factors
]
3216 self
._inner
_product
_matrix
= matrix
.block_diagonal(
3217 [J
._inner
_product
_matrix
for J
in factors
]
3219 self
._inner
_product
_matrix
._cache
= {'hermitian': True}
3220 self
._inner
_product
_matrix
.set_immutable()
3222 # Building the multiplication table is a bit more tricky
3223 # because we have to embed the entries of the factors'
3224 # multiplication tables into the product EJA.
3226 self
._multiplication
_table
= [ [zed
for j
in range(i
+1)]
3229 # Keep track of an offset that tallies the dimensions of all
3230 # previous factors. If the second factor is dim=2 and if the
3231 # first one is dim=3, then we want to skip the first 3x3 block
3232 # when copying the multiplication table for the second factor.
3235 phi_f
= self
.cartesian_embedding(f
)
3236 factor_dim
= factors
[f
].dimension()
3237 for i
in range(factor_dim
):
3238 for j
in range(i
+1):
3239 f_ij
= factors
[f
]._multiplication
_table
[i
][j
]
3241 self
._multiplication
_table
[offset
+i
][offset
+j
] = e
3242 offset
+= factor_dim
3244 self
.rank
.set_cache(sum(J
.rank() for J
in factors
))
3245 ones
= tuple(J
.one().to_matrix() for J
in factors
)
3246 self
.one
.set_cache(self(ones
))
3248 def _sets_keys(self
):
3253 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
3254 ....: RealSymmetricEJA)
3258 The superclass uses ``_sets_keys()`` to implement its
3259 ``cartesian_factors()`` method::
3261 sage: J1 = RealSymmetricEJA(2,
3263 ....: orthonormalize=False,
3265 sage: J2 = ComplexHermitianEJA(2,field=QQ,orthonormalize=False)
3266 sage: J = cartesian_product([J1,J2])
3267 sage: x = sum(i*J.gens()[i] for i in range(len(J.gens())))
3268 sage: x.cartesian_factors()
3269 (a1 + 2*a2, 3*b0 + 4*b1 + 5*b2 + 6*b3)
3272 # Copy/pasted from CombinatorialFreeModule_CartesianProduct,
3273 # but returning a tuple instead of a list.
3274 return tuple(range(len(self
.cartesian_factors())))
3276 def cartesian_factors(self
):
3277 # Copy/pasted from CombinatorialFreeModule_CartesianProduct.
3280 def cartesian_factor(self
, i
):
3282 Return the ``i``th factor of this algebra.
3284 return self
._sets
[i
]
3287 # Copy/pasted from CombinatorialFreeModule_CartesianProduct.
3288 from sage
.categories
.cartesian_product
import cartesian_product
3289 return cartesian_product
.symbol
.join("%s" % factor
3290 for factor
in self
._sets
)
3294 def cartesian_projection(self
, i
):
3298 sage: from mjo.eja.eja_algebra import (random_eja,
3299 ....: JordanSpinEJA,
3301 ....: RealSymmetricEJA,
3302 ....: ComplexHermitianEJA)
3306 The projection morphisms are Euclidean Jordan algebra
3309 sage: J1 = HadamardEJA(2)
3310 sage: J2 = RealSymmetricEJA(2)
3311 sage: J = cartesian_product([J1,J2])
3312 sage: J.cartesian_projection(0)
3313 Linear operator between finite-dimensional Euclidean Jordan
3314 algebras represented by the matrix:
3317 Domain: Euclidean Jordan algebra of dimension 2 over Algebraic
3318 Real Field (+) Euclidean Jordan algebra of dimension 3 over
3319 Algebraic Real Field
3320 Codomain: Euclidean Jordan algebra of dimension 2 over Algebraic
3322 sage: J.cartesian_projection(1)
3323 Linear operator between finite-dimensional Euclidean Jordan
3324 algebras represented by the matrix:
3328 Domain: Euclidean Jordan algebra of dimension 2 over Algebraic
3329 Real Field (+) Euclidean Jordan algebra of dimension 3 over
3330 Algebraic Real Field
3331 Codomain: Euclidean Jordan algebra of dimension 3 over Algebraic
3334 The projections work the way you'd expect on the vector
3335 representation of an element::
3337 sage: J1 = JordanSpinEJA(2)
3338 sage: J2 = ComplexHermitianEJA(2)
3339 sage: J = cartesian_product([J1,J2])
3340 sage: pi_left = J.cartesian_projection(0)
3341 sage: pi_right = J.cartesian_projection(1)
3342 sage: pi_left(J.one()).to_vector()
3344 sage: pi_right(J.one()).to_vector()
3346 sage: J.one().to_vector()
3351 The answer never changes::
3353 sage: set_random_seed()
3354 sage: J1 = random_eja()
3355 sage: J2 = random_eja()
3356 sage: J = cartesian_product([J1,J2])
3357 sage: P0 = J.cartesian_projection(0)
3358 sage: P1 = J.cartesian_projection(0)
3363 offset
= sum( self
.cartesian_factor(k
).dimension()
3365 Ji
= self
.cartesian_factor(i
)
3366 Pi
= self
._module
_morphism
(lambda j
: Ji
.monomial(j
- offset
),
3369 return FiniteDimensionalEJAOperator(self
,Ji
,Pi
.matrix())
3372 def cartesian_embedding(self
, i
):
3376 sage: from mjo.eja.eja_algebra import (random_eja,
3377 ....: JordanSpinEJA,
3379 ....: RealSymmetricEJA)
3383 The embedding morphisms are Euclidean Jordan algebra
3386 sage: J1 = HadamardEJA(2)
3387 sage: J2 = RealSymmetricEJA(2)
3388 sage: J = cartesian_product([J1,J2])
3389 sage: J.cartesian_embedding(0)
3390 Linear operator between finite-dimensional Euclidean Jordan
3391 algebras represented by the matrix:
3397 Domain: Euclidean Jordan algebra of dimension 2 over
3398 Algebraic Real Field
3399 Codomain: Euclidean Jordan algebra of dimension 2 over
3400 Algebraic Real Field (+) Euclidean Jordan algebra of
3401 dimension 3 over Algebraic Real Field
3402 sage: J.cartesian_embedding(1)
3403 Linear operator between finite-dimensional Euclidean Jordan
3404 algebras represented by the matrix:
3410 Domain: Euclidean Jordan algebra of dimension 3 over
3411 Algebraic Real Field
3412 Codomain: Euclidean Jordan algebra of dimension 2 over
3413 Algebraic Real Field (+) Euclidean Jordan algebra of
3414 dimension 3 over Algebraic Real Field
3416 The embeddings work the way you'd expect on the vector
3417 representation of an element::
3419 sage: J1 = JordanSpinEJA(3)
3420 sage: J2 = RealSymmetricEJA(2)
3421 sage: J = cartesian_product([J1,J2])
3422 sage: iota_left = J.cartesian_embedding(0)
3423 sage: iota_right = J.cartesian_embedding(1)
3424 sage: iota_left(J1.zero()) == J.zero()
3426 sage: iota_right(J2.zero()) == J.zero()
3428 sage: J1.one().to_vector()
3430 sage: iota_left(J1.one()).to_vector()
3432 sage: J2.one().to_vector()
3434 sage: iota_right(J2.one()).to_vector()
3436 sage: J.one().to_vector()
3441 The answer never changes::
3443 sage: set_random_seed()
3444 sage: J1 = random_eja()
3445 sage: J2 = random_eja()
3446 sage: J = cartesian_product([J1,J2])
3447 sage: E0 = J.cartesian_embedding(0)
3448 sage: E1 = J.cartesian_embedding(0)
3452 Composing a projection with the corresponding inclusion should
3453 produce the identity map, and mismatching them should produce
3456 sage: set_random_seed()
3457 sage: J1 = random_eja()
3458 sage: J2 = random_eja()
3459 sage: J = cartesian_product([J1,J2])
3460 sage: iota_left = J.cartesian_embedding(0)
3461 sage: iota_right = J.cartesian_embedding(1)
3462 sage: pi_left = J.cartesian_projection(0)
3463 sage: pi_right = J.cartesian_projection(1)
3464 sage: pi_left*iota_left == J1.one().operator()
3466 sage: pi_right*iota_right == J2.one().operator()
3468 sage: (pi_left*iota_right).is_zero()
3470 sage: (pi_right*iota_left).is_zero()
3474 offset
= sum( self
.cartesian_factor(k
).dimension()
3476 Ji
= self
.cartesian_factor(i
)
3477 Ei
= Ji
._module
_morphism
(lambda j
: self
.monomial(j
+ offset
),
3479 return FiniteDimensionalEJAOperator(Ji
,self
,Ei
.matrix())
3483 FiniteDimensionalEJA
.CartesianProduct
= CartesianProductEJA
3485 class RationalBasisCartesianProductEJA(CartesianProductEJA
,
3488 A separate class for products of algebras for which we know a
3493 sage: from mjo.eja.eja_algebra import (FiniteDimensionalEJA,
3495 ....: JordanSpinEJA,
3496 ....: RealSymmetricEJA)
3500 This gives us fast characteristic polynomial computations in
3501 product algebras, too::
3504 sage: J1 = JordanSpinEJA(2)
3505 sage: J2 = RealSymmetricEJA(3)
3506 sage: J = cartesian_product([J1,J2])
3507 sage: J.characteristic_polynomial_of().degree()
3514 The ``cartesian_product()`` function only uses the first factor to
3515 decide where the result will live; thus we have to be careful to
3516 check that all factors do indeed have a ``rational_algebra()`` method
3517 before we construct an algebra that claims to have a rational basis::
3519 sage: J1 = HadamardEJA(2)
3520 sage: jp = lambda X,Y: X*Y
3521 sage: ip = lambda X,Y: X[0,0]*Y[0,0]
3522 sage: b1 = matrix(QQ, [[1]])
3523 sage: J2 = FiniteDimensionalEJA((b1,), jp, ip)
3524 sage: cartesian_product([J2,J1]) # factor one not RationalBasisEJA
3525 Euclidean Jordan algebra of dimension 1 over Algebraic Real
3526 Field (+) Euclidean Jordan algebra of dimension 2 over Algebraic
3528 sage: cartesian_product([J1,J2]) # factor one is RationalBasisEJA
3529 Traceback (most recent call last):
3531 ValueError: factor not a RationalBasisEJA
3534 def __init__(self
, algebras
, **kwargs
):
3535 if not all( hasattr(r
, "rational_algebra") for r
in algebras
):
3536 raise ValueError("factor not a RationalBasisEJA")
3538 CartesianProductEJA
.__init
__(self
, algebras
, **kwargs
)
3541 def rational_algebra(self
):
3542 if self
.base_ring() is QQ
:
3545 return cartesian_product([
3546 r
.rational_algebra() for r
in self
.cartesian_factors()
3550 RationalBasisEJA
.CartesianProduct
= RationalBasisCartesianProductEJA
3552 def random_eja(max_dimension
=None, *args
, **kwargs
):
3557 sage: from mjo.eja.eja_algebra import random_eja
3561 sage: set_random_seed()
3562 sage: n = ZZ.random_element(1,5)
3563 sage: J = random_eja(max_dimension=n, field=QQ, orthonormalize=False)
3564 sage: J.dimension() <= n
3568 # Use the ConcreteEJA default as the total upper bound (regardless
3569 # of any whether or not any individual factors set a lower limit).
3570 if max_dimension
is None:
3571 max_dimension
= ConcreteEJA
._max
_random
_instance
_dimension
()
3572 J1
= ConcreteEJA
.random_instance(max_dimension
, *args
, **kwargs
)
3575 # Roll the dice to see if we attempt a Cartesian product.
3576 dice_roll
= ZZ
.random_element(len(ConcreteEJA
.__subclasses
__()) + 1)
3577 new_max_dimension
= max_dimension
- J1
.dimension()
3578 if new_max_dimension
== 0 or dice_roll
!= 0:
3579 # If it's already as big as we're willing to tolerate, just
3580 # return it and don't worry about Cartesian products.
3583 # Use random_eja() again so we can get more than two factors
3584 # if the sub-call also Decides on a cartesian product.
3585 J2
= random_eja(new_max_dimension
, *args
, **kwargs
)
3586 return cartesian_product([J1
,J2
])