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.
57 sage: from mjo.eja.eja_algebra import random_eja
62 Euclidean Jordan algebra of dimension...
65 from sage
.algebras
.quatalg
.quaternion_algebra
import QuaternionAlgebra
66 from sage
.categories
.magmatic_algebras
import MagmaticAlgebras
67 from sage
.categories
.sets_cat
import cartesian_product
68 from sage
.combinat
.free_module
import CombinatorialFreeModule
69 from sage
.matrix
.constructor
import matrix
70 from sage
.matrix
.matrix_space
import MatrixSpace
71 from sage
.misc
.cachefunc
import cached_method
72 from sage
.misc
.table
import table
73 from sage
.modules
.free_module
import FreeModule
, VectorSpace
74 from sage
.rings
.all
import (ZZ
, QQ
, AA
, QQbar
, RR
, RLF
, CLF
,
77 from mjo
.eja
.eja_element
import FiniteDimensionalEJAElement
78 from mjo
.eja
.eja_operator
import FiniteDimensionalEJAOperator
79 from mjo
.eja
.eja_utils
import _all2list
, _mat2vec
81 class FiniteDimensionalEJA(CombinatorialFreeModule
):
83 A finite-dimensional Euclidean Jordan algebra.
87 - ``basis`` -- a tuple; a tuple of basis elements in "matrix
88 form," which must be the same form as the arguments to
89 ``jordan_product`` and ``inner_product``. In reality, "matrix
90 form" can be either vectors, matrices, or a Cartesian product
91 (ordered tuple) of vectors or matrices. All of these would
92 ideally be vector spaces in sage with no special-casing
93 needed; but in reality we turn vectors into column-matrices
94 and Cartesian products `(a,b)` into column matrices
95 `(a,b)^{T}` after converting `a` and `b` themselves.
97 - ``jordan_product`` -- a function; afunction of two ``basis``
98 elements (in matrix form) that returns their jordan product,
99 also in matrix form; this will be applied to ``basis`` to
100 compute a multiplication table for the algebra.
102 - ``inner_product`` -- a function; a function of two ``basis``
103 elements (in matrix form) that returns their inner
104 product. This will be applied to ``basis`` to compute an
105 inner-product table (basically a matrix) for this algebra.
107 - ``field`` -- a subfield of the reals (default: ``AA``); the scalar
108 field for the algebra.
110 - ``orthonormalize`` -- boolean (default: ``True``); whether or
111 not to orthonormalize the basis. Doing so is expensive and
112 generally rules out using the rationals as your ``field``, but
113 is required for spectral decompositions.
117 sage: from mjo.eja.eja_algebra import random_eja
121 We should compute that an element subalgebra is associative even
122 if we circumvent the element method::
124 sage: set_random_seed()
125 sage: J = random_eja(field=QQ,orthonormalize=False)
126 sage: x = J.random_element()
127 sage: A = x.subalgebra_generated_by(orthonormalize=False)
128 sage: basis = tuple(b.superalgebra_element() for b in A.basis())
129 sage: J.subalgebra(basis, orthonormalize=False).is_associative()
133 Element
= FiniteDimensionalEJAElement
142 cartesian_product
=False,
150 if not field
.is_subring(RR
):
151 # Note: this does return true for the real algebraic
152 # field, the rationals, and any quadratic field where
153 # we've specified a real embedding.
154 raise ValueError("scalar field is not real")
157 # Check commutativity of the Jordan and inner-products.
158 # This has to be done before we build the multiplication
159 # and inner-product tables/matrices, because we take
160 # advantage of symmetry in the process.
161 if not all( jordan_product(bi
,bj
) == jordan_product(bj
,bi
)
164 raise ValueError("Jordan product is not commutative")
166 if not all( inner_product(bi
,bj
) == inner_product(bj
,bi
)
169 raise ValueError("inner-product is not commutative")
172 category
= MagmaticAlgebras(field
).FiniteDimensional()
173 category
= category
.WithBasis().Unital().Commutative()
175 if associative
is None:
176 # We should figure it out. As with check_axioms, we have to do
177 # this without the help of the _jordan_product_is_associative()
178 # method because we need to know the category before we
179 # initialize the algebra.
180 associative
= all( jordan_product(jordan_product(bi
,bj
),bk
)
182 jordan_product(bi
,jordan_product(bj
,bk
))
188 # Element subalgebras can take advantage of this.
189 category
= category
.Associative()
190 if cartesian_product
:
191 # Use join() here because otherwise we only get the
192 # "Cartesian product of..." and not the things themselves.
193 category
= category
.join([category
,
194 category
.CartesianProducts()])
196 # Call the superclass constructor so that we can use its from_vector()
197 # method to build our multiplication table.
198 CombinatorialFreeModule
.__init
__(self
,
205 # Now comes all of the hard work. We'll be constructing an
206 # ambient vector space V that our (vectorized) basis lives in,
207 # as well as a subspace W of V spanned by those (vectorized)
208 # basis elements. The W-coordinates are the coefficients that
209 # we see in things like x = 1*b1 + 2*b2.
214 degree
= len(_all2list(basis
[0]))
216 # Build an ambient space that fits our matrix basis when
217 # written out as "long vectors."
218 V
= VectorSpace(field
, degree
)
220 # The matrix that will hole the orthonormal -> unorthonormal
221 # coordinate transformation.
222 self
._deortho
_matrix
= None
225 # Save a copy of the un-orthonormalized basis for later.
226 # Convert it to ambient V (vector) coordinates while we're
227 # at it, because we'd have to do it later anyway.
228 deortho_vector_basis
= tuple( V(_all2list(b
)) for b
in basis
)
230 from mjo
.eja
.eja_utils
import gram_schmidt
231 basis
= tuple(gram_schmidt(basis
, inner_product
))
233 # Save the (possibly orthonormalized) matrix basis for
235 self
._matrix
_basis
= basis
237 # Now create the vector space for the algebra, which will have
238 # its own set of non-ambient coordinates (in terms of the
240 vector_basis
= tuple( V(_all2list(b
)) for b
in basis
)
241 W
= V
.span_of_basis( vector_basis
, check
=check_axioms
)
244 # Now "W" is the vector space of our algebra coordinates. The
245 # variables "X1", "X2",... refer to the entries of vectors in
246 # W. Thus to convert back and forth between the orthonormal
247 # coordinates and the given ones, we need to stick the original
249 U
= V
.span_of_basis( deortho_vector_basis
, check
=check_axioms
)
250 self
._deortho
_matrix
= matrix( U
.coordinate_vector(q
)
251 for q
in vector_basis
)
254 # Now we actually compute the multiplication and inner-product
255 # tables/matrices using the possibly-orthonormalized basis.
256 self
._inner
_product
_matrix
= matrix
.identity(field
, n
)
257 self
._multiplication
_table
= [ [0 for j
in range(i
+1)]
260 # Note: the Jordan and inner-products are defined in terms
261 # of the ambient basis. It's important that their arguments
262 # are in ambient coordinates as well.
265 # ortho basis w.r.t. ambient coords
269 # The jordan product returns a matrixy answer, so we
270 # have to convert it to the algebra coordinates.
271 elt
= jordan_product(q_i
, q_j
)
272 elt
= W
.coordinate_vector(V(_all2list(elt
)))
273 self
._multiplication
_table
[i
][j
] = self
.from_vector(elt
)
275 if not orthonormalize
:
276 # If we're orthonormalizing the basis with respect
277 # to an inner-product, then the inner-product
278 # matrix with respect to the resulting basis is
279 # just going to be the identity.
280 ip
= inner_product(q_i
, q_j
)
281 self
._inner
_product
_matrix
[i
,j
] = ip
282 self
._inner
_product
_matrix
[j
,i
] = ip
284 self
._inner
_product
_matrix
._cache
= {'hermitian': True}
285 self
._inner
_product
_matrix
.set_immutable()
288 if not self
._is
_jordanian
():
289 raise ValueError("Jordan identity does not hold")
290 if not self
._inner
_product
_is
_associative
():
291 raise ValueError("inner product is not associative")
294 def _coerce_map_from_base_ring(self
):
296 Disable the map from the base ring into the algebra.
298 Performing a nonsense conversion like this automatically
299 is counterpedagogical. The fallback is to try the usual
300 element constructor, which should also fail.
304 sage: from mjo.eja.eja_algebra import random_eja
308 sage: set_random_seed()
309 sage: J = random_eja()
311 Traceback (most recent call last):
313 ValueError: not an element of this algebra
319 def product_on_basis(self
, i
, j
):
321 Returns the Jordan product of the `i` and `j`th basis elements.
323 This completely defines the Jordan product on the algebra, and
324 is used direclty by our superclass machinery to implement
329 sage: from mjo.eja.eja_algebra import random_eja
333 sage: set_random_seed()
334 sage: J = random_eja()
335 sage: n = J.dimension()
338 sage: bi_bj = J.zero()*J.zero()
340 ....: i = ZZ.random_element(n)
341 ....: j = ZZ.random_element(n)
342 ....: bi = J.monomial(i)
343 ....: bj = J.monomial(j)
344 ....: bi_bj = J.product_on_basis(i,j)
349 # We only stored the lower-triangular portion of the
350 # multiplication table.
352 return self
._multiplication
_table
[i
][j
]
354 return self
._multiplication
_table
[j
][i
]
356 def inner_product(self
, x
, y
):
358 The inner product associated with this Euclidean Jordan algebra.
360 Defaults to the trace inner product, but can be overridden by
361 subclasses if they are sure that the necessary properties are
366 sage: from mjo.eja.eja_algebra import (random_eja,
368 ....: BilinearFormEJA)
372 Our inner product is "associative," which means the following for
373 a symmetric bilinear form::
375 sage: set_random_seed()
376 sage: J = random_eja()
377 sage: x,y,z = J.random_elements(3)
378 sage: (x*y).inner_product(z) == y.inner_product(x*z)
383 Ensure that this is the usual inner product for the algebras
386 sage: set_random_seed()
387 sage: J = HadamardEJA.random_instance()
388 sage: x,y = J.random_elements(2)
389 sage: actual = x.inner_product(y)
390 sage: expected = x.to_vector().inner_product(y.to_vector())
391 sage: actual == expected
394 Ensure that this is one-half of the trace inner-product in a
395 BilinearFormEJA that isn't just the reals (when ``n`` isn't
396 one). This is in Faraut and Koranyi, and also my "On the
399 sage: set_random_seed()
400 sage: J = BilinearFormEJA.random_instance()
401 sage: n = J.dimension()
402 sage: x = J.random_element()
403 sage: y = J.random_element()
404 sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
408 B
= self
._inner
_product
_matrix
409 return (B
*x
.to_vector()).inner_product(y
.to_vector())
412 def is_associative(self
):
414 Return whether or not this algebra's Jordan product is associative.
418 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
422 sage: J = ComplexHermitianEJA(3, field=QQ, orthonormalize=False)
423 sage: J.is_associative()
425 sage: x = sum(J.gens())
426 sage: A = x.subalgebra_generated_by(orthonormalize=False)
427 sage: A.is_associative()
431 return "Associative" in self
.category().axioms()
433 def _is_commutative(self
):
435 Whether or not this algebra's multiplication table is commutative.
437 This method should of course always return ``True``, unless
438 this algebra was constructed with ``check_axioms=False`` and
439 passed an invalid multiplication table.
441 return all( x
*y
== y
*x
for x
in self
.gens() for y
in self
.gens() )
443 def _is_jordanian(self
):
445 Whether or not this algebra's multiplication table respects the
446 Jordan identity `(x^{2})(xy) = x(x^{2}y)`.
448 We only check one arrangement of `x` and `y`, so for a
449 ``True`` result to be truly true, you should also check
450 :meth:`_is_commutative`. This method should of course always
451 return ``True``, unless this algebra was constructed with
452 ``check_axioms=False`` and passed an invalid multiplication table.
454 return all( (self
.monomial(i
)**2)*(self
.monomial(i
)*self
.monomial(j
))
456 (self
.monomial(i
))*((self
.monomial(i
)**2)*self
.monomial(j
))
457 for i
in range(self
.dimension())
458 for j
in range(self
.dimension()) )
460 def _jordan_product_is_associative(self
):
462 Return whether or not this algebra's Jordan product is
463 associative; that is, whether or not `x*(y*z) = (x*y)*z`
466 This method should agree with :meth:`is_associative` unless
467 you lied about the value of the ``associative`` parameter
468 when you constructed the algebra.
472 sage: from mjo.eja.eja_algebra import (random_eja,
473 ....: RealSymmetricEJA,
474 ....: ComplexHermitianEJA,
475 ....: QuaternionHermitianEJA)
479 sage: J = RealSymmetricEJA(4, orthonormalize=False)
480 sage: J._jordan_product_is_associative()
482 sage: x = sum(J.gens())
483 sage: A = x.subalgebra_generated_by()
484 sage: A._jordan_product_is_associative()
489 sage: J = ComplexHermitianEJA(2,field=QQ,orthonormalize=False)
490 sage: J._jordan_product_is_associative()
492 sage: x = sum(J.gens())
493 sage: A = x.subalgebra_generated_by(orthonormalize=False)
494 sage: A._jordan_product_is_associative()
499 sage: J = QuaternionHermitianEJA(2)
500 sage: J._jordan_product_is_associative()
502 sage: x = sum(J.gens())
503 sage: A = x.subalgebra_generated_by()
504 sage: A._jordan_product_is_associative()
509 The values we've presupplied to the constructors agree with
512 sage: set_random_seed()
513 sage: J = random_eja()
514 sage: J.is_associative() == J._jordan_product_is_associative()
520 # Used to check whether or not something is zero.
523 # I don't know of any examples that make this magnitude
524 # necessary because I don't know how to make an
525 # associative algebra when the element subalgebra
526 # construction is unreliable (as it is over RDF; we can't
527 # find the degree of an element because we can't compute
528 # the rank of a matrix). But even multiplication of floats
529 # is non-associative, so *some* epsilon is needed... let's
530 # just take the one from _inner_product_is_associative?
533 for i
in range(self
.dimension()):
534 for j
in range(self
.dimension()):
535 for k
in range(self
.dimension()):
539 diff
= (x
*y
)*z
- x
*(y
*z
)
541 if diff
.norm() > epsilon
:
546 def _inner_product_is_associative(self
):
548 Return whether or not this algebra's inner product `B` is
549 associative; that is, whether or not `B(xy,z) = B(x,yz)`.
551 This method should of course always return ``True``, unless
552 this algebra was constructed with ``check_axioms=False`` and
553 passed an invalid Jordan or inner-product.
557 # Used to check whether or not something is zero.
560 # This choice is sufficient to allow the construction of
561 # QuaternionHermitianEJA(2, field=RDF) with check_axioms=True.
564 for i
in range(self
.dimension()):
565 for j
in range(self
.dimension()):
566 for k
in range(self
.dimension()):
570 diff
= (x
*y
).inner_product(z
) - x
.inner_product(y
*z
)
572 if diff
.abs() > epsilon
:
577 def _element_constructor_(self
, elt
):
579 Construct an element of this algebra from its vector or matrix
582 This gets called only after the parent element _call_ method
583 fails to find a coercion for the argument.
587 sage: from mjo.eja.eja_algebra import (random_eja,
590 ....: RealSymmetricEJA)
594 The identity in `S^n` is converted to the identity in the EJA::
596 sage: J = RealSymmetricEJA(3)
597 sage: I = matrix.identity(QQ,3)
598 sage: J(I) == J.one()
601 This skew-symmetric matrix can't be represented in the EJA::
603 sage: J = RealSymmetricEJA(3)
604 sage: A = matrix(QQ,3, lambda i,j: i-j)
606 Traceback (most recent call last):
608 ValueError: not an element of this algebra
610 Tuples work as well, provided that the matrix basis for the
611 algebra consists of them::
613 sage: J1 = HadamardEJA(3)
614 sage: J2 = RealSymmetricEJA(2)
615 sage: J = cartesian_product([J1,J2])
616 sage: J( (J1.matrix_basis()[1], J2.matrix_basis()[2]) )
621 Ensure that we can convert any element back and forth
622 faithfully between its matrix and algebra representations::
624 sage: set_random_seed()
625 sage: J = random_eja()
626 sage: x = J.random_element()
627 sage: J(x.to_matrix()) == x
630 We cannot coerce elements between algebras just because their
631 matrix representations are compatible::
633 sage: J1 = HadamardEJA(3)
634 sage: J2 = JordanSpinEJA(3)
636 Traceback (most recent call last):
638 ValueError: not an element of this algebra
640 Traceback (most recent call last):
642 ValueError: not an element of this algebra
644 msg
= "not an element of this algebra"
645 if elt
in self
.base_ring():
646 # Ensure that no base ring -> algebra coercion is performed
647 # by this method. There's some stupidity in sage that would
648 # otherwise propagate to this method; for example, sage thinks
649 # that the integer 3 belongs to the space of 2-by-2 matrices.
650 raise ValueError(msg
)
653 # Try to convert a vector into a column-matrix...
655 except (AttributeError, TypeError):
656 # and ignore failure, because we weren't really expecting
657 # a vector as an argument anyway.
660 if elt
not in self
.matrix_space():
661 raise ValueError(msg
)
663 # Thanks for nothing! Matrix spaces aren't vector spaces in
664 # Sage, so we have to figure out its matrix-basis coordinates
665 # ourselves. We use the basis space's ring instead of the
666 # element's ring because the basis space might be an algebraic
667 # closure whereas the base ring of the 3-by-3 identity matrix
668 # could be QQ instead of QQbar.
670 # And, we also have to handle Cartesian product bases (when
671 # the matrix basis consists of tuples) here. The "good news"
672 # is that we're already converting everything to long vectors,
673 # and that strategy works for tuples as well.
675 # We pass check=False because the matrix basis is "guaranteed"
676 # to be linearly independent... right? Ha ha.
678 V
= VectorSpace(self
.base_ring(), len(elt
))
679 W
= V
.span_of_basis( (V(_all2list(s
)) for s
in self
.matrix_basis()),
683 coords
= W
.coordinate_vector(V(elt
))
684 except ArithmeticError: # vector is not in free module
685 raise ValueError(msg
)
687 return self
.from_vector(coords
)
691 Return a string representation of ``self``.
695 sage: from mjo.eja.eja_algebra import JordanSpinEJA
699 Ensure that it says what we think it says::
701 sage: JordanSpinEJA(2, field=AA)
702 Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
703 sage: JordanSpinEJA(3, field=RDF)
704 Euclidean Jordan algebra of dimension 3 over Real Double Field
707 fmt
= "Euclidean Jordan algebra of dimension {} over {}"
708 return fmt
.format(self
.dimension(), self
.base_ring())
712 def characteristic_polynomial_of(self
):
714 Return the algebra's "characteristic polynomial of" function,
715 which is itself a multivariate polynomial that, when evaluated
716 at the coordinates of some algebra element, returns that
717 element's characteristic polynomial.
719 The resulting polynomial has `n+1` variables, where `n` is the
720 dimension of this algebra. The first `n` variables correspond to
721 the coordinates of an algebra element: when evaluated at the
722 coordinates of an algebra element with respect to a certain
723 basis, the result is a univariate polynomial (in the one
724 remaining variable ``t``), namely the characteristic polynomial
729 sage: from mjo.eja.eja_algebra import JordanSpinEJA, TrivialEJA
733 The characteristic polynomial in the spin algebra is given in
734 Alizadeh, Example 11.11::
736 sage: J = JordanSpinEJA(3)
737 sage: p = J.characteristic_polynomial_of(); p
738 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
739 sage: xvec = J.one().to_vector()
743 By definition, the characteristic polynomial is a monic
744 degree-zero polynomial in a rank-zero algebra. Note that
745 Cayley-Hamilton is indeed satisfied since the polynomial
746 ``1`` evaluates to the identity element of the algebra on
749 sage: J = TrivialEJA()
750 sage: J.characteristic_polynomial_of()
757 # The list of coefficient polynomials a_0, a_1, a_2, ..., a_(r-1).
758 a
= self
._charpoly
_coefficients
()
760 # We go to a bit of trouble here to reorder the
761 # indeterminates, so that it's easier to evaluate the
762 # characteristic polynomial at x's coordinates and get back
763 # something in terms of t, which is what we want.
764 S
= PolynomialRing(self
.base_ring(),'t')
768 S
= PolynomialRing(S
, R
.variable_names())
771 return (t
**r
+ sum( a
[k
]*(t
**k
) for k
in range(r
) ))
773 def coordinate_polynomial_ring(self
):
775 The multivariate polynomial ring in which this algebra's
776 :meth:`characteristic_polynomial_of` lives.
780 sage: from mjo.eja.eja_algebra import (HadamardEJA,
781 ....: RealSymmetricEJA)
785 sage: J = HadamardEJA(2)
786 sage: J.coordinate_polynomial_ring()
787 Multivariate Polynomial Ring in X1, X2...
788 sage: J = RealSymmetricEJA(3,field=QQ,orthonormalize=False)
789 sage: J.coordinate_polynomial_ring()
790 Multivariate Polynomial Ring in X1, X2, X3, X4, X5, X6...
793 var_names
= tuple( "X%d" % z
for z
in range(1, self
.dimension()+1) )
794 return PolynomialRing(self
.base_ring(), var_names
)
796 def inner_product(self
, x
, y
):
798 The inner product associated with this Euclidean Jordan algebra.
800 Defaults to the trace inner product, but can be overridden by
801 subclasses if they are sure that the necessary properties are
806 sage: from mjo.eja.eja_algebra import (random_eja,
808 ....: BilinearFormEJA)
812 Our inner product is "associative," which means the following for
813 a symmetric bilinear form::
815 sage: set_random_seed()
816 sage: J = random_eja()
817 sage: x,y,z = J.random_elements(3)
818 sage: (x*y).inner_product(z) == y.inner_product(x*z)
823 Ensure that this is the usual inner product for the algebras
826 sage: set_random_seed()
827 sage: J = HadamardEJA.random_instance()
828 sage: x,y = J.random_elements(2)
829 sage: actual = x.inner_product(y)
830 sage: expected = x.to_vector().inner_product(y.to_vector())
831 sage: actual == expected
834 Ensure that this is one-half of the trace inner-product in a
835 BilinearFormEJA that isn't just the reals (when ``n`` isn't
836 one). This is in Faraut and Koranyi, and also my "On the
839 sage: set_random_seed()
840 sage: J = BilinearFormEJA.random_instance()
841 sage: n = J.dimension()
842 sage: x = J.random_element()
843 sage: y = J.random_element()
844 sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
847 B
= self
._inner
_product
_matrix
848 return (B
*x
.to_vector()).inner_product(y
.to_vector())
851 def is_trivial(self
):
853 Return whether or not this algebra is trivial.
855 A trivial algebra contains only the zero element.
859 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
864 sage: J = ComplexHermitianEJA(3)
870 sage: J = TrivialEJA()
875 return self
.dimension() == 0
878 def multiplication_table(self
):
880 Return a visual representation of this algebra's multiplication
881 table (on basis elements).
885 sage: from mjo.eja.eja_algebra import JordanSpinEJA
889 sage: J = JordanSpinEJA(4)
890 sage: J.multiplication_table()
891 +----++----+----+----+----+
892 | * || b0 | b1 | b2 | b3 |
893 +====++====+====+====+====+
894 | b0 || b0 | b1 | b2 | b3 |
895 +----++----+----+----+----+
896 | b1 || b1 | b0 | 0 | 0 |
897 +----++----+----+----+----+
898 | b2 || b2 | 0 | b0 | 0 |
899 +----++----+----+----+----+
900 | b3 || b3 | 0 | 0 | b0 |
901 +----++----+----+----+----+
905 # Prepend the header row.
906 M
= [["*"] + list(self
.gens())]
908 # And to each subsequent row, prepend an entry that belongs to
909 # the left-side "header column."
910 M
+= [ [self
.monomial(i
)] + [ self
.monomial(i
)*self
.monomial(j
)
914 return table(M
, header_row
=True, header_column
=True, frame
=True)
917 def matrix_basis(self
):
919 Return an (often more natural) representation of this algebras
920 basis as an ordered tuple of matrices.
922 Every finite-dimensional Euclidean Jordan Algebra is a, up to
923 Jordan isomorphism, a direct sum of five simple
924 algebras---four of which comprise Hermitian matrices. And the
925 last type of algebra can of course be thought of as `n`-by-`1`
926 column matrices (ambiguusly called column vectors) to avoid
927 special cases. As a result, matrices (and column vectors) are
928 a natural representation format for Euclidean Jordan algebra
931 But, when we construct an algebra from a basis of matrices,
932 those matrix representations are lost in favor of coordinate
933 vectors *with respect to* that basis. We could eventually
934 convert back if we tried hard enough, but having the original
935 representations handy is valuable enough that we simply store
936 them and return them from this method.
938 Why implement this for non-matrix algebras? Avoiding special
939 cases for the :class:`BilinearFormEJA` pays with simplicity in
940 its own right. But mainly, we would like to be able to assume
941 that elements of a :class:`CartesianProductEJA` can be displayed
942 nicely, without having to have special classes for direct sums
943 one of whose components was a matrix algebra.
947 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
948 ....: RealSymmetricEJA)
952 sage: J = RealSymmetricEJA(2)
954 Finite family {0: b0, 1: b1, 2: b2}
955 sage: J.matrix_basis()
957 [1 0] [ 0 0.7071067811865475?] [0 0]
958 [0 0], [0.7071067811865475? 0], [0 1]
963 sage: J = JordanSpinEJA(2)
965 Finite family {0: b0, 1: b1}
966 sage: J.matrix_basis()
972 return self
._matrix
_basis
975 def matrix_space(self
):
977 Return the matrix space in which this algebra's elements live, if
978 we think of them as matrices (including column vectors of the
981 "By default" this will be an `n`-by-`1` column-matrix space,
982 except when the algebra is trivial. There it's `n`-by-`n`
983 (where `n` is zero), to ensure that two elements of the matrix
984 space (empty matrices) can be multiplied. For algebras of
985 matrices, this returns the space in which their
986 real embeddings live.
990 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
992 ....: QuaternionHermitianEJA,
997 By default, the matrix representation is just a column-matrix
998 equivalent to the vector representation::
1000 sage: J = JordanSpinEJA(3)
1001 sage: J.matrix_space()
1002 Full MatrixSpace of 3 by 1 dense matrices over Algebraic
1005 The matrix representation in the trivial algebra is
1006 zero-by-zero instead of the usual `n`-by-one::
1008 sage: J = TrivialEJA()
1009 sage: J.matrix_space()
1010 Full MatrixSpace of 0 by 0 dense matrices over Algebraic
1013 The matrix space for complex/quaternion Hermitian matrix EJA
1014 is the space in which their real-embeddings live, not the
1015 original complex/quaternion matrix space::
1017 sage: J = ComplexHermitianEJA(2,field=QQ,orthonormalize=False)
1018 sage: J.matrix_space()
1019 Full MatrixSpace of 4 by 4 dense matrices over Rational Field
1020 sage: J = QuaternionHermitianEJA(1,field=QQ,orthonormalize=False)
1021 sage: J.matrix_space()
1022 Full MatrixSpace of 4 by 4 dense matrices over Rational Field
1025 if self
.is_trivial():
1026 return MatrixSpace(self
.base_ring(), 0)
1028 return self
.matrix_basis()[0].parent()
1034 Return the unit element of this algebra.
1038 sage: from mjo.eja.eja_algebra import (HadamardEJA,
1043 We can compute unit element in the Hadamard EJA::
1045 sage: J = HadamardEJA(5)
1047 b0 + b1 + b2 + b3 + b4
1049 The unit element in the Hadamard EJA is inherited in the
1050 subalgebras generated by its elements::
1052 sage: J = HadamardEJA(5)
1054 b0 + b1 + b2 + b3 + b4
1055 sage: x = sum(J.gens())
1056 sage: A = x.subalgebra_generated_by(orthonormalize=False)
1059 sage: A.one().superalgebra_element()
1060 b0 + b1 + b2 + b3 + b4
1064 The identity element acts like the identity, regardless of
1065 whether or not we orthonormalize::
1067 sage: set_random_seed()
1068 sage: J = random_eja()
1069 sage: x = J.random_element()
1070 sage: J.one()*x == x and x*J.one() == x
1072 sage: A = x.subalgebra_generated_by()
1073 sage: y = A.random_element()
1074 sage: A.one()*y == y and y*A.one() == y
1079 sage: set_random_seed()
1080 sage: J = random_eja(field=QQ, orthonormalize=False)
1081 sage: x = J.random_element()
1082 sage: J.one()*x == x and x*J.one() == x
1084 sage: A = x.subalgebra_generated_by(orthonormalize=False)
1085 sage: y = A.random_element()
1086 sage: A.one()*y == y and y*A.one() == y
1089 The matrix of the unit element's operator is the identity,
1090 regardless of the base field and whether or not we
1093 sage: set_random_seed()
1094 sage: J = random_eja()
1095 sage: actual = J.one().operator().matrix()
1096 sage: expected = matrix.identity(J.base_ring(), J.dimension())
1097 sage: actual == expected
1099 sage: x = J.random_element()
1100 sage: A = x.subalgebra_generated_by()
1101 sage: actual = A.one().operator().matrix()
1102 sage: expected = matrix.identity(A.base_ring(), A.dimension())
1103 sage: actual == expected
1108 sage: set_random_seed()
1109 sage: J = random_eja(field=QQ, orthonormalize=False)
1110 sage: actual = J.one().operator().matrix()
1111 sage: expected = matrix.identity(J.base_ring(), J.dimension())
1112 sage: actual == expected
1114 sage: x = J.random_element()
1115 sage: A = x.subalgebra_generated_by(orthonormalize=False)
1116 sage: actual = A.one().operator().matrix()
1117 sage: expected = matrix.identity(A.base_ring(), A.dimension())
1118 sage: actual == expected
1121 Ensure that the cached unit element (often precomputed by
1122 hand) agrees with the computed one::
1124 sage: set_random_seed()
1125 sage: J = random_eja()
1126 sage: cached = J.one()
1127 sage: J.one.clear_cache()
1128 sage: J.one() == cached
1133 sage: set_random_seed()
1134 sage: J = random_eja(field=QQ, orthonormalize=False)
1135 sage: cached = J.one()
1136 sage: J.one.clear_cache()
1137 sage: J.one() == cached
1141 # We can brute-force compute the matrices of the operators
1142 # that correspond to the basis elements of this algebra.
1143 # If some linear combination of those basis elements is the
1144 # algebra identity, then the same linear combination of
1145 # their matrices has to be the identity matrix.
1147 # Of course, matrices aren't vectors in sage, so we have to
1148 # appeal to the "long vectors" isometry.
1149 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
1151 # Now we use basic linear algebra to find the coefficients,
1152 # of the matrices-as-vectors-linear-combination, which should
1153 # work for the original algebra basis too.
1154 A
= matrix(self
.base_ring(), oper_vecs
)
1156 # We used the isometry on the left-hand side already, but we
1157 # still need to do it for the right-hand side. Recall that we
1158 # wanted something that summed to the identity matrix.
1159 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
1161 # Now if there's an identity element in the algebra, this
1162 # should work. We solve on the left to avoid having to
1163 # transpose the matrix "A".
1164 return self
.from_vector(A
.solve_left(b
))
1167 def peirce_decomposition(self
, c
):
1169 The Peirce decomposition of this algebra relative to the
1172 In the future, this can be extended to a complete system of
1173 orthogonal idempotents.
1177 - ``c`` -- an idempotent of this algebra.
1181 A triple (J0, J5, J1) containing two subalgebras and one subspace
1184 - ``J0`` -- the algebra on the eigenspace of ``c.operator()``
1185 corresponding to the eigenvalue zero.
1187 - ``J5`` -- the eigenspace (NOT a subalgebra) of ``c.operator()``
1188 corresponding to the eigenvalue one-half.
1190 - ``J1`` -- the algebra on the eigenspace of ``c.operator()``
1191 corresponding to the eigenvalue one.
1193 These are the only possible eigenspaces for that operator, and this
1194 algebra is a direct sum of them. The spaces ``J0`` and ``J1`` are
1195 orthogonal, and are subalgebras of this algebra with the appropriate
1200 sage: from mjo.eja.eja_algebra import random_eja, RealSymmetricEJA
1204 The canonical example comes from the symmetric matrices, which
1205 decompose into diagonal and off-diagonal parts::
1207 sage: J = RealSymmetricEJA(3)
1208 sage: C = matrix(QQ, [ [1,0,0],
1212 sage: J0,J5,J1 = J.peirce_decomposition(c)
1214 Euclidean Jordan algebra of dimension 1...
1216 Vector space of degree 6 and dimension 2...
1218 Euclidean Jordan algebra of dimension 3...
1219 sage: J0.one().to_matrix()
1223 sage: orig_df = AA.options.display_format
1224 sage: AA.options.display_format = 'radical'
1225 sage: J.from_vector(J5.basis()[0]).to_matrix()
1229 sage: J.from_vector(J5.basis()[1]).to_matrix()
1233 sage: AA.options.display_format = orig_df
1234 sage: J1.one().to_matrix()
1241 Every algebra decomposes trivially with respect to its identity
1244 sage: set_random_seed()
1245 sage: J = random_eja()
1246 sage: J0,J5,J1 = J.peirce_decomposition(J.one())
1247 sage: J0.dimension() == 0 and J5.dimension() == 0
1249 sage: J1.superalgebra() == J and J1.dimension() == J.dimension()
1252 The decomposition is into eigenspaces, and its components are
1253 therefore necessarily orthogonal. Moreover, the identity
1254 elements in the two subalgebras are the projections onto their
1255 respective subspaces of the superalgebra's identity element::
1257 sage: set_random_seed()
1258 sage: J = random_eja()
1259 sage: x = J.random_element()
1260 sage: if not J.is_trivial():
1261 ....: while x.is_nilpotent():
1262 ....: x = J.random_element()
1263 sage: c = x.subalgebra_idempotent()
1264 sage: J0,J5,J1 = J.peirce_decomposition(c)
1266 sage: for (w,y,z) in zip(J0.basis(), J5.basis(), J1.basis()):
1267 ....: w = w.superalgebra_element()
1268 ....: y = J.from_vector(y)
1269 ....: z = z.superalgebra_element()
1270 ....: ipsum += w.inner_product(y).abs()
1271 ....: ipsum += w.inner_product(z).abs()
1272 ....: ipsum += y.inner_product(z).abs()
1275 sage: J1(c) == J1.one()
1277 sage: J0(J.one() - c) == J0.one()
1281 if not c
.is_idempotent():
1282 raise ValueError("element is not idempotent: %s" % c
)
1284 # Default these to what they should be if they turn out to be
1285 # trivial, because eigenspaces_left() won't return eigenvalues
1286 # corresponding to trivial spaces (e.g. it returns only the
1287 # eigenspace corresponding to lambda=1 if you take the
1288 # decomposition relative to the identity element).
1289 trivial
= self
.subalgebra(())
1290 J0
= trivial
# eigenvalue zero
1291 J5
= VectorSpace(self
.base_ring(), 0) # eigenvalue one-half
1292 J1
= trivial
# eigenvalue one
1294 for (eigval
, eigspace
) in c
.operator().matrix().right_eigenspaces():
1295 if eigval
== ~
(self
.base_ring()(2)):
1298 gens
= tuple( self
.from_vector(b
) for b
in eigspace
.basis() )
1299 subalg
= self
.subalgebra(gens
, check_axioms
=False)
1305 raise ValueError("unexpected eigenvalue: %s" % eigval
)
1310 def random_element(self
, thorough
=False):
1312 Return a random element of this algebra.
1314 Our algebra superclass method only returns a linear
1315 combination of at most two basis elements. We instead
1316 want the vector space "random element" method that
1317 returns a more diverse selection.
1321 - ``thorough`` -- (boolean; default False) whether or not we
1322 should generate irrational coefficients for the random
1323 element when our base ring is irrational; this slows the
1324 algebra operations to a crawl, but any truly random method
1328 # For a general base ring... maybe we can trust this to do the
1329 # right thing? Unlikely, but.
1330 V
= self
.vector_space()
1331 v
= V
.random_element()
1333 if self
.base_ring() is AA
:
1334 # The "random element" method of the algebraic reals is
1335 # stupid at the moment, and only returns integers between
1336 # -2 and 2, inclusive:
1338 # https://trac.sagemath.org/ticket/30875
1340 # Instead, we implement our own "random vector" method,
1341 # and then coerce that into the algebra. We use the vector
1342 # space degree here instead of the dimension because a
1343 # subalgebra could (for example) be spanned by only two
1344 # vectors, each with five coordinates. We need to
1345 # generate all five coordinates.
1347 v
*= QQbar
.random_element().real()
1349 v
*= QQ
.random_element()
1351 return self
.from_vector(V
.coordinate_vector(v
))
1353 def random_elements(self
, count
, thorough
=False):
1355 Return ``count`` random elements as a tuple.
1359 - ``thorough`` -- (boolean; default False) whether or not we
1360 should generate irrational coefficients for the random
1361 elements when our base ring is irrational; this slows the
1362 algebra operations to a crawl, but any truly random method
1367 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1371 sage: J = JordanSpinEJA(3)
1372 sage: x,y,z = J.random_elements(3)
1373 sage: all( [ x in J, y in J, z in J ])
1375 sage: len( J.random_elements(10) ) == 10
1379 return tuple( self
.random_element(thorough
)
1380 for idx
in range(count
) )
1384 def _charpoly_coefficients(self
):
1386 The `r` polynomial coefficients of the "characteristic polynomial
1391 sage: from mjo.eja.eja_algebra import random_eja
1395 The theory shows that these are all homogeneous polynomials of
1398 sage: set_random_seed()
1399 sage: J = random_eja()
1400 sage: all(p.is_homogeneous() for p in J._charpoly_coefficients())
1404 n
= self
.dimension()
1405 R
= self
.coordinate_polynomial_ring()
1407 F
= R
.fraction_field()
1410 # From a result in my book, these are the entries of the
1411 # basis representation of L_x.
1412 return sum( vars[k
]*self
.monomial(k
).operator().matrix()[i
,j
]
1415 L_x
= matrix(F
, n
, n
, L_x_i_j
)
1418 if self
.rank
.is_in_cache():
1420 # There's no need to pad the system with redundant
1421 # columns if we *know* they'll be redundant.
1424 # Compute an extra power in case the rank is equal to
1425 # the dimension (otherwise, we would stop at x^(r-1)).
1426 x_powers
= [ (L_x
**k
)*self
.one().to_vector()
1427 for k
in range(n
+1) ]
1428 A
= matrix
.column(F
, x_powers
[:n
])
1429 AE
= A
.extended_echelon_form()
1436 # The theory says that only the first "r" coefficients are
1437 # nonzero, and they actually live in the original polynomial
1438 # ring and not the fraction field. We negate them because in
1439 # the actual characteristic polynomial, they get moved to the
1440 # other side where x^r lives. We don't bother to trim A_rref
1441 # down to a square matrix and solve the resulting system,
1442 # because the upper-left r-by-r portion of A_rref is
1443 # guaranteed to be the identity matrix, so e.g.
1445 # A_rref.solve_right(Y)
1447 # would just be returning Y.
1448 return (-E
*b
)[:r
].change_ring(R
)
1453 Return the rank of this EJA.
1455 This is a cached method because we know the rank a priori for
1456 all of the algebras we can construct. Thus we can avoid the
1457 expensive ``_charpoly_coefficients()`` call unless we truly
1458 need to compute the whole characteristic polynomial.
1462 sage: from mjo.eja.eja_algebra import (HadamardEJA,
1463 ....: JordanSpinEJA,
1464 ....: RealSymmetricEJA,
1465 ....: ComplexHermitianEJA,
1466 ....: QuaternionHermitianEJA,
1471 The rank of the Jordan spin algebra is always two::
1473 sage: JordanSpinEJA(2).rank()
1475 sage: JordanSpinEJA(3).rank()
1477 sage: JordanSpinEJA(4).rank()
1480 The rank of the `n`-by-`n` Hermitian real, complex, or
1481 quaternion matrices is `n`::
1483 sage: RealSymmetricEJA(4).rank()
1485 sage: ComplexHermitianEJA(3).rank()
1487 sage: QuaternionHermitianEJA(2).rank()
1492 Ensure that every EJA that we know how to construct has a
1493 positive integer rank, unless the algebra is trivial in
1494 which case its rank will be zero::
1496 sage: set_random_seed()
1497 sage: J = random_eja()
1501 sage: r > 0 or (r == 0 and J.is_trivial())
1504 Ensure that computing the rank actually works, since the ranks
1505 of all simple algebras are known and will be cached by default::
1507 sage: set_random_seed() # long time
1508 sage: J = random_eja() # long time
1509 sage: cached = J.rank() # long time
1510 sage: J.rank.clear_cache() # long time
1511 sage: J.rank() == cached # long time
1515 return len(self
._charpoly
_coefficients
())
1518 def subalgebra(self
, basis
, **kwargs
):
1520 Create a subalgebra of this algebra from the given basis.
1522 from mjo
.eja
.eja_subalgebra
import FiniteDimensionalEJASubalgebra
1523 return FiniteDimensionalEJASubalgebra(self
, basis
, **kwargs
)
1526 def vector_space(self
):
1528 Return the vector space that underlies this algebra.
1532 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1536 sage: J = RealSymmetricEJA(2)
1537 sage: J.vector_space()
1538 Vector space of dimension 3 over...
1541 return self
.zero().to_vector().parent().ambient_vector_space()
1545 class RationalBasisEJA(FiniteDimensionalEJA
):
1547 Algebras whose supplied basis elements have all rational entries.
1551 sage: from mjo.eja.eja_algebra import BilinearFormEJA
1555 The supplied basis is orthonormalized by default::
1557 sage: B = matrix(QQ, [[1, 0, 0], [0, 25, -32], [0, -32, 41]])
1558 sage: J = BilinearFormEJA(B)
1559 sage: J.matrix_basis()
1576 # Abuse the check_field parameter to check that the entries of
1577 # out basis (in ambient coordinates) are in the field QQ.
1578 # Use _all2list to get the vector coordinates of octonion
1579 # entries and not the octonions themselves (which are not
1581 if not all( all(b_i
in QQ
for b_i
in _all2list(b
))
1583 raise TypeError("basis not rational")
1585 super().__init
__(basis
,
1589 check_field
=check_field
,
1592 self
._rational
_algebra
= None
1594 # There's no point in constructing the extra algebra if this
1595 # one is already rational.
1597 # Note: the same Jordan and inner-products work here,
1598 # because they are necessarily defined with respect to
1599 # ambient coordinates and not any particular basis.
1600 self
._rational
_algebra
= FiniteDimensionalEJA(
1605 associative
=self
.is_associative(),
1606 orthonormalize
=False,
1611 def _charpoly_coefficients(self
):
1615 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
1616 ....: JordanSpinEJA)
1620 The base ring of the resulting polynomial coefficients is what
1621 it should be, and not the rationals (unless the algebra was
1622 already over the rationals)::
1624 sage: J = JordanSpinEJA(3)
1625 sage: J._charpoly_coefficients()
1626 (X1^2 - X2^2 - X3^2, -2*X1)
1627 sage: a0 = J._charpoly_coefficients()[0]
1629 Algebraic Real Field
1630 sage: a0.base_ring()
1631 Algebraic Real Field
1634 if self
._rational
_algebra
is None:
1635 # There's no need to construct *another* algebra over the
1636 # rationals if this one is already over the
1637 # rationals. Likewise, if we never orthonormalized our
1638 # basis, we might as well just use the given one.
1639 return super()._charpoly
_coefficients
()
1641 # Do the computation over the rationals. The answer will be
1642 # the same, because all we've done is a change of basis.
1643 # Then, change back from QQ to our real base ring
1644 a
= ( a_i
.change_ring(self
.base_ring())
1645 for a_i
in self
._rational
_algebra
._charpoly
_coefficients
() )
1647 if self
._deortho
_matrix
is None:
1648 # This can happen if our base ring was, say, AA and we
1649 # chose not to (or didn't need to) orthonormalize. It's
1650 # still faster to do the computations over QQ even if
1651 # the numbers in the boxes stay the same.
1654 # Otherwise, convert the coordinate variables back to the
1655 # deorthonormalized ones.
1656 R
= self
.coordinate_polynomial_ring()
1657 from sage
.modules
.free_module_element
import vector
1658 X
= vector(R
, R
.gens())
1659 BX
= self
._deortho
_matrix
*X
1661 subs_dict
= { X[i]: BX[i] for i in range(len(X)) }
1662 return tuple( a_i
.subs(subs_dict
) for a_i
in a
)
1664 class ConcreteEJA(FiniteDimensionalEJA
):
1666 A class for the Euclidean Jordan algebras that we know by name.
1668 These are the Jordan algebras whose basis, multiplication table,
1669 rank, and so on are known a priori. More to the point, they are
1670 the Euclidean Jordan algebras for which we are able to conjure up
1671 a "random instance."
1675 sage: from mjo.eja.eja_algebra import ConcreteEJA
1679 Our basis is normalized with respect to the algebra's inner
1680 product, unless we specify otherwise::
1682 sage: set_random_seed()
1683 sage: J = ConcreteEJA.random_instance()
1684 sage: all( b.norm() == 1 for b in J.gens() )
1687 Since our basis is orthonormal with respect to the algebra's inner
1688 product, and since we know that this algebra is an EJA, any
1689 left-multiplication operator's matrix will be symmetric because
1690 natural->EJA basis representation is an isometry and within the
1691 EJA the operator is self-adjoint by the Jordan axiom::
1693 sage: set_random_seed()
1694 sage: J = ConcreteEJA.random_instance()
1695 sage: x = J.random_element()
1696 sage: x.operator().is_self_adjoint()
1701 def _max_random_instance_size():
1703 Return an integer "size" that is an upper bound on the size of
1704 this algebra when it is used in a random test
1705 case. Unfortunately, the term "size" is ambiguous -- when
1706 dealing with `R^n` under either the Hadamard or Jordan spin
1707 product, the "size" refers to the dimension `n`. When dealing
1708 with a matrix algebra (real symmetric or complex/quaternion
1709 Hermitian), it refers to the size of the matrix, which is far
1710 less than the dimension of the underlying vector space.
1712 This method must be implemented in each subclass.
1714 raise NotImplementedError
1717 def random_instance(cls
, *args
, **kwargs
):
1719 Return a random instance of this type of algebra.
1721 This method should be implemented in each subclass.
1723 from sage
.misc
.prandom
import choice
1724 eja_class
= choice(cls
.__subclasses
__())
1726 # These all bubble up to the RationalBasisEJA superclass
1727 # constructor, so any (kw)args valid there are also valid
1729 return eja_class
.random_instance(*args
, **kwargs
)
1734 def jordan_product(X
,Y
):
1735 return (X
*Y
+ Y
*X
)/2
1738 def trace_inner_product(X
,Y
):
1740 A trace inner-product for matrices that aren't embedded in the
1741 reals. It takes MATRICES as arguments, not EJA elements.
1743 return (X
*Y
).trace().real()
1745 class RealEmbeddedMatrixEJA(MatrixEJA
):
1747 def dimension_over_reals():
1749 The dimension of this matrix's base ring over the reals.
1751 The reals are dimension one over themselves, obviously; that's
1752 just `\mathbb{R}^{1}`. Likewise, the complex numbers `a + bi`
1753 have dimension two. Finally, the quaternions have dimension
1754 four over the reals.
1756 This is used to determine the size of the matrix returned from
1757 :meth:`real_embed`, among other things.
1759 raise NotImplementedError
1762 def real_embed(cls
,M
):
1764 Embed the matrix ``M`` into a space of real matrices.
1766 The matrix ``M`` can have entries in any field at the moment:
1767 the real numbers, complex numbers, or quaternions. And although
1768 they are not a field, we can probably support octonions at some
1769 point, too. This function returns a real matrix that "acts like"
1770 the original with respect to matrix multiplication; i.e.
1772 real_embed(M*N) = real_embed(M)*real_embed(N)
1775 if M
.ncols() != M
.nrows():
1776 raise ValueError("the matrix 'M' must be square")
1781 def real_unembed(cls
,M
):
1783 The inverse of :meth:`real_embed`.
1785 if M
.ncols() != M
.nrows():
1786 raise ValueError("the matrix 'M' must be square")
1787 if not ZZ(M
.nrows()).mod(cls
.dimension_over_reals()).is_zero():
1788 raise ValueError("the matrix 'M' must be a real embedding")
1793 def trace_inner_product(cls
,X
,Y
):
1795 Compute the trace inner-product of two real-embeddings.
1799 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
1800 ....: QuaternionHermitianEJA)
1804 sage: set_random_seed()
1805 sage: J = ComplexHermitianEJA.random_instance()
1806 sage: x,y = J.random_elements(2)
1807 sage: Xe = x.to_matrix()
1808 sage: Ye = y.to_matrix()
1809 sage: X = J.real_unembed(Xe)
1810 sage: Y = J.real_unembed(Ye)
1811 sage: expected = (X*Y).trace().real()
1812 sage: actual = J.trace_inner_product(Xe,Ye)
1813 sage: actual == expected
1818 sage: set_random_seed()
1819 sage: J = QuaternionHermitianEJA.random_instance()
1820 sage: x,y = J.random_elements(2)
1821 sage: Xe = x.to_matrix()
1822 sage: Ye = y.to_matrix()
1823 sage: X = J.real_unembed(Xe)
1824 sage: Y = J.real_unembed(Ye)
1825 sage: expected = (X*Y).trace().coefficient_tuple()[0]
1826 sage: actual = J.trace_inner_product(Xe,Ye)
1827 sage: actual == expected
1831 # This does in fact compute the real part of the trace.
1832 # If we compute the trace of e.g. a complex matrix M,
1833 # then we do so by adding up its diagonal entries --
1834 # call them z_1 through z_n. The real embedding of z_1
1835 # will be a 2-by-2 REAL matrix [a, b; -b, a] whose trace
1836 # as a REAL matrix will be 2*a = 2*Re(z_1). And so forth.
1837 return (X
*Y
).trace()/cls
.dimension_over_reals()
1839 class RealSymmetricEJA(RationalBasisEJA
, ConcreteEJA
, MatrixEJA
):
1841 The rank-n simple EJA consisting of real symmetric n-by-n
1842 matrices, the usual symmetric Jordan product, and the trace inner
1843 product. It has dimension `(n^2 + n)/2` over the reals.
1847 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1851 sage: J = RealSymmetricEJA(2)
1852 sage: b0, b1, b2 = J.gens()
1860 In theory, our "field" can be any subfield of the reals::
1862 sage: RealSymmetricEJA(2, field=RDF, check_axioms=True)
1863 Euclidean Jordan algebra of dimension 3 over Real Double Field
1864 sage: RealSymmetricEJA(2, field=RR, check_axioms=True)
1865 Euclidean Jordan algebra of dimension 3 over Real Field with
1866 53 bits of precision
1870 The dimension of this algebra is `(n^2 + n) / 2`::
1872 sage: set_random_seed()
1873 sage: n_max = RealSymmetricEJA._max_random_instance_size()
1874 sage: n = ZZ.random_element(1, n_max)
1875 sage: J = RealSymmetricEJA(n)
1876 sage: J.dimension() == (n^2 + n)/2
1879 The Jordan multiplication is what we think it is::
1881 sage: set_random_seed()
1882 sage: J = RealSymmetricEJA.random_instance()
1883 sage: x,y = J.random_elements(2)
1884 sage: actual = (x*y).to_matrix()
1885 sage: X = x.to_matrix()
1886 sage: Y = y.to_matrix()
1887 sage: expected = (X*Y + Y*X)/2
1888 sage: actual == expected
1890 sage: J(expected) == x*y
1893 We can change the generator prefix::
1895 sage: RealSymmetricEJA(3, prefix='q').gens()
1896 (q0, q1, q2, q3, q4, q5)
1898 We can construct the (trivial) algebra of rank zero::
1900 sage: RealSymmetricEJA(0)
1901 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1905 def _denormalized_basis(cls
, n
, field
):
1907 Return a basis for the space of real symmetric n-by-n matrices.
1911 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1915 sage: set_random_seed()
1916 sage: n = ZZ.random_element(1,5)
1917 sage: B = RealSymmetricEJA._denormalized_basis(n,ZZ)
1918 sage: all( M.is_symmetric() for M in B)
1922 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1926 for j
in range(i
+1):
1927 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
1931 Sij
= Eij
+ Eij
.transpose()
1937 def _max_random_instance_size():
1938 return 4 # Dimension 10
1941 def random_instance(cls
, **kwargs
):
1943 Return a random instance of this type of algebra.
1945 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1946 return cls(n
, **kwargs
)
1948 def __init__(self
, n
, field
=AA
, **kwargs
):
1949 # We know this is a valid EJA, but will double-check
1950 # if the user passes check_axioms=True.
1951 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
1957 super().__init
__(self
._denormalized
_basis
(n
,field
),
1958 self
.jordan_product
,
1959 self
.trace_inner_product
,
1961 associative
=associative
,
1964 # TODO: this could be factored out somehow, but is left here
1965 # because the MatrixEJA is not presently a subclass of the
1966 # FDEJA class that defines rank() and one().
1967 self
.rank
.set_cache(n
)
1968 idV
= self
.matrix_space().one()
1969 self
.one
.set_cache(self(idV
))
1973 class ComplexMatrixEJA(RealEmbeddedMatrixEJA
):
1974 # A manual dictionary-cache for the complex_extension() method,
1975 # since apparently @classmethods can't also be @cached_methods.
1976 _complex_extension
= {}
1979 def complex_extension(cls
,field
):
1981 The complex field that we embed/unembed, as an extension
1982 of the given ``field``.
1984 if field
in cls
._complex
_extension
:
1985 return cls
._complex
_extension
[field
]
1987 # Sage doesn't know how to adjoin the complex "i" (the root of
1988 # x^2 + 1) to a field in a general way. Here, we just enumerate
1989 # all of the cases that I have cared to support so far.
1991 # Sage doesn't know how to embed AA into QQbar, i.e. how
1992 # to adjoin sqrt(-1) to AA.
1994 elif not field
.is_exact():
1996 F
= field
.complex_field()
1998 # Works for QQ and... maybe some other fields.
1999 R
= PolynomialRing(field
, 'z')
2001 F
= field
.extension(z
**2 + 1, 'I', embedding
=CLF(-1).sqrt())
2003 cls
._complex
_extension
[field
] = F
2007 def dimension_over_reals():
2011 def real_embed(cls
,M
):
2013 Embed the n-by-n complex matrix ``M`` into the space of real
2014 matrices of size 2n-by-2n via the map the sends each entry `z = a +
2015 bi` to the block matrix ``[[a,b],[-b,a]]``.
2019 sage: from mjo.eja.eja_algebra import ComplexMatrixEJA
2023 sage: F = QuadraticField(-1, 'I')
2024 sage: x1 = F(4 - 2*i)
2025 sage: x2 = F(1 + 2*i)
2028 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
2029 sage: ComplexMatrixEJA.real_embed(M)
2038 Embedding is a homomorphism (isomorphism, in fact)::
2040 sage: set_random_seed()
2041 sage: n = ZZ.random_element(3)
2042 sage: F = QuadraticField(-1, 'I')
2043 sage: X = random_matrix(F, n)
2044 sage: Y = random_matrix(F, n)
2045 sage: Xe = ComplexMatrixEJA.real_embed(X)
2046 sage: Ye = ComplexMatrixEJA.real_embed(Y)
2047 sage: XYe = ComplexMatrixEJA.real_embed(X*Y)
2052 super().real_embed(M
)
2055 # We don't need any adjoined elements...
2056 field
= M
.base_ring().base_ring()
2062 blocks
.append(matrix(field
, 2, [ [ a
, b
],
2065 return matrix
.block(field
, n
, blocks
)
2069 def real_unembed(cls
,M
):
2071 The inverse of _embed_complex_matrix().
2075 sage: from mjo.eja.eja_algebra import ComplexMatrixEJA
2079 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
2080 ....: [-2, 1, -4, 3],
2081 ....: [ 9, 10, 11, 12],
2082 ....: [-10, 9, -12, 11] ])
2083 sage: ComplexMatrixEJA.real_unembed(A)
2085 [ 10*I + 9 12*I + 11]
2089 Unembedding is the inverse of embedding::
2091 sage: set_random_seed()
2092 sage: F = QuadraticField(-1, 'I')
2093 sage: M = random_matrix(F, 3)
2094 sage: Me = ComplexMatrixEJA.real_embed(M)
2095 sage: ComplexMatrixEJA.real_unembed(Me) == M
2099 super().real_unembed(M
)
2101 d
= cls
.dimension_over_reals()
2102 F
= cls
.complex_extension(M
.base_ring())
2105 # Go top-left to bottom-right (reading order), converting every
2106 # 2-by-2 block we see to a single complex element.
2108 for k
in range(n
/d
):
2109 for j
in range(n
/d
):
2110 submat
= M
[d
*k
:d
*k
+d
,d
*j
:d
*j
+d
]
2111 if submat
[0,0] != submat
[1,1]:
2112 raise ValueError('bad on-diagonal submatrix')
2113 if submat
[0,1] != -submat
[1,0]:
2114 raise ValueError('bad off-diagonal submatrix')
2115 z
= submat
[0,0] + submat
[0,1]*i
2118 return matrix(F
, n
/d
, elements
)
2121 class ComplexHermitianEJA(RationalBasisEJA
, ConcreteEJA
, ComplexMatrixEJA
):
2123 The rank-n simple EJA consisting of complex Hermitian n-by-n
2124 matrices over the real numbers, the usual symmetric Jordan product,
2125 and the real-part-of-trace inner product. It has dimension `n^2` over
2130 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
2134 In theory, our "field" can be any subfield of the reals::
2136 sage: ComplexHermitianEJA(2, field=RDF, check_axioms=True)
2137 Euclidean Jordan algebra of dimension 4 over Real Double Field
2138 sage: ComplexHermitianEJA(2, field=RR, check_axioms=True)
2139 Euclidean Jordan algebra of dimension 4 over Real Field with
2140 53 bits of precision
2144 The dimension of this algebra is `n^2`::
2146 sage: set_random_seed()
2147 sage: n_max = ComplexHermitianEJA._max_random_instance_size()
2148 sage: n = ZZ.random_element(1, n_max)
2149 sage: J = ComplexHermitianEJA(n)
2150 sage: J.dimension() == n^2
2153 The Jordan multiplication is what we think it is::
2155 sage: set_random_seed()
2156 sage: J = ComplexHermitianEJA.random_instance()
2157 sage: x,y = J.random_elements(2)
2158 sage: actual = (x*y).to_matrix()
2159 sage: X = x.to_matrix()
2160 sage: Y = y.to_matrix()
2161 sage: expected = (X*Y + Y*X)/2
2162 sage: actual == expected
2164 sage: J(expected) == x*y
2167 We can change the generator prefix::
2169 sage: ComplexHermitianEJA(2, prefix='z').gens()
2172 We can construct the (trivial) algebra of rank zero::
2174 sage: ComplexHermitianEJA(0)
2175 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2180 def _denormalized_basis(cls
, n
, field
):
2182 Returns a basis for the space of complex Hermitian n-by-n matrices.
2184 Why do we embed these? Basically, because all of numerical linear
2185 algebra assumes that you're working with vectors consisting of `n`
2186 entries from a field and scalars from the same field. There's no way
2187 to tell SageMath that (for example) the vectors contain complex
2188 numbers, while the scalar field is real.
2192 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
2196 sage: set_random_seed()
2197 sage: n = ZZ.random_element(1,5)
2198 sage: B = ComplexHermitianEJA._denormalized_basis(n,ZZ)
2199 sage: all( M.is_symmetric() for M in B)
2203 R
= PolynomialRing(ZZ
, 'z')
2205 F
= ZZ
.extension(z
**2 + 1, 'I')
2208 # This is like the symmetric case, but we need to be careful:
2210 # * We want conjugate-symmetry, not just symmetry.
2211 # * The diagonal will (as a result) be real.
2214 Eij
= matrix
.zero(F
,n
)
2216 for j
in range(i
+1):
2220 Sij
= cls
.real_embed(Eij
)
2223 # The second one has a minus because it's conjugated.
2224 Eij
[j
,i
] = 1 # Eij = Eij + Eij.transpose()
2225 Sij_real
= cls
.real_embed(Eij
)
2227 # Eij = I*Eij - I*Eij.transpose()
2230 Sij_imag
= cls
.real_embed(Eij
)
2236 # Since we embedded the entries, we can drop back to the
2237 # desired real "field" instead of the extension "F".
2238 return tuple( s
.change_ring(field
) for s
in S
)
2241 def __init__(self
, n
, field
=AA
, **kwargs
):
2242 # We know this is a valid EJA, but will double-check
2243 # if the user passes check_axioms=True.
2244 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2250 super().__init
__(self
._denormalized
_basis
(n
,field
),
2251 self
.jordan_product
,
2252 self
.trace_inner_product
,
2254 associative
=associative
,
2256 # TODO: this could be factored out somehow, but is left here
2257 # because the MatrixEJA is not presently a subclass of the
2258 # FDEJA class that defines rank() and one().
2259 self
.rank
.set_cache(n
)
2260 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
2261 self
.one
.set_cache(self(idV
))
2264 def _max_random_instance_size():
2265 return 3 # Dimension 9
2268 def random_instance(cls
, **kwargs
):
2270 Return a random instance of this type of algebra.
2272 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2273 return cls(n
, **kwargs
)
2275 class QuaternionMatrixEJA(RealEmbeddedMatrixEJA
):
2277 # A manual dictionary-cache for the quaternion_extension() method,
2278 # since apparently @classmethods can't also be @cached_methods.
2279 _quaternion_extension
= {}
2282 def quaternion_extension(cls
,field
):
2284 The quaternion field that we embed/unembed, as an extension
2285 of the given ``field``.
2287 if field
in cls
._quaternion
_extension
:
2288 return cls
._quaternion
_extension
[field
]
2290 Q
= QuaternionAlgebra(field
,-1,-1)
2292 cls
._quaternion
_extension
[field
] = Q
2296 def dimension_over_reals():
2300 def real_embed(cls
,M
):
2302 Embed the n-by-n quaternion matrix ``M`` into the space of real
2303 matrices of size 4n-by-4n by first sending each quaternion entry `z
2304 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
2305 c+di],[-c + di, a-bi]]`, and then embedding those into a real
2310 sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA
2314 sage: Q = QuaternionAlgebra(QQ,-1,-1)
2315 sage: i,j,k = Q.gens()
2316 sage: x = 1 + 2*i + 3*j + 4*k
2317 sage: M = matrix(Q, 1, [[x]])
2318 sage: QuaternionMatrixEJA.real_embed(M)
2324 Embedding is a homomorphism (isomorphism, in fact)::
2326 sage: set_random_seed()
2327 sage: n = ZZ.random_element(2)
2328 sage: Q = QuaternionAlgebra(QQ,-1,-1)
2329 sage: X = random_matrix(Q, n)
2330 sage: Y = random_matrix(Q, n)
2331 sage: Xe = QuaternionMatrixEJA.real_embed(X)
2332 sage: Ye = QuaternionMatrixEJA.real_embed(Y)
2333 sage: XYe = QuaternionMatrixEJA.real_embed(X*Y)
2338 super().real_embed(M
)
2339 quaternions
= M
.base_ring()
2342 F
= QuadraticField(-1, 'I')
2347 t
= z
.coefficient_tuple()
2352 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
2353 [-c
+ d
*i
, a
- b
*i
]])
2354 realM
= ComplexMatrixEJA
.real_embed(cplxM
)
2355 blocks
.append(realM
)
2357 # We should have real entries by now, so use the realest field
2358 # we've got for the return value.
2359 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
2364 def real_unembed(cls
,M
):
2366 The inverse of _embed_quaternion_matrix().
2370 sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA
2374 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
2375 ....: [-2, 1, -4, 3],
2376 ....: [-3, 4, 1, -2],
2377 ....: [-4, -3, 2, 1]])
2378 sage: QuaternionMatrixEJA.real_unembed(M)
2379 [1 + 2*i + 3*j + 4*k]
2383 Unembedding is the inverse of embedding::
2385 sage: set_random_seed()
2386 sage: Q = QuaternionAlgebra(QQ, -1, -1)
2387 sage: M = random_matrix(Q, 3)
2388 sage: Me = QuaternionMatrixEJA.real_embed(M)
2389 sage: QuaternionMatrixEJA.real_unembed(Me) == M
2393 super().real_unembed(M
)
2395 d
= cls
.dimension_over_reals()
2397 # Use the base ring of the matrix to ensure that its entries can be
2398 # multiplied by elements of the quaternion algebra.
2399 Q
= cls
.quaternion_extension(M
.base_ring())
2402 # Go top-left to bottom-right (reading order), converting every
2403 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
2406 for l
in range(n
/d
):
2407 for m
in range(n
/d
):
2408 submat
= ComplexMatrixEJA
.real_unembed(
2409 M
[d
*l
:d
*l
+d
,d
*m
:d
*m
+d
] )
2410 if submat
[0,0] != submat
[1,1].conjugate():
2411 raise ValueError('bad on-diagonal submatrix')
2412 if submat
[0,1] != -submat
[1,0].conjugate():
2413 raise ValueError('bad off-diagonal submatrix')
2414 z
= submat
[0,0].real()
2415 z
+= submat
[0,0].imag()*i
2416 z
+= submat
[0,1].real()*j
2417 z
+= submat
[0,1].imag()*k
2420 return matrix(Q
, n
/d
, elements
)
2423 class QuaternionHermitianEJA(RationalBasisEJA
,
2425 QuaternionMatrixEJA
):
2427 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
2428 matrices, the usual symmetric Jordan product, and the
2429 real-part-of-trace inner product. It has dimension `2n^2 - n` over
2434 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2438 In theory, our "field" can be any subfield of the reals::
2440 sage: QuaternionHermitianEJA(2, field=RDF, check_axioms=True)
2441 Euclidean Jordan algebra of dimension 6 over Real Double Field
2442 sage: QuaternionHermitianEJA(2, field=RR, check_axioms=True)
2443 Euclidean Jordan algebra of dimension 6 over Real Field with
2444 53 bits of precision
2448 The dimension of this algebra is `2*n^2 - n`::
2450 sage: set_random_seed()
2451 sage: n_max = QuaternionHermitianEJA._max_random_instance_size()
2452 sage: n = ZZ.random_element(1, n_max)
2453 sage: J = QuaternionHermitianEJA(n)
2454 sage: J.dimension() == 2*(n^2) - n
2457 The Jordan multiplication is what we think it is::
2459 sage: set_random_seed()
2460 sage: J = QuaternionHermitianEJA.random_instance()
2461 sage: x,y = J.random_elements(2)
2462 sage: actual = (x*y).to_matrix()
2463 sage: X = x.to_matrix()
2464 sage: Y = y.to_matrix()
2465 sage: expected = (X*Y + Y*X)/2
2466 sage: actual == expected
2468 sage: J(expected) == x*y
2471 We can change the generator prefix::
2473 sage: QuaternionHermitianEJA(2, prefix='a').gens()
2474 (a0, a1, a2, a3, a4, a5)
2476 We can construct the (trivial) algebra of rank zero::
2478 sage: QuaternionHermitianEJA(0)
2479 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2483 def _denormalized_basis(cls
, n
, field
):
2485 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
2487 Why do we embed these? Basically, because all of numerical
2488 linear algebra assumes that you're working with vectors consisting
2489 of `n` entries from a field and scalars from the same field. There's
2490 no way to tell SageMath that (for example) the vectors contain
2491 complex numbers, while the scalar field is real.
2495 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2499 sage: set_random_seed()
2500 sage: n = ZZ.random_element(1,5)
2501 sage: B = QuaternionHermitianEJA._denormalized_basis(n,ZZ)
2502 sage: all( M.is_symmetric() for M in B )
2506 Q
= QuaternionAlgebra(QQ
,-1,-1)
2509 # This is like the symmetric case, but we need to be careful:
2511 # * We want conjugate-symmetry, not just symmetry.
2512 # * The diagonal will (as a result) be real.
2515 Eij
= matrix
.zero(Q
,n
)
2517 for j
in range(i
+1):
2521 Sij
= cls
.real_embed(Eij
)
2524 # The second, third, and fourth ones have a minus
2525 # because they're conjugated.
2526 # Eij = Eij + Eij.transpose()
2528 Sij_real
= cls
.real_embed(Eij
)
2530 # Eij = I*(Eij - Eij.transpose())
2533 Sij_I
= cls
.real_embed(Eij
)
2535 # Eij = J*(Eij - Eij.transpose())
2538 Sij_J
= cls
.real_embed(Eij
)
2540 # Eij = K*(Eij - Eij.transpose())
2543 Sij_K
= cls
.real_embed(Eij
)
2549 # Since we embedded the entries, we can drop back to the
2550 # desired real "field" instead of the quaternion algebra "Q".
2551 return tuple( s
.change_ring(field
) for s
in S
)
2554 def __init__(self
, n
, field
=AA
, **kwargs
):
2555 # We know this is a valid EJA, but will double-check
2556 # if the user passes check_axioms=True.
2557 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2563 super().__init
__(self
._denormalized
_basis
(n
,field
),
2564 self
.jordan_product
,
2565 self
.trace_inner_product
,
2567 associative
=associative
,
2570 # TODO: this could be factored out somehow, but is left here
2571 # because the MatrixEJA is not presently a subclass of the
2572 # FDEJA class that defines rank() and one().
2573 self
.rank
.set_cache(n
)
2574 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
2575 self
.one
.set_cache(self(idV
))
2579 def _max_random_instance_size():
2581 The maximum rank of a random QuaternionHermitianEJA.
2583 return 2 # Dimension 6
2586 def random_instance(cls
, **kwargs
):
2588 Return a random instance of this type of algebra.
2590 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2591 return cls(n
, **kwargs
)
2593 class OctonionHermitianEJA(RationalBasisEJA
, ConcreteEJA
, MatrixEJA
):
2597 sage: from mjo.eja.eja_algebra import (FiniteDimensionalEJA,
2598 ....: OctonionHermitianEJA)
2602 The 3-by-3 algebra satisfies the axioms of an EJA::
2604 sage: OctonionHermitianEJA(3, # long time
2605 ....: field=QQ, # long time
2606 ....: orthonormalize=False, # long time
2607 ....: check_axioms=True) # long time
2608 Euclidean Jordan algebra of dimension 27 over Rational Field
2610 After a change-of-basis, the 2-by-2 algebra has the same
2611 multiplication table as the ten-dimensional Jordan spin algebra::
2613 sage: b = OctonionHermitianEJA._denormalized_basis(2,QQ)
2614 sage: basis = (b[0] + b[9],) + b[1:9] + (b[0] - b[9],)
2615 sage: jp = OctonionHermitianEJA.jordan_product
2616 sage: ip = OctonionHermitianEJA.trace_inner_product
2617 sage: J = FiniteDimensionalEJA(basis,
2621 ....: orthonormalize=False)
2622 sage: J.multiplication_table()
2623 +----++----+----+----+----+----+----+----+----+----+----+
2624 | * || b0 | b1 | b2 | b3 | b4 | b5 | b6 | b7 | b8 | b9 |
2625 +====++====+====+====+====+====+====+====+====+====+====+
2626 | b0 || b0 | b1 | b2 | b3 | b4 | b5 | b6 | b7 | b8 | b9 |
2627 +----++----+----+----+----+----+----+----+----+----+----+
2628 | b1 || b1 | b0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
2629 +----++----+----+----+----+----+----+----+----+----+----+
2630 | b2 || b2 | 0 | b0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
2631 +----++----+----+----+----+----+----+----+----+----+----+
2632 | b3 || b3 | 0 | 0 | b0 | 0 | 0 | 0 | 0 | 0 | 0 |
2633 +----++----+----+----+----+----+----+----+----+----+----+
2634 | b4 || b4 | 0 | 0 | 0 | b0 | 0 | 0 | 0 | 0 | 0 |
2635 +----++----+----+----+----+----+----+----+----+----+----+
2636 | b5 || b5 | 0 | 0 | 0 | 0 | b0 | 0 | 0 | 0 | 0 |
2637 +----++----+----+----+----+----+----+----+----+----+----+
2638 | b6 || b6 | 0 | 0 | 0 | 0 | 0 | b0 | 0 | 0 | 0 |
2639 +----++----+----+----+----+----+----+----+----+----+----+
2640 | b7 || b7 | 0 | 0 | 0 | 0 | 0 | 0 | b0 | 0 | 0 |
2641 +----++----+----+----+----+----+----+----+----+----+----+
2642 | b8 || b8 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | b0 | 0 |
2643 +----++----+----+----+----+----+----+----+----+----+----+
2644 | b9 || b9 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | b0 |
2645 +----++----+----+----+----+----+----+----+----+----+----+
2649 We can actually construct the 27-dimensional Albert algebra,
2650 and we get the right unit element if we recompute it::
2652 sage: J = OctonionHermitianEJA(3, # long time
2653 ....: field=QQ, # long time
2654 ....: orthonormalize=False) # long time
2655 sage: J.one.clear_cache() # long time
2656 sage: J.one() # long time
2658 sage: J.one().to_matrix() # long time
2667 The 2-by-2 algebra is isomorphic to the ten-dimensional Jordan
2668 spin algebra, but just to be sure, we recompute its rank::
2670 sage: J = OctonionHermitianEJA(2, # long time
2671 ....: field=QQ, # long time
2672 ....: orthonormalize=False) # long time
2673 sage: J.rank.clear_cache() # long time
2674 sage: J.rank() # long time
2679 def _max_random_instance_size():
2681 The maximum rank of a random QuaternionHermitianEJA.
2683 return 1 # Dimension 1
2686 def random_instance(cls
, **kwargs
):
2688 Return a random instance of this type of algebra.
2690 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2691 return cls(n
, **kwargs
)
2693 def __init__(self
, n
, field
=AA
, **kwargs
):
2695 # Otherwise we don't get an EJA.
2696 raise ValueError("n cannot exceed 3")
2698 # We know this is a valid EJA, but will double-check
2699 # if the user passes check_axioms=True.
2700 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2702 super().__init
__(self
._denormalized
_basis
(n
,field
),
2703 self
.jordan_product
,
2704 self
.trace_inner_product
,
2708 # TODO: this could be factored out somehow, but is left here
2709 # because the MatrixEJA is not presently a subclass of the
2710 # FDEJA class that defines rank() and one().
2711 self
.rank
.set_cache(n
)
2712 idV
= self
.matrix_space().one()
2713 self
.one
.set_cache(self(idV
))
2717 def _denormalized_basis(cls
, n
, field
):
2719 Returns a basis for the space of octonion Hermitian n-by-n
2724 sage: from mjo.eja.eja_algebra import OctonionHermitianEJA
2728 sage: B = OctonionHermitianEJA._denormalized_basis(3,QQ)
2729 sage: all( M.is_hermitian() for M in B )
2735 from mjo
.octonions
import OctonionMatrixAlgebra
2736 MS
= OctonionMatrixAlgebra(n
, scalars
=field
)
2737 es
= MS
.entry_algebra().gens()
2741 for j
in range(i
+1):
2743 E_ii
= MS
.monomial( (i
,j
,es
[0]) )
2747 E_ij
= MS
.monomial( (i
,j
,e
) )
2749 # If the conjugate has a negative sign in front
2750 # of it, (j,i,ec) won't be a monomial!
2751 if (j
,i
,ec
) in MS
.indices():
2752 E_ij
+= MS
.monomial( (j
,i
,ec
) )
2754 E_ij
-= MS
.monomial( (j
,i
,-ec
) )
2757 return tuple( basis
)
2760 def trace_inner_product(X
,Y
):
2762 The octonions don't know that the reals are embedded in them,
2763 so we have to take the e0 component ourselves.
2767 sage: from mjo.eja.eja_algebra import OctonionHermitianEJA
2771 sage: J = OctonionHermitianEJA(2,field=QQ,orthonormalize=False)
2772 sage: I = J.one().to_matrix()
2773 sage: J.trace_inner_product(I, -I)
2777 return (X
*Y
).trace().real().coefficient(0)
2780 class AlbertEJA(OctonionHermitianEJA
):
2782 The Albert algebra is the algebra of three-by-three Hermitian
2783 matrices whose entries are octonions.
2787 sage: from mjo.eja.eja_algebra import AlbertEJA
2791 sage: AlbertEJA(field=QQ, orthonormalize=False)
2792 Euclidean Jordan algebra of dimension 27 over Rational Field
2793 sage: AlbertEJA() # long time
2794 Euclidean Jordan algebra of dimension 27 over Algebraic Real Field
2797 def __init__(self
, *args
, **kwargs
):
2798 super().__init
__(3, *args
, **kwargs
)
2801 class HadamardEJA(RationalBasisEJA
, ConcreteEJA
):
2803 Return the Euclidean Jordan algebra on `R^n` with the Hadamard
2804 (pointwise real-number multiplication) Jordan product and the
2805 usual inner-product.
2807 This is nothing more than the Cartesian product of ``n`` copies of
2808 the one-dimensional Jordan spin algebra, and is the most common
2809 example of a non-simple Euclidean Jordan algebra.
2813 sage: from mjo.eja.eja_algebra import HadamardEJA
2817 This multiplication table can be verified by hand::
2819 sage: J = HadamardEJA(3)
2820 sage: b0,b1,b2 = J.gens()
2836 We can change the generator prefix::
2838 sage: HadamardEJA(3, prefix='r').gens()
2841 def __init__(self
, n
, field
=AA
, **kwargs
):
2843 jordan_product
= lambda x
,y
: x
2844 inner_product
= lambda x
,y
: x
2846 def jordan_product(x
,y
):
2848 return P( xi
*yi
for (xi
,yi
) in zip(x
,y
) )
2850 def inner_product(x
,y
):
2853 # New defaults for keyword arguments. Don't orthonormalize
2854 # because our basis is already orthonormal with respect to our
2855 # inner-product. Don't check the axioms, because we know this
2856 # is a valid EJA... but do double-check if the user passes
2857 # check_axioms=True. Note: we DON'T override the "check_field"
2858 # default here, because the user can pass in a field!
2859 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2860 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2862 column_basis
= tuple( b
.column()
2863 for b
in FreeModule(field
, n
).basis() )
2864 super().__init
__(column_basis
,
2870 self
.rank
.set_cache(n
)
2873 self
.one
.set_cache( self
.zero() )
2875 self
.one
.set_cache( sum(self
.gens()) )
2878 def _max_random_instance_size():
2880 The maximum dimension of a random HadamardEJA.
2885 def random_instance(cls
, **kwargs
):
2887 Return a random instance of this type of algebra.
2889 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2890 return cls(n
, **kwargs
)
2893 class BilinearFormEJA(RationalBasisEJA
, ConcreteEJA
):
2895 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2896 with the half-trace inner product and jordan product ``x*y =
2897 (<Bx,y>,y_bar>, x0*y_bar + y0*x_bar)`` where `B = 1 \times B22` is
2898 a symmetric positive-definite "bilinear form" matrix. Its
2899 dimension is the size of `B`, and it has rank two in dimensions
2900 larger than two. It reduces to the ``JordanSpinEJA`` when `B` is
2901 the identity matrix of order ``n``.
2903 We insist that the one-by-one upper-left identity block of `B` be
2904 passed in as well so that we can be passed a matrix of size zero
2905 to construct a trivial algebra.
2909 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
2910 ....: JordanSpinEJA)
2914 When no bilinear form is specified, the identity matrix is used,
2915 and the resulting algebra is the Jordan spin algebra::
2917 sage: B = matrix.identity(AA,3)
2918 sage: J0 = BilinearFormEJA(B)
2919 sage: J1 = JordanSpinEJA(3)
2920 sage: J0.multiplication_table() == J0.multiplication_table()
2923 An error is raised if the matrix `B` does not correspond to a
2924 positive-definite bilinear form::
2926 sage: B = matrix.random(QQ,2,3)
2927 sage: J = BilinearFormEJA(B)
2928 Traceback (most recent call last):
2930 ValueError: bilinear form is not positive-definite
2931 sage: B = matrix.zero(QQ,3)
2932 sage: J = BilinearFormEJA(B)
2933 Traceback (most recent call last):
2935 ValueError: bilinear form is not positive-definite
2939 We can create a zero-dimensional algebra::
2941 sage: B = matrix.identity(AA,0)
2942 sage: J = BilinearFormEJA(B)
2946 We can check the multiplication condition given in the Jordan, von
2947 Neumann, and Wigner paper (and also discussed on my "On the
2948 symmetry..." paper). Note that this relies heavily on the standard
2949 choice of basis, as does anything utilizing the bilinear form
2950 matrix. We opt not to orthonormalize the basis, because if we
2951 did, we would have to normalize the `s_{i}` in a similar manner::
2953 sage: set_random_seed()
2954 sage: n = ZZ.random_element(5)
2955 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2956 sage: B11 = matrix.identity(QQ,1)
2957 sage: B22 = M.transpose()*M
2958 sage: B = block_matrix(2,2,[ [B11,0 ],
2960 sage: J = BilinearFormEJA(B, orthonormalize=False)
2961 sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
2962 sage: V = J.vector_space()
2963 sage: sis = [ J( V([0] + (M.inverse()*ei).list()).column() )
2964 ....: for ei in eis ]
2965 sage: actual = [ sis[i]*sis[j]
2966 ....: for i in range(n-1)
2967 ....: for j in range(n-1) ]
2968 sage: expected = [ J.one() if i == j else J.zero()
2969 ....: for i in range(n-1)
2970 ....: for j in range(n-1) ]
2971 sage: actual == expected
2975 def __init__(self
, B
, field
=AA
, **kwargs
):
2976 # The matrix "B" is supplied by the user in most cases,
2977 # so it makes sense to check whether or not its positive-
2978 # definite unless we are specifically asked not to...
2979 if ("check_axioms" not in kwargs
) or kwargs
["check_axioms"]:
2980 if not B
.is_positive_definite():
2981 raise ValueError("bilinear form is not positive-definite")
2983 # However, all of the other data for this EJA is computed
2984 # by us in manner that guarantees the axioms are
2985 # satisfied. So, again, unless we are specifically asked to
2986 # verify things, we'll skip the rest of the checks.
2987 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2989 def inner_product(x
,y
):
2990 return (y
.T
*B
*x
)[0,0]
2992 def jordan_product(x
,y
):
2998 z0
= inner_product(y
,x
)
2999 zbar
= y0
*xbar
+ x0
*ybar
3000 return P([z0
] + zbar
.list())
3003 column_basis
= tuple( b
.column()
3004 for b
in FreeModule(field
, n
).basis() )
3006 # TODO: I haven't actually checked this, but it seems legit.
3011 super().__init
__(column_basis
,
3015 associative
=associative
,
3018 # The rank of this algebra is two, unless we're in a
3019 # one-dimensional ambient space (because the rank is bounded
3020 # by the ambient dimension).
3021 self
.rank
.set_cache(min(n
,2))
3024 self
.one
.set_cache( self
.zero() )
3026 self
.one
.set_cache( self
.monomial(0) )
3029 def _max_random_instance_size():
3031 The maximum dimension of a random BilinearFormEJA.
3036 def random_instance(cls
, **kwargs
):
3038 Return a random instance of this algebra.
3040 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
3042 B
= matrix
.identity(ZZ
, n
)
3043 return cls(B
, **kwargs
)
3045 B11
= matrix
.identity(ZZ
, 1)
3046 M
= matrix
.random(ZZ
, n
-1)
3047 I
= matrix
.identity(ZZ
, n
-1)
3049 while alpha
.is_zero():
3050 alpha
= ZZ
.random_element().abs()
3051 B22
= M
.transpose()*M
+ alpha
*I
3053 from sage
.matrix
.special
import block_matrix
3054 B
= block_matrix(2,2, [ [B11
, ZZ(0) ],
3057 return cls(B
, **kwargs
)
3060 class JordanSpinEJA(BilinearFormEJA
):
3062 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
3063 with the usual inner product and jordan product ``x*y =
3064 (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
3069 sage: from mjo.eja.eja_algebra import JordanSpinEJA
3073 This multiplication table can be verified by hand::
3075 sage: J = JordanSpinEJA(4)
3076 sage: b0,b1,b2,b3 = J.gens()
3092 We can change the generator prefix::
3094 sage: JordanSpinEJA(2, prefix='B').gens()
3099 Ensure that we have the usual inner product on `R^n`::
3101 sage: set_random_seed()
3102 sage: J = JordanSpinEJA.random_instance()
3103 sage: x,y = J.random_elements(2)
3104 sage: actual = x.inner_product(y)
3105 sage: expected = x.to_vector().inner_product(y.to_vector())
3106 sage: actual == expected
3110 def __init__(self
, n
, *args
, **kwargs
):
3111 # This is a special case of the BilinearFormEJA with the
3112 # identity matrix as its bilinear form.
3113 B
= matrix
.identity(ZZ
, n
)
3115 # Don't orthonormalize because our basis is already
3116 # orthonormal with respect to our inner-product.
3117 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
3119 # But also don't pass check_field=False here, because the user
3120 # can pass in a field!
3121 super().__init
__(B
, *args
, **kwargs
)
3124 def _max_random_instance_size():
3126 The maximum dimension of a random JordanSpinEJA.
3131 def random_instance(cls
, **kwargs
):
3133 Return a random instance of this type of algebra.
3135 Needed here to override the implementation for ``BilinearFormEJA``.
3137 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
3138 return cls(n
, **kwargs
)
3141 class TrivialEJA(RationalBasisEJA
, ConcreteEJA
):
3143 The trivial Euclidean Jordan algebra consisting of only a zero element.
3147 sage: from mjo.eja.eja_algebra import TrivialEJA
3151 sage: J = TrivialEJA()
3158 sage: 7*J.one()*12*J.one()
3160 sage: J.one().inner_product(J.one())
3162 sage: J.one().norm()
3164 sage: J.one().subalgebra_generated_by()
3165 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
3170 def __init__(self
, **kwargs
):
3171 jordan_product
= lambda x
,y
: x
3172 inner_product
= lambda x
,y
: 0
3175 # New defaults for keyword arguments
3176 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
3177 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
3179 super().__init
__(basis
,
3185 # The rank is zero using my definition, namely the dimension of the
3186 # largest subalgebra generated by any element.
3187 self
.rank
.set_cache(0)
3188 self
.one
.set_cache( self
.zero() )
3191 def random_instance(cls
, **kwargs
):
3192 # We don't take a "size" argument so the superclass method is
3193 # inappropriate for us.
3194 return cls(**kwargs
)
3197 class CartesianProductEJA(FiniteDimensionalEJA
):
3199 The external (orthogonal) direct sum of two or more Euclidean
3200 Jordan algebras. Every Euclidean Jordan algebra decomposes into an
3201 orthogonal direct sum of simple Euclidean Jordan algebras which is
3202 then isometric to a Cartesian product, so no generality is lost by
3203 providing only this construction.
3207 sage: from mjo.eja.eja_algebra import (random_eja,
3208 ....: CartesianProductEJA,
3210 ....: JordanSpinEJA,
3211 ....: RealSymmetricEJA)
3215 The Jordan product is inherited from our factors and implemented by
3216 our CombinatorialFreeModule Cartesian product superclass::
3218 sage: set_random_seed()
3219 sage: J1 = HadamardEJA(2)
3220 sage: J2 = RealSymmetricEJA(2)
3221 sage: J = cartesian_product([J1,J2])
3222 sage: x,y = J.random_elements(2)
3226 The ability to retrieve the original factors is implemented by our
3227 CombinatorialFreeModule Cartesian product superclass::
3229 sage: J1 = HadamardEJA(2, field=QQ)
3230 sage: J2 = JordanSpinEJA(3, field=QQ)
3231 sage: J = cartesian_product([J1,J2])
3232 sage: J.cartesian_factors()
3233 (Euclidean Jordan algebra of dimension 2 over Rational Field,
3234 Euclidean Jordan algebra of dimension 3 over Rational Field)
3236 You can provide more than two factors::
3238 sage: J1 = HadamardEJA(2)
3239 sage: J2 = JordanSpinEJA(3)
3240 sage: J3 = RealSymmetricEJA(3)
3241 sage: cartesian_product([J1,J2,J3])
3242 Euclidean Jordan algebra of dimension 2 over Algebraic Real
3243 Field (+) Euclidean Jordan algebra of dimension 3 over Algebraic
3244 Real Field (+) Euclidean Jordan algebra of dimension 6 over
3245 Algebraic Real Field
3247 Rank is additive on a Cartesian product::
3249 sage: J1 = HadamardEJA(1)
3250 sage: J2 = RealSymmetricEJA(2)
3251 sage: J = cartesian_product([J1,J2])
3252 sage: J1.rank.clear_cache()
3253 sage: J2.rank.clear_cache()
3254 sage: J.rank.clear_cache()
3257 sage: J.rank() == J1.rank() + J2.rank()
3260 The same rank computation works over the rationals, with whatever
3263 sage: J1 = HadamardEJA(1, field=QQ, orthonormalize=False)
3264 sage: J2 = RealSymmetricEJA(2, field=QQ, orthonormalize=False)
3265 sage: J = cartesian_product([J1,J2])
3266 sage: J1.rank.clear_cache()
3267 sage: J2.rank.clear_cache()
3268 sage: J.rank.clear_cache()
3271 sage: J.rank() == J1.rank() + J2.rank()
3274 The product algebra will be associative if and only if all of its
3275 components are associative::
3277 sage: J1 = HadamardEJA(2)
3278 sage: J1.is_associative()
3280 sage: J2 = HadamardEJA(3)
3281 sage: J2.is_associative()
3283 sage: J3 = RealSymmetricEJA(3)
3284 sage: J3.is_associative()
3286 sage: CP1 = cartesian_product([J1,J2])
3287 sage: CP1.is_associative()
3289 sage: CP2 = cartesian_product([J1,J3])
3290 sage: CP2.is_associative()
3293 Cartesian products of Cartesian products work::
3295 sage: J1 = JordanSpinEJA(1)
3296 sage: J2 = JordanSpinEJA(1)
3297 sage: J3 = JordanSpinEJA(1)
3298 sage: J = cartesian_product([J1,cartesian_product([J2,J3])])
3299 sage: J.multiplication_table()
3300 +----++----+----+----+
3301 | * || b0 | b1 | b2 |
3302 +====++====+====+====+
3303 | b0 || b0 | 0 | 0 |
3304 +----++----+----+----+
3305 | b1 || 0 | b1 | 0 |
3306 +----++----+----+----+
3307 | b2 || 0 | 0 | b2 |
3308 +----++----+----+----+
3309 sage: HadamardEJA(3).multiplication_table()
3310 +----++----+----+----+
3311 | * || b0 | b1 | b2 |
3312 +====++====+====+====+
3313 | b0 || b0 | 0 | 0 |
3314 +----++----+----+----+
3315 | b1 || 0 | b1 | 0 |
3316 +----++----+----+----+
3317 | b2 || 0 | 0 | b2 |
3318 +----++----+----+----+
3322 All factors must share the same base field::
3324 sage: J1 = HadamardEJA(2, field=QQ)
3325 sage: J2 = RealSymmetricEJA(2)
3326 sage: CartesianProductEJA((J1,J2))
3327 Traceback (most recent call last):
3329 ValueError: all factors must share the same base field
3331 The cached unit element is the same one that would be computed::
3333 sage: set_random_seed() # long time
3334 sage: J1 = random_eja() # long time
3335 sage: J2 = random_eja() # long time
3336 sage: J = cartesian_product([J1,J2]) # long time
3337 sage: actual = J.one() # long time
3338 sage: J.one.clear_cache() # long time
3339 sage: expected = J.one() # long time
3340 sage: actual == expected # long time
3344 Element
= FiniteDimensionalEJAElement
3347 def __init__(self
, factors
, **kwargs
):
3352 self
._sets
= factors
3354 field
= factors
[0].base_ring()
3355 if not all( J
.base_ring() == field
for J
in factors
):
3356 raise ValueError("all factors must share the same base field")
3358 associative
= all( f
.is_associative() for f
in factors
)
3360 MS
= self
.matrix_space()
3364 for b
in factors
[i
].matrix_basis():
3369 basis
= tuple( MS(b
) for b
in basis
)
3371 # Define jordan/inner products that operate on that matrix_basis.
3372 def jordan_product(x
,y
):
3374 (factors
[i
](x
[i
])*factors
[i
](y
[i
])).to_matrix()
3378 def inner_product(x
, y
):
3380 factors
[i
](x
[i
]).inner_product(factors
[i
](y
[i
]))
3384 # There's no need to check the field since it already came
3385 # from an EJA. Likewise the axioms are guaranteed to be
3386 # satisfied, unless the guy writing this class sucks.
3388 # If you want the basis to be orthonormalized, orthonormalize
3390 FiniteDimensionalEJA
.__init
__(self
,
3395 orthonormalize
=False,
3396 associative
=associative
,
3397 cartesian_product
=True,
3401 ones
= tuple(J
.one().to_matrix() for J
in factors
)
3402 self
.one
.set_cache(self(ones
))
3403 self
.rank
.set_cache(sum(J
.rank() for J
in factors
))
3405 def cartesian_factors(self
):
3406 # Copy/pasted from CombinatorialFreeModule_CartesianProduct.
3409 def cartesian_factor(self
, i
):
3411 Return the ``i``th factor of this algebra.
3413 return self
._sets
[i
]
3416 # Copy/pasted from CombinatorialFreeModule_CartesianProduct.
3417 from sage
.categories
.cartesian_product
import cartesian_product
3418 return cartesian_product
.symbol
.join("%s" % factor
3419 for factor
in self
._sets
)
3421 def matrix_space(self
):
3423 Return the space that our matrix basis lives in as a Cartesian
3426 We don't simply use the ``cartesian_product()`` functor here
3427 because it acts differently on SageMath MatrixSpaces and our
3428 custom MatrixAlgebras, which are CombinatorialFreeModules. We
3429 always want the result to be represented (and indexed) as
3434 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
3436 ....: OctonionHermitianEJA,
3437 ....: RealSymmetricEJA)
3441 sage: J1 = HadamardEJA(1)
3442 sage: J2 = RealSymmetricEJA(2)
3443 sage: J = cartesian_product([J1,J2])
3444 sage: J.matrix_space()
3445 The Cartesian product of (Full MatrixSpace of 1 by 1 dense
3446 matrices over Algebraic Real Field, Full MatrixSpace of 2
3447 by 2 dense matrices over Algebraic Real Field)
3451 sage: J1 = ComplexHermitianEJA(1)
3452 sage: J2 = ComplexHermitianEJA(1)
3453 sage: J = cartesian_product([J1,J2])
3454 sage: J.one().to_matrix()[0]
3457 sage: J.one().to_matrix()[1]
3463 sage: J1 = OctonionHermitianEJA(1)
3464 sage: J2 = OctonionHermitianEJA(1)
3465 sage: J = cartesian_product([J1,J2])
3466 sage: J.one().to_matrix()[0]
3470 sage: J.one().to_matrix()[1]
3476 scalars
= self
.cartesian_factor(0).base_ring()
3478 # This category isn't perfect, but is good enough for what we
3480 cat
= MagmaticAlgebras(scalars
).FiniteDimensional().WithBasis()
3481 cat
= cat
.Unital().CartesianProducts()
3482 factors
= tuple( J
.matrix_space() for J
in self
.cartesian_factors() )
3484 from sage
.sets
.cartesian_product
import CartesianProduct
3485 return CartesianProduct(factors
, cat
)
3489 def cartesian_projection(self
, i
):
3493 sage: from mjo.eja.eja_algebra import (random_eja,
3494 ....: JordanSpinEJA,
3496 ....: RealSymmetricEJA,
3497 ....: ComplexHermitianEJA)
3501 The projection morphisms are Euclidean Jordan algebra
3504 sage: J1 = HadamardEJA(2)
3505 sage: J2 = RealSymmetricEJA(2)
3506 sage: J = cartesian_product([J1,J2])
3507 sage: J.cartesian_projection(0)
3508 Linear operator between finite-dimensional Euclidean Jordan
3509 algebras represented by the matrix:
3512 Domain: Euclidean Jordan algebra of dimension 2 over Algebraic
3513 Real Field (+) Euclidean Jordan algebra of dimension 3 over
3514 Algebraic Real Field
3515 Codomain: Euclidean Jordan algebra of dimension 2 over Algebraic
3517 sage: J.cartesian_projection(1)
3518 Linear operator between finite-dimensional Euclidean Jordan
3519 algebras represented by the matrix:
3523 Domain: Euclidean Jordan algebra of dimension 2 over Algebraic
3524 Real Field (+) Euclidean Jordan algebra of dimension 3 over
3525 Algebraic Real Field
3526 Codomain: Euclidean Jordan algebra of dimension 3 over Algebraic
3529 The projections work the way you'd expect on the vector
3530 representation of an element::
3532 sage: J1 = JordanSpinEJA(2)
3533 sage: J2 = ComplexHermitianEJA(2)
3534 sage: J = cartesian_product([J1,J2])
3535 sage: pi_left = J.cartesian_projection(0)
3536 sage: pi_right = J.cartesian_projection(1)
3537 sage: pi_left(J.one()).to_vector()
3539 sage: pi_right(J.one()).to_vector()
3541 sage: J.one().to_vector()
3546 The answer never changes::
3548 sage: set_random_seed()
3549 sage: J1 = random_eja()
3550 sage: J2 = random_eja()
3551 sage: J = cartesian_product([J1,J2])
3552 sage: P0 = J.cartesian_projection(0)
3553 sage: P1 = J.cartesian_projection(0)
3558 offset
= sum( self
.cartesian_factor(k
).dimension()
3560 Ji
= self
.cartesian_factor(i
)
3561 Pi
= self
._module
_morphism
(lambda j
: Ji
.monomial(j
- offset
),
3564 return FiniteDimensionalEJAOperator(self
,Ji
,Pi
.matrix())
3567 def cartesian_embedding(self
, i
):
3571 sage: from mjo.eja.eja_algebra import (random_eja,
3572 ....: JordanSpinEJA,
3574 ....: RealSymmetricEJA)
3578 The embedding morphisms are Euclidean Jordan algebra
3581 sage: J1 = HadamardEJA(2)
3582 sage: J2 = RealSymmetricEJA(2)
3583 sage: J = cartesian_product([J1,J2])
3584 sage: J.cartesian_embedding(0)
3585 Linear operator between finite-dimensional Euclidean Jordan
3586 algebras represented by the matrix:
3592 Domain: Euclidean Jordan algebra of dimension 2 over
3593 Algebraic Real Field
3594 Codomain: Euclidean Jordan algebra of dimension 2 over
3595 Algebraic Real Field (+) Euclidean Jordan algebra of
3596 dimension 3 over Algebraic Real Field
3597 sage: J.cartesian_embedding(1)
3598 Linear operator between finite-dimensional Euclidean Jordan
3599 algebras represented by the matrix:
3605 Domain: Euclidean Jordan algebra of dimension 3 over
3606 Algebraic Real Field
3607 Codomain: Euclidean Jordan algebra of dimension 2 over
3608 Algebraic Real Field (+) Euclidean Jordan algebra of
3609 dimension 3 over Algebraic Real Field
3611 The embeddings work the way you'd expect on the vector
3612 representation of an element::
3614 sage: J1 = JordanSpinEJA(3)
3615 sage: J2 = RealSymmetricEJA(2)
3616 sage: J = cartesian_product([J1,J2])
3617 sage: iota_left = J.cartesian_embedding(0)
3618 sage: iota_right = J.cartesian_embedding(1)
3619 sage: iota_left(J1.zero()) == J.zero()
3621 sage: iota_right(J2.zero()) == J.zero()
3623 sage: J1.one().to_vector()
3625 sage: iota_left(J1.one()).to_vector()
3627 sage: J2.one().to_vector()
3629 sage: iota_right(J2.one()).to_vector()
3631 sage: J.one().to_vector()
3636 The answer never changes::
3638 sage: set_random_seed()
3639 sage: J1 = random_eja()
3640 sage: J2 = random_eja()
3641 sage: J = cartesian_product([J1,J2])
3642 sage: E0 = J.cartesian_embedding(0)
3643 sage: E1 = J.cartesian_embedding(0)
3647 Composing a projection with the corresponding inclusion should
3648 produce the identity map, and mismatching them should produce
3651 sage: set_random_seed()
3652 sage: J1 = random_eja()
3653 sage: J2 = random_eja()
3654 sage: J = cartesian_product([J1,J2])
3655 sage: iota_left = J.cartesian_embedding(0)
3656 sage: iota_right = J.cartesian_embedding(1)
3657 sage: pi_left = J.cartesian_projection(0)
3658 sage: pi_right = J.cartesian_projection(1)
3659 sage: pi_left*iota_left == J1.one().operator()
3661 sage: pi_right*iota_right == J2.one().operator()
3663 sage: (pi_left*iota_right).is_zero()
3665 sage: (pi_right*iota_left).is_zero()
3669 offset
= sum( self
.cartesian_factor(k
).dimension()
3671 Ji
= self
.cartesian_factor(i
)
3672 Ei
= Ji
._module
_morphism
(lambda j
: self
.monomial(j
+ offset
),
3674 return FiniteDimensionalEJAOperator(Ji
,self
,Ei
.matrix())
3678 FiniteDimensionalEJA
.CartesianProduct
= CartesianProductEJA
3680 class RationalBasisCartesianProductEJA(CartesianProductEJA
,
3683 A separate class for products of algebras for which we know a
3688 sage: from mjo.eja.eja_algebra import (HadamardEJA,
3689 ....: JordanSpinEJA,
3690 ....: OctonionHermitianEJA,
3691 ....: RealSymmetricEJA)
3695 This gives us fast characteristic polynomial computations in
3696 product algebras, too::
3699 sage: J1 = JordanSpinEJA(2)
3700 sage: J2 = RealSymmetricEJA(3)
3701 sage: J = cartesian_product([J1,J2])
3702 sage: J.characteristic_polynomial_of().degree()
3709 The ``cartesian_product()`` function only uses the first factor to
3710 decide where the result will live; thus we have to be careful to
3711 check that all factors do indeed have a `_rational_algebra` member
3712 before we try to access it::
3714 sage: J1 = OctonionHermitianEJA(1) # no rational basis
3715 sage: J2 = HadamardEJA(2)
3716 sage: cartesian_product([J1,J2])
3717 Euclidean Jordan algebra of dimension 1 over Algebraic Real Field
3718 (+) Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
3719 sage: cartesian_product([J2,J1])
3720 Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
3721 (+) Euclidean Jordan algebra of dimension 1 over Algebraic Real Field
3724 def __init__(self
, algebras
, **kwargs
):
3725 CartesianProductEJA
.__init
__(self
, algebras
, **kwargs
)
3727 self
._rational
_algebra
= None
3728 if self
.vector_space().base_field() is not QQ
:
3729 if all( hasattr(r
, "_rational_algebra") for r
in algebras
):
3730 self
._rational
_algebra
= cartesian_product([
3731 r
._rational
_algebra
for r
in algebras
3735 RationalBasisEJA
.CartesianProduct
= RationalBasisCartesianProductEJA
3737 def random_eja(*args
, **kwargs
):
3738 J1
= ConcreteEJA
.random_instance(*args
, **kwargs
)
3740 # This might make Cartesian products appear roughly as often as
3741 # any other ConcreteEJA.
3742 if ZZ
.random_element(len(ConcreteEJA
.__subclasses
__()) + 1) == 0:
3743 # Use random_eja() again so we can get more than two factors.
3744 J2
= random_eja(*args
, **kwargs
)
3745 J
= cartesian_product([J1
,J2
])