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 if not all( all(b_i
in QQ
for b_i
in b
.list()) for b
in basis
):
1579 raise TypeError("basis not rational")
1581 super().__init
__(basis
,
1585 check_field
=check_field
,
1588 self
._rational
_algebra
= None
1590 # There's no point in constructing the extra algebra if this
1591 # one is already rational.
1593 # Note: the same Jordan and inner-products work here,
1594 # because they are necessarily defined with respect to
1595 # ambient coordinates and not any particular basis.
1596 self
._rational
_algebra
= FiniteDimensionalEJA(
1601 associative
=self
.is_associative(),
1602 orthonormalize
=False,
1607 def _charpoly_coefficients(self
):
1611 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
1612 ....: JordanSpinEJA)
1616 The base ring of the resulting polynomial coefficients is what
1617 it should be, and not the rationals (unless the algebra was
1618 already over the rationals)::
1620 sage: J = JordanSpinEJA(3)
1621 sage: J._charpoly_coefficients()
1622 (X1^2 - X2^2 - X3^2, -2*X1)
1623 sage: a0 = J._charpoly_coefficients()[0]
1625 Algebraic Real Field
1626 sage: a0.base_ring()
1627 Algebraic Real Field
1630 if self
._rational
_algebra
is None:
1631 # There's no need to construct *another* algebra over the
1632 # rationals if this one is already over the
1633 # rationals. Likewise, if we never orthonormalized our
1634 # basis, we might as well just use the given one.
1635 return super()._charpoly
_coefficients
()
1637 # Do the computation over the rationals. The answer will be
1638 # the same, because all we've done is a change of basis.
1639 # Then, change back from QQ to our real base ring
1640 a
= ( a_i
.change_ring(self
.base_ring())
1641 for a_i
in self
._rational
_algebra
._charpoly
_coefficients
() )
1643 if self
._deortho
_matrix
is None:
1644 # This can happen if our base ring was, say, AA and we
1645 # chose not to (or didn't need to) orthonormalize. It's
1646 # still faster to do the computations over QQ even if
1647 # the numbers in the boxes stay the same.
1650 # Otherwise, convert the coordinate variables back to the
1651 # deorthonormalized ones.
1652 R
= self
.coordinate_polynomial_ring()
1653 from sage
.modules
.free_module_element
import vector
1654 X
= vector(R
, R
.gens())
1655 BX
= self
._deortho
_matrix
*X
1657 subs_dict
= { X[i]: BX[i] for i in range(len(X)) }
1658 return tuple( a_i
.subs(subs_dict
) for a_i
in a
)
1660 class ConcreteEJA(FiniteDimensionalEJA
):
1662 A class for the Euclidean Jordan algebras that we know by name.
1664 These are the Jordan algebras whose basis, multiplication table,
1665 rank, and so on are known a priori. More to the point, they are
1666 the Euclidean Jordan algebras for which we are able to conjure up
1667 a "random instance."
1671 sage: from mjo.eja.eja_algebra import ConcreteEJA
1675 Our basis is normalized with respect to the algebra's inner
1676 product, unless we specify otherwise::
1678 sage: set_random_seed()
1679 sage: J = ConcreteEJA.random_instance()
1680 sage: all( b.norm() == 1 for b in J.gens() )
1683 Since our basis is orthonormal with respect to the algebra's inner
1684 product, and since we know that this algebra is an EJA, any
1685 left-multiplication operator's matrix will be symmetric because
1686 natural->EJA basis representation is an isometry and within the
1687 EJA the operator is self-adjoint by the Jordan axiom::
1689 sage: set_random_seed()
1690 sage: J = ConcreteEJA.random_instance()
1691 sage: x = J.random_element()
1692 sage: x.operator().is_self_adjoint()
1697 def _max_random_instance_size():
1699 Return an integer "size" that is an upper bound on the size of
1700 this algebra when it is used in a random test
1701 case. Unfortunately, the term "size" is ambiguous -- when
1702 dealing with `R^n` under either the Hadamard or Jordan spin
1703 product, the "size" refers to the dimension `n`. When dealing
1704 with a matrix algebra (real symmetric or complex/quaternion
1705 Hermitian), it refers to the size of the matrix, which is far
1706 less than the dimension of the underlying vector space.
1708 This method must be implemented in each subclass.
1710 raise NotImplementedError
1713 def random_instance(cls
, *args
, **kwargs
):
1715 Return a random instance of this type of algebra.
1717 This method should be implemented in each subclass.
1719 from sage
.misc
.prandom
import choice
1720 eja_class
= choice(cls
.__subclasses
__())
1722 # These all bubble up to the RationalBasisEJA superclass
1723 # constructor, so any (kw)args valid there are also valid
1725 return eja_class
.random_instance(*args
, **kwargs
)
1730 def jordan_product(X
,Y
):
1731 return (X
*Y
+ Y
*X
)/2
1734 def trace_inner_product(X
,Y
):
1736 A trace inner-product for matrices that aren't embedded in the
1737 reals. It takes MATRICES as arguments, not EJA elements.
1739 return (X
*Y
).trace().real()
1741 class RealEmbeddedMatrixEJA(MatrixEJA
):
1743 def dimension_over_reals():
1745 The dimension of this matrix's base ring over the reals.
1747 The reals are dimension one over themselves, obviously; that's
1748 just `\mathbb{R}^{1}`. Likewise, the complex numbers `a + bi`
1749 have dimension two. Finally, the quaternions have dimension
1750 four over the reals.
1752 This is used to determine the size of the matrix returned from
1753 :meth:`real_embed`, among other things.
1755 raise NotImplementedError
1758 def real_embed(cls
,M
):
1760 Embed the matrix ``M`` into a space of real matrices.
1762 The matrix ``M`` can have entries in any field at the moment:
1763 the real numbers, complex numbers, or quaternions. And although
1764 they are not a field, we can probably support octonions at some
1765 point, too. This function returns a real matrix that "acts like"
1766 the original with respect to matrix multiplication; i.e.
1768 real_embed(M*N) = real_embed(M)*real_embed(N)
1771 if M
.ncols() != M
.nrows():
1772 raise ValueError("the matrix 'M' must be square")
1777 def real_unembed(cls
,M
):
1779 The inverse of :meth:`real_embed`.
1781 if M
.ncols() != M
.nrows():
1782 raise ValueError("the matrix 'M' must be square")
1783 if not ZZ(M
.nrows()).mod(cls
.dimension_over_reals()).is_zero():
1784 raise ValueError("the matrix 'M' must be a real embedding")
1789 def trace_inner_product(cls
,X
,Y
):
1791 Compute the trace inner-product of two real-embeddings.
1795 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
1796 ....: QuaternionHermitianEJA)
1800 sage: set_random_seed()
1801 sage: J = ComplexHermitianEJA.random_instance()
1802 sage: x,y = J.random_elements(2)
1803 sage: Xe = x.to_matrix()
1804 sage: Ye = y.to_matrix()
1805 sage: X = J.real_unembed(Xe)
1806 sage: Y = J.real_unembed(Ye)
1807 sage: expected = (X*Y).trace().real()
1808 sage: actual = J.trace_inner_product(Xe,Ye)
1809 sage: actual == expected
1814 sage: set_random_seed()
1815 sage: J = QuaternionHermitianEJA.random_instance()
1816 sage: x,y = J.random_elements(2)
1817 sage: Xe = x.to_matrix()
1818 sage: Ye = y.to_matrix()
1819 sage: X = J.real_unembed(Xe)
1820 sage: Y = J.real_unembed(Ye)
1821 sage: expected = (X*Y).trace().coefficient_tuple()[0]
1822 sage: actual = J.trace_inner_product(Xe,Ye)
1823 sage: actual == expected
1827 # This does in fact compute the real part of the trace.
1828 # If we compute the trace of e.g. a complex matrix M,
1829 # then we do so by adding up its diagonal entries --
1830 # call them z_1 through z_n. The real embedding of z_1
1831 # will be a 2-by-2 REAL matrix [a, b; -b, a] whose trace
1832 # as a REAL matrix will be 2*a = 2*Re(z_1). And so forth.
1833 return (X
*Y
).trace()/cls
.dimension_over_reals()
1835 class RealSymmetricEJA(ConcreteEJA
, RationalBasisEJA
, MatrixEJA
):
1837 The rank-n simple EJA consisting of real symmetric n-by-n
1838 matrices, the usual symmetric Jordan product, and the trace inner
1839 product. It has dimension `(n^2 + n)/2` over the reals.
1843 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1847 sage: J = RealSymmetricEJA(2)
1848 sage: b0, b1, b2 = J.gens()
1856 In theory, our "field" can be any subfield of the reals::
1858 sage: RealSymmetricEJA(2, field=RDF, check_axioms=True)
1859 Euclidean Jordan algebra of dimension 3 over Real Double Field
1860 sage: RealSymmetricEJA(2, field=RR, check_axioms=True)
1861 Euclidean Jordan algebra of dimension 3 over Real Field with
1862 53 bits of precision
1866 The dimension of this algebra is `(n^2 + n) / 2`::
1868 sage: set_random_seed()
1869 sage: n_max = RealSymmetricEJA._max_random_instance_size()
1870 sage: n = ZZ.random_element(1, n_max)
1871 sage: J = RealSymmetricEJA(n)
1872 sage: J.dimension() == (n^2 + n)/2
1875 The Jordan multiplication is what we think it is::
1877 sage: set_random_seed()
1878 sage: J = RealSymmetricEJA.random_instance()
1879 sage: x,y = J.random_elements(2)
1880 sage: actual = (x*y).to_matrix()
1881 sage: X = x.to_matrix()
1882 sage: Y = y.to_matrix()
1883 sage: expected = (X*Y + Y*X)/2
1884 sage: actual == expected
1886 sage: J(expected) == x*y
1889 We can change the generator prefix::
1891 sage: RealSymmetricEJA(3, prefix='q').gens()
1892 (q0, q1, q2, q3, q4, q5)
1894 We can construct the (trivial) algebra of rank zero::
1896 sage: RealSymmetricEJA(0)
1897 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1901 def _denormalized_basis(cls
, n
, field
):
1903 Return a basis for the space of real symmetric n-by-n matrices.
1907 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1911 sage: set_random_seed()
1912 sage: n = ZZ.random_element(1,5)
1913 sage: B = RealSymmetricEJA._denormalized_basis(n,ZZ)
1914 sage: all( M.is_symmetric() for M in B)
1918 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1922 for j
in range(i
+1):
1923 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
1927 Sij
= Eij
+ Eij
.transpose()
1933 def _max_random_instance_size():
1934 return 4 # Dimension 10
1937 def random_instance(cls
, **kwargs
):
1939 Return a random instance of this type of algebra.
1941 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1942 return cls(n
, **kwargs
)
1944 def __init__(self
, n
, field
=AA
, **kwargs
):
1945 # We know this is a valid EJA, but will double-check
1946 # if the user passes check_axioms=True.
1947 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
1953 super().__init
__(self
._denormalized
_basis
(n
,field
),
1954 self
.jordan_product
,
1955 self
.trace_inner_product
,
1957 associative
=associative
,
1960 # TODO: this could be factored out somehow, but is left here
1961 # because the MatrixEJA is not presently a subclass of the
1962 # FDEJA class that defines rank() and one().
1963 self
.rank
.set_cache(n
)
1964 idV
= self
.matrix_space().one()
1965 self
.one
.set_cache(self(idV
))
1969 class ComplexMatrixEJA(RealEmbeddedMatrixEJA
):
1970 # A manual dictionary-cache for the complex_extension() method,
1971 # since apparently @classmethods can't also be @cached_methods.
1972 _complex_extension
= {}
1975 def complex_extension(cls
,field
):
1977 The complex field that we embed/unembed, as an extension
1978 of the given ``field``.
1980 if field
in cls
._complex
_extension
:
1981 return cls
._complex
_extension
[field
]
1983 # Sage doesn't know how to adjoin the complex "i" (the root of
1984 # x^2 + 1) to a field in a general way. Here, we just enumerate
1985 # all of the cases that I have cared to support so far.
1987 # Sage doesn't know how to embed AA into QQbar, i.e. how
1988 # to adjoin sqrt(-1) to AA.
1990 elif not field
.is_exact():
1992 F
= field
.complex_field()
1994 # Works for QQ and... maybe some other fields.
1995 R
= PolynomialRing(field
, 'z')
1997 F
= field
.extension(z
**2 + 1, 'I', embedding
=CLF(-1).sqrt())
1999 cls
._complex
_extension
[field
] = F
2003 def dimension_over_reals():
2007 def real_embed(cls
,M
):
2009 Embed the n-by-n complex matrix ``M`` into the space of real
2010 matrices of size 2n-by-2n via the map the sends each entry `z = a +
2011 bi` to the block matrix ``[[a,b],[-b,a]]``.
2015 sage: from mjo.eja.eja_algebra import ComplexMatrixEJA
2019 sage: F = QuadraticField(-1, 'I')
2020 sage: x1 = F(4 - 2*i)
2021 sage: x2 = F(1 + 2*i)
2024 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
2025 sage: ComplexMatrixEJA.real_embed(M)
2034 Embedding is a homomorphism (isomorphism, in fact)::
2036 sage: set_random_seed()
2037 sage: n = ZZ.random_element(3)
2038 sage: F = QuadraticField(-1, 'I')
2039 sage: X = random_matrix(F, n)
2040 sage: Y = random_matrix(F, n)
2041 sage: Xe = ComplexMatrixEJA.real_embed(X)
2042 sage: Ye = ComplexMatrixEJA.real_embed(Y)
2043 sage: XYe = ComplexMatrixEJA.real_embed(X*Y)
2048 super().real_embed(M
)
2051 # We don't need any adjoined elements...
2052 field
= M
.base_ring().base_ring()
2058 blocks
.append(matrix(field
, 2, [ [ a
, b
],
2061 return matrix
.block(field
, n
, blocks
)
2065 def real_unembed(cls
,M
):
2067 The inverse of _embed_complex_matrix().
2071 sage: from mjo.eja.eja_algebra import ComplexMatrixEJA
2075 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
2076 ....: [-2, 1, -4, 3],
2077 ....: [ 9, 10, 11, 12],
2078 ....: [-10, 9, -12, 11] ])
2079 sage: ComplexMatrixEJA.real_unembed(A)
2081 [ 10*I + 9 12*I + 11]
2085 Unembedding is the inverse of embedding::
2087 sage: set_random_seed()
2088 sage: F = QuadraticField(-1, 'I')
2089 sage: M = random_matrix(F, 3)
2090 sage: Me = ComplexMatrixEJA.real_embed(M)
2091 sage: ComplexMatrixEJA.real_unembed(Me) == M
2095 super().real_unembed(M
)
2097 d
= cls
.dimension_over_reals()
2098 F
= cls
.complex_extension(M
.base_ring())
2101 # Go top-left to bottom-right (reading order), converting every
2102 # 2-by-2 block we see to a single complex element.
2104 for k
in range(n
/d
):
2105 for j
in range(n
/d
):
2106 submat
= M
[d
*k
:d
*k
+d
,d
*j
:d
*j
+d
]
2107 if submat
[0,0] != submat
[1,1]:
2108 raise ValueError('bad on-diagonal submatrix')
2109 if submat
[0,1] != -submat
[1,0]:
2110 raise ValueError('bad off-diagonal submatrix')
2111 z
= submat
[0,0] + submat
[0,1]*i
2114 return matrix(F
, n
/d
, elements
)
2117 class ComplexHermitianEJA(ConcreteEJA
, RationalBasisEJA
, ComplexMatrixEJA
):
2119 The rank-n simple EJA consisting of complex Hermitian n-by-n
2120 matrices over the real numbers, the usual symmetric Jordan product,
2121 and the real-part-of-trace inner product. It has dimension `n^2` over
2126 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
2130 In theory, our "field" can be any subfield of the reals::
2132 sage: ComplexHermitianEJA(2, field=RDF, check_axioms=True)
2133 Euclidean Jordan algebra of dimension 4 over Real Double Field
2134 sage: ComplexHermitianEJA(2, field=RR, check_axioms=True)
2135 Euclidean Jordan algebra of dimension 4 over Real Field with
2136 53 bits of precision
2140 The dimension of this algebra is `n^2`::
2142 sage: set_random_seed()
2143 sage: n_max = ComplexHermitianEJA._max_random_instance_size()
2144 sage: n = ZZ.random_element(1, n_max)
2145 sage: J = ComplexHermitianEJA(n)
2146 sage: J.dimension() == n^2
2149 The Jordan multiplication is what we think it is::
2151 sage: set_random_seed()
2152 sage: J = ComplexHermitianEJA.random_instance()
2153 sage: x,y = J.random_elements(2)
2154 sage: actual = (x*y).to_matrix()
2155 sage: X = x.to_matrix()
2156 sage: Y = y.to_matrix()
2157 sage: expected = (X*Y + Y*X)/2
2158 sage: actual == expected
2160 sage: J(expected) == x*y
2163 We can change the generator prefix::
2165 sage: ComplexHermitianEJA(2, prefix='z').gens()
2168 We can construct the (trivial) algebra of rank zero::
2170 sage: ComplexHermitianEJA(0)
2171 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2176 def _denormalized_basis(cls
, n
, field
):
2178 Returns a basis for the space of complex Hermitian n-by-n matrices.
2180 Why do we embed these? Basically, because all of numerical linear
2181 algebra assumes that you're working with vectors consisting of `n`
2182 entries from a field and scalars from the same field. There's no way
2183 to tell SageMath that (for example) the vectors contain complex
2184 numbers, while the scalar field is real.
2188 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
2192 sage: set_random_seed()
2193 sage: n = ZZ.random_element(1,5)
2194 sage: B = ComplexHermitianEJA._denormalized_basis(n,ZZ)
2195 sage: all( M.is_symmetric() for M in B)
2199 R
= PolynomialRing(ZZ
, 'z')
2201 F
= ZZ
.extension(z
**2 + 1, 'I')
2204 # This is like the symmetric case, but we need to be careful:
2206 # * We want conjugate-symmetry, not just symmetry.
2207 # * The diagonal will (as a result) be real.
2210 Eij
= matrix
.zero(F
,n
)
2212 for j
in range(i
+1):
2216 Sij
= cls
.real_embed(Eij
)
2219 # The second one has a minus because it's conjugated.
2220 Eij
[j
,i
] = 1 # Eij = Eij + Eij.transpose()
2221 Sij_real
= cls
.real_embed(Eij
)
2223 # Eij = I*Eij - I*Eij.transpose()
2226 Sij_imag
= cls
.real_embed(Eij
)
2232 # Since we embedded the entries, we can drop back to the
2233 # desired real "field" instead of the extension "F".
2234 return tuple( s
.change_ring(field
) for s
in S
)
2237 def __init__(self
, n
, field
=AA
, **kwargs
):
2238 # We know this is a valid EJA, but will double-check
2239 # if the user passes check_axioms=True.
2240 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2246 super().__init
__(self
._denormalized
_basis
(n
,field
),
2247 self
.jordan_product
,
2248 self
.trace_inner_product
,
2250 associative
=associative
,
2252 # TODO: this could be factored out somehow, but is left here
2253 # because the MatrixEJA is not presently a subclass of the
2254 # FDEJA class that defines rank() and one().
2255 self
.rank
.set_cache(n
)
2256 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
2257 self
.one
.set_cache(self(idV
))
2260 def _max_random_instance_size():
2261 return 3 # Dimension 9
2264 def random_instance(cls
, **kwargs
):
2266 Return a random instance of this type of algebra.
2268 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2269 return cls(n
, **kwargs
)
2271 class QuaternionMatrixEJA(RealEmbeddedMatrixEJA
):
2273 # A manual dictionary-cache for the quaternion_extension() method,
2274 # since apparently @classmethods can't also be @cached_methods.
2275 _quaternion_extension
= {}
2278 def quaternion_extension(cls
,field
):
2280 The quaternion field that we embed/unembed, as an extension
2281 of the given ``field``.
2283 if field
in cls
._quaternion
_extension
:
2284 return cls
._quaternion
_extension
[field
]
2286 Q
= QuaternionAlgebra(field
,-1,-1)
2288 cls
._quaternion
_extension
[field
] = Q
2292 def dimension_over_reals():
2296 def real_embed(cls
,M
):
2298 Embed the n-by-n quaternion matrix ``M`` into the space of real
2299 matrices of size 4n-by-4n by first sending each quaternion entry `z
2300 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
2301 c+di],[-c + di, a-bi]]`, and then embedding those into a real
2306 sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA
2310 sage: Q = QuaternionAlgebra(QQ,-1,-1)
2311 sage: i,j,k = Q.gens()
2312 sage: x = 1 + 2*i + 3*j + 4*k
2313 sage: M = matrix(Q, 1, [[x]])
2314 sage: QuaternionMatrixEJA.real_embed(M)
2320 Embedding is a homomorphism (isomorphism, in fact)::
2322 sage: set_random_seed()
2323 sage: n = ZZ.random_element(2)
2324 sage: Q = QuaternionAlgebra(QQ,-1,-1)
2325 sage: X = random_matrix(Q, n)
2326 sage: Y = random_matrix(Q, n)
2327 sage: Xe = QuaternionMatrixEJA.real_embed(X)
2328 sage: Ye = QuaternionMatrixEJA.real_embed(Y)
2329 sage: XYe = QuaternionMatrixEJA.real_embed(X*Y)
2334 super().real_embed(M
)
2335 quaternions
= M
.base_ring()
2338 F
= QuadraticField(-1, 'I')
2343 t
= z
.coefficient_tuple()
2348 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
2349 [-c
+ d
*i
, a
- b
*i
]])
2350 realM
= ComplexMatrixEJA
.real_embed(cplxM
)
2351 blocks
.append(realM
)
2353 # We should have real entries by now, so use the realest field
2354 # we've got for the return value.
2355 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
2360 def real_unembed(cls
,M
):
2362 The inverse of _embed_quaternion_matrix().
2366 sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA
2370 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
2371 ....: [-2, 1, -4, 3],
2372 ....: [-3, 4, 1, -2],
2373 ....: [-4, -3, 2, 1]])
2374 sage: QuaternionMatrixEJA.real_unembed(M)
2375 [1 + 2*i + 3*j + 4*k]
2379 Unembedding is the inverse of embedding::
2381 sage: set_random_seed()
2382 sage: Q = QuaternionAlgebra(QQ, -1, -1)
2383 sage: M = random_matrix(Q, 3)
2384 sage: Me = QuaternionMatrixEJA.real_embed(M)
2385 sage: QuaternionMatrixEJA.real_unembed(Me) == M
2389 super().real_unembed(M
)
2391 d
= cls
.dimension_over_reals()
2393 # Use the base ring of the matrix to ensure that its entries can be
2394 # multiplied by elements of the quaternion algebra.
2395 Q
= cls
.quaternion_extension(M
.base_ring())
2398 # Go top-left to bottom-right (reading order), converting every
2399 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
2402 for l
in range(n
/d
):
2403 for m
in range(n
/d
):
2404 submat
= ComplexMatrixEJA
.real_unembed(
2405 M
[d
*l
:d
*l
+d
,d
*m
:d
*m
+d
] )
2406 if submat
[0,0] != submat
[1,1].conjugate():
2407 raise ValueError('bad on-diagonal submatrix')
2408 if submat
[0,1] != -submat
[1,0].conjugate():
2409 raise ValueError('bad off-diagonal submatrix')
2410 z
= submat
[0,0].real()
2411 z
+= submat
[0,0].imag()*i
2412 z
+= submat
[0,1].real()*j
2413 z
+= submat
[0,1].imag()*k
2416 return matrix(Q
, n
/d
, elements
)
2419 class QuaternionHermitianEJA(ConcreteEJA
,
2421 QuaternionMatrixEJA
):
2423 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
2424 matrices, the usual symmetric Jordan product, and the
2425 real-part-of-trace inner product. It has dimension `2n^2 - n` over
2430 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2434 In theory, our "field" can be any subfield of the reals::
2436 sage: QuaternionHermitianEJA(2, field=RDF, check_axioms=True)
2437 Euclidean Jordan algebra of dimension 6 over Real Double Field
2438 sage: QuaternionHermitianEJA(2, field=RR, check_axioms=True)
2439 Euclidean Jordan algebra of dimension 6 over Real Field with
2440 53 bits of precision
2444 The dimension of this algebra is `2*n^2 - n`::
2446 sage: set_random_seed()
2447 sage: n_max = QuaternionHermitianEJA._max_random_instance_size()
2448 sage: n = ZZ.random_element(1, n_max)
2449 sage: J = QuaternionHermitianEJA(n)
2450 sage: J.dimension() == 2*(n^2) - n
2453 The Jordan multiplication is what we think it is::
2455 sage: set_random_seed()
2456 sage: J = QuaternionHermitianEJA.random_instance()
2457 sage: x,y = J.random_elements(2)
2458 sage: actual = (x*y).to_matrix()
2459 sage: X = x.to_matrix()
2460 sage: Y = y.to_matrix()
2461 sage: expected = (X*Y + Y*X)/2
2462 sage: actual == expected
2464 sage: J(expected) == x*y
2467 We can change the generator prefix::
2469 sage: QuaternionHermitianEJA(2, prefix='a').gens()
2470 (a0, a1, a2, a3, a4, a5)
2472 We can construct the (trivial) algebra of rank zero::
2474 sage: QuaternionHermitianEJA(0)
2475 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2479 def _denormalized_basis(cls
, n
, field
):
2481 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
2483 Why do we embed these? Basically, because all of numerical
2484 linear algebra assumes that you're working with vectors consisting
2485 of `n` entries from a field and scalars from the same field. There's
2486 no way to tell SageMath that (for example) the vectors contain
2487 complex numbers, while the scalar field is real.
2491 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2495 sage: set_random_seed()
2496 sage: n = ZZ.random_element(1,5)
2497 sage: B = QuaternionHermitianEJA._denormalized_basis(n,ZZ)
2498 sage: all( M.is_symmetric() for M in B )
2502 Q
= QuaternionAlgebra(QQ
,-1,-1)
2505 # This is like the symmetric case, but we need to be careful:
2507 # * We want conjugate-symmetry, not just symmetry.
2508 # * The diagonal will (as a result) be real.
2511 Eij
= matrix
.zero(Q
,n
)
2513 for j
in range(i
+1):
2517 Sij
= cls
.real_embed(Eij
)
2520 # The second, third, and fourth ones have a minus
2521 # because they're conjugated.
2522 # Eij = Eij + Eij.transpose()
2524 Sij_real
= cls
.real_embed(Eij
)
2526 # Eij = I*(Eij - Eij.transpose())
2529 Sij_I
= cls
.real_embed(Eij
)
2531 # Eij = J*(Eij - Eij.transpose())
2534 Sij_J
= cls
.real_embed(Eij
)
2536 # Eij = K*(Eij - Eij.transpose())
2539 Sij_K
= cls
.real_embed(Eij
)
2545 # Since we embedded the entries, we can drop back to the
2546 # desired real "field" instead of the quaternion algebra "Q".
2547 return tuple( s
.change_ring(field
) for s
in S
)
2550 def __init__(self
, n
, field
=AA
, **kwargs
):
2551 # We know this is a valid EJA, but will double-check
2552 # if the user passes check_axioms=True.
2553 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2559 super().__init
__(self
._denormalized
_basis
(n
,field
),
2560 self
.jordan_product
,
2561 self
.trace_inner_product
,
2563 associative
=associative
,
2566 # TODO: this could be factored out somehow, but is left here
2567 # because the MatrixEJA is not presently a subclass of the
2568 # FDEJA class that defines rank() and one().
2569 self
.rank
.set_cache(n
)
2570 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
2571 self
.one
.set_cache(self(idV
))
2575 def _max_random_instance_size():
2577 The maximum rank of a random QuaternionHermitianEJA.
2579 return 2 # Dimension 6
2582 def random_instance(cls
, **kwargs
):
2584 Return a random instance of this type of algebra.
2586 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2587 return cls(n
, **kwargs
)
2589 class OctonionHermitianEJA(ConcreteEJA
, MatrixEJA
):
2593 sage: from mjo.eja.eja_algebra import (FiniteDimensionalEJA,
2594 ....: OctonionHermitianEJA)
2598 The 3-by-3 algebra satisfies the axioms of an EJA::
2600 sage: OctonionHermitianEJA(3, # long time
2601 ....: field=QQ, # long time
2602 ....: orthonormalize=False, # long time
2603 ....: check_axioms=True) # long time
2604 Euclidean Jordan algebra of dimension 27 over Rational Field
2606 After a change-of-basis, the 2-by-2 algebra has the same
2607 multiplication table as the ten-dimensional Jordan spin algebra::
2609 sage: b = OctonionHermitianEJA._denormalized_basis(2,QQ)
2610 sage: basis = (b[0] + b[9],) + b[1:9] + (b[0] - b[9],)
2611 sage: jp = OctonionHermitianEJA.jordan_product
2612 sage: ip = OctonionHermitianEJA.trace_inner_product
2613 sage: J = FiniteDimensionalEJA(basis,
2617 ....: orthonormalize=False)
2618 sage: J.multiplication_table()
2619 +----++----+----+----+----+----+----+----+----+----+----+
2620 | * || b0 | b1 | b2 | b3 | b4 | b5 | b6 | b7 | b8 | b9 |
2621 +====++====+====+====+====+====+====+====+====+====+====+
2622 | b0 || b0 | b1 | b2 | b3 | b4 | b5 | b6 | b7 | b8 | b9 |
2623 +----++----+----+----+----+----+----+----+----+----+----+
2624 | b1 || b1 | b0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
2625 +----++----+----+----+----+----+----+----+----+----+----+
2626 | b2 || b2 | 0 | b0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
2627 +----++----+----+----+----+----+----+----+----+----+----+
2628 | b3 || b3 | 0 | 0 | b0 | 0 | 0 | 0 | 0 | 0 | 0 |
2629 +----++----+----+----+----+----+----+----+----+----+----+
2630 | b4 || b4 | 0 | 0 | 0 | b0 | 0 | 0 | 0 | 0 | 0 |
2631 +----++----+----+----+----+----+----+----+----+----+----+
2632 | b5 || b5 | 0 | 0 | 0 | 0 | b0 | 0 | 0 | 0 | 0 |
2633 +----++----+----+----+----+----+----+----+----+----+----+
2634 | b6 || b6 | 0 | 0 | 0 | 0 | 0 | b0 | 0 | 0 | 0 |
2635 +----++----+----+----+----+----+----+----+----+----+----+
2636 | b7 || b7 | 0 | 0 | 0 | 0 | 0 | 0 | b0 | 0 | 0 |
2637 +----++----+----+----+----+----+----+----+----+----+----+
2638 | b8 || b8 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | b0 | 0 |
2639 +----++----+----+----+----+----+----+----+----+----+----+
2640 | b9 || b9 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | b0 |
2641 +----++----+----+----+----+----+----+----+----+----+----+
2645 We can actually construct the 27-dimensional Albert algebra,
2646 and we get the right unit element if we recompute it::
2648 sage: J = OctonionHermitianEJA(3, # long time
2649 ....: field=QQ, # long time
2650 ....: orthonormalize=False) # long time
2651 sage: J.one.clear_cache() # long time
2652 sage: J.one() # long time
2654 sage: J.one().to_matrix() # long time
2663 The 2-by-2 algebra is isomorphic to the ten-dimensional Jordan
2664 spin algebra, but just to be sure, we recompute its rank::
2666 sage: J = OctonionHermitianEJA(2, # long time
2667 ....: field=QQ, # long time
2668 ....: orthonormalize=False) # long time
2669 sage: J.rank.clear_cache() # long time
2670 sage: J.rank() # long time
2675 def _max_random_instance_size():
2677 The maximum rank of a random QuaternionHermitianEJA.
2679 return 1 # Dimension 1
2682 def random_instance(cls
, **kwargs
):
2684 Return a random instance of this type of algebra.
2686 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2687 return cls(n
, **kwargs
)
2689 def __init__(self
, n
, field
=AA
, **kwargs
):
2691 # Otherwise we don't get an EJA.
2692 raise ValueError("n cannot exceed 3")
2694 # We know this is a valid EJA, but will double-check
2695 # if the user passes check_axioms=True.
2696 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2698 super().__init
__(self
._denormalized
_basis
(n
,field
),
2699 self
.jordan_product
,
2700 self
.trace_inner_product
,
2704 # TODO: this could be factored out somehow, but is left here
2705 # because the MatrixEJA is not presently a subclass of the
2706 # FDEJA class that defines rank() and one().
2707 self
.rank
.set_cache(n
)
2708 idV
= self
.matrix_space().one()
2709 self
.one
.set_cache(self(idV
))
2713 def _denormalized_basis(cls
, n
, field
):
2715 Returns a basis for the space of octonion Hermitian n-by-n
2720 sage: from mjo.eja.eja_algebra import OctonionHermitianEJA
2724 sage: B = OctonionHermitianEJA._denormalized_basis(3,QQ)
2725 sage: all( M.is_hermitian() for M in B )
2731 from mjo
.octonions
import OctonionMatrixAlgebra
2732 MS
= OctonionMatrixAlgebra(n
, scalars
=field
)
2733 es
= MS
.entry_algebra().gens()
2737 for j
in range(i
+1):
2739 E_ii
= MS
.monomial( (i
,j
,es
[0]) )
2743 E_ij
= MS
.monomial( (i
,j
,e
) )
2745 # If the conjugate has a negative sign in front
2746 # of it, (j,i,ec) won't be a monomial!
2747 if (j
,i
,ec
) in MS
.indices():
2748 E_ij
+= MS
.monomial( (j
,i
,ec
) )
2750 E_ij
-= MS
.monomial( (j
,i
,-ec
) )
2753 return tuple( basis
)
2756 def trace_inner_product(X
,Y
):
2758 The octonions don't know that the reals are embedded in them,
2759 so we have to take the e0 component ourselves.
2763 sage: from mjo.eja.eja_algebra import OctonionHermitianEJA
2767 sage: J = OctonionHermitianEJA(2,field=QQ,orthonormalize=False)
2768 sage: I = J.one().to_matrix()
2769 sage: J.trace_inner_product(I, -I)
2773 return (X
*Y
).trace().real().coefficient(0)
2776 class AlbertEJA(OctonionHermitianEJA
):
2778 The Albert algebra is the algebra of three-by-three Hermitian
2779 matrices whose entries are octonions.
2783 from mjo.eja.eja_algebra import AlbertEJA
2787 sage: AlbertEJA(field=QQ, orthonormalize=False)
2788 Euclidean Jordan algebra of dimension 27 over Rational Field
2789 sage: AlbertEJA() # long time
2790 Euclidean Jordan algebra of dimension 27 over Algebraic Real Field
2793 def __init__(self
, *args
, **kwargs
):
2794 super().__init
__(3, *args
, **kwargs
)
2797 class HadamardEJA(ConcreteEJA
, RationalBasisEJA
):
2799 Return the Euclidean Jordan Algebra corresponding to the set
2800 `R^n` under the Hadamard product.
2802 Note: this is nothing more than the Cartesian product of ``n``
2803 copies of the spin algebra. Once Cartesian product algebras
2804 are implemented, this can go.
2808 sage: from mjo.eja.eja_algebra import HadamardEJA
2812 This multiplication table can be verified by hand::
2814 sage: J = HadamardEJA(3)
2815 sage: b0,b1,b2 = J.gens()
2831 We can change the generator prefix::
2833 sage: HadamardEJA(3, prefix='r').gens()
2837 def __init__(self
, n
, field
=AA
, **kwargs
):
2839 jordan_product
= lambda x
,y
: x
2840 inner_product
= lambda x
,y
: x
2842 def jordan_product(x
,y
):
2844 return P( xi
*yi
for (xi
,yi
) in zip(x
,y
) )
2846 def inner_product(x
,y
):
2849 # New defaults for keyword arguments. Don't orthonormalize
2850 # because our basis is already orthonormal with respect to our
2851 # inner-product. Don't check the axioms, because we know this
2852 # is a valid EJA... but do double-check if the user passes
2853 # check_axioms=True. Note: we DON'T override the "check_field"
2854 # default here, because the user can pass in a field!
2855 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2856 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2858 column_basis
= tuple( b
.column()
2859 for b
in FreeModule(field
, n
).basis() )
2860 super().__init
__(column_basis
,
2866 self
.rank
.set_cache(n
)
2869 self
.one
.set_cache( self
.zero() )
2871 self
.one
.set_cache( sum(self
.gens()) )
2874 def _max_random_instance_size():
2876 The maximum dimension of a random HadamardEJA.
2881 def random_instance(cls
, **kwargs
):
2883 Return a random instance of this type of algebra.
2885 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2886 return cls(n
, **kwargs
)
2889 class BilinearFormEJA(ConcreteEJA
, RationalBasisEJA
):
2891 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2892 with the half-trace inner product and jordan product ``x*y =
2893 (<Bx,y>,y_bar>, x0*y_bar + y0*x_bar)`` where `B = 1 \times B22` is
2894 a symmetric positive-definite "bilinear form" matrix. Its
2895 dimension is the size of `B`, and it has rank two in dimensions
2896 larger than two. It reduces to the ``JordanSpinEJA`` when `B` is
2897 the identity matrix of order ``n``.
2899 We insist that the one-by-one upper-left identity block of `B` be
2900 passed in as well so that we can be passed a matrix of size zero
2901 to construct a trivial algebra.
2905 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
2906 ....: JordanSpinEJA)
2910 When no bilinear form is specified, the identity matrix is used,
2911 and the resulting algebra is the Jordan spin algebra::
2913 sage: B = matrix.identity(AA,3)
2914 sage: J0 = BilinearFormEJA(B)
2915 sage: J1 = JordanSpinEJA(3)
2916 sage: J0.multiplication_table() == J0.multiplication_table()
2919 An error is raised if the matrix `B` does not correspond to a
2920 positive-definite bilinear form::
2922 sage: B = matrix.random(QQ,2,3)
2923 sage: J = BilinearFormEJA(B)
2924 Traceback (most recent call last):
2926 ValueError: bilinear form is not positive-definite
2927 sage: B = matrix.zero(QQ,3)
2928 sage: J = BilinearFormEJA(B)
2929 Traceback (most recent call last):
2931 ValueError: bilinear form is not positive-definite
2935 We can create a zero-dimensional algebra::
2937 sage: B = matrix.identity(AA,0)
2938 sage: J = BilinearFormEJA(B)
2942 We can check the multiplication condition given in the Jordan, von
2943 Neumann, and Wigner paper (and also discussed on my "On the
2944 symmetry..." paper). Note that this relies heavily on the standard
2945 choice of basis, as does anything utilizing the bilinear form
2946 matrix. We opt not to orthonormalize the basis, because if we
2947 did, we would have to normalize the `s_{i}` in a similar manner::
2949 sage: set_random_seed()
2950 sage: n = ZZ.random_element(5)
2951 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2952 sage: B11 = matrix.identity(QQ,1)
2953 sage: B22 = M.transpose()*M
2954 sage: B = block_matrix(2,2,[ [B11,0 ],
2956 sage: J = BilinearFormEJA(B, orthonormalize=False)
2957 sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
2958 sage: V = J.vector_space()
2959 sage: sis = [ J( V([0] + (M.inverse()*ei).list()).column() )
2960 ....: for ei in eis ]
2961 sage: actual = [ sis[i]*sis[j]
2962 ....: for i in range(n-1)
2963 ....: for j in range(n-1) ]
2964 sage: expected = [ J.one() if i == j else J.zero()
2965 ....: for i in range(n-1)
2966 ....: for j in range(n-1) ]
2967 sage: actual == expected
2971 def __init__(self
, B
, field
=AA
, **kwargs
):
2972 # The matrix "B" is supplied by the user in most cases,
2973 # so it makes sense to check whether or not its positive-
2974 # definite unless we are specifically asked not to...
2975 if ("check_axioms" not in kwargs
) or kwargs
["check_axioms"]:
2976 if not B
.is_positive_definite():
2977 raise ValueError("bilinear form is not positive-definite")
2979 # However, all of the other data for this EJA is computed
2980 # by us in manner that guarantees the axioms are
2981 # satisfied. So, again, unless we are specifically asked to
2982 # verify things, we'll skip the rest of the checks.
2983 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2985 def inner_product(x
,y
):
2986 return (y
.T
*B
*x
)[0,0]
2988 def jordan_product(x
,y
):
2994 z0
= inner_product(y
,x
)
2995 zbar
= y0
*xbar
+ x0
*ybar
2996 return P([z0
] + zbar
.list())
2999 column_basis
= tuple( b
.column()
3000 for b
in FreeModule(field
, n
).basis() )
3002 # TODO: I haven't actually checked this, but it seems legit.
3007 super().__init
__(column_basis
,
3011 associative
=associative
,
3014 # The rank of this algebra is two, unless we're in a
3015 # one-dimensional ambient space (because the rank is bounded
3016 # by the ambient dimension).
3017 self
.rank
.set_cache(min(n
,2))
3020 self
.one
.set_cache( self
.zero() )
3022 self
.one
.set_cache( self
.monomial(0) )
3025 def _max_random_instance_size():
3027 The maximum dimension of a random BilinearFormEJA.
3032 def random_instance(cls
, **kwargs
):
3034 Return a random instance of this algebra.
3036 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
3038 B
= matrix
.identity(ZZ
, n
)
3039 return cls(B
, **kwargs
)
3041 B11
= matrix
.identity(ZZ
, 1)
3042 M
= matrix
.random(ZZ
, n
-1)
3043 I
= matrix
.identity(ZZ
, n
-1)
3045 while alpha
.is_zero():
3046 alpha
= ZZ
.random_element().abs()
3047 B22
= M
.transpose()*M
+ alpha
*I
3049 from sage
.matrix
.special
import block_matrix
3050 B
= block_matrix(2,2, [ [B11
, ZZ(0) ],
3053 return cls(B
, **kwargs
)
3056 class JordanSpinEJA(BilinearFormEJA
):
3058 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
3059 with the usual inner product and jordan product ``x*y =
3060 (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
3065 sage: from mjo.eja.eja_algebra import JordanSpinEJA
3069 This multiplication table can be verified by hand::
3071 sage: J = JordanSpinEJA(4)
3072 sage: b0,b1,b2,b3 = J.gens()
3088 We can change the generator prefix::
3090 sage: JordanSpinEJA(2, prefix='B').gens()
3095 Ensure that we have the usual inner product on `R^n`::
3097 sage: set_random_seed()
3098 sage: J = JordanSpinEJA.random_instance()
3099 sage: x,y = J.random_elements(2)
3100 sage: actual = x.inner_product(y)
3101 sage: expected = x.to_vector().inner_product(y.to_vector())
3102 sage: actual == expected
3106 def __init__(self
, n
, *args
, **kwargs
):
3107 # This is a special case of the BilinearFormEJA with the
3108 # identity matrix as its bilinear form.
3109 B
= matrix
.identity(ZZ
, n
)
3111 # Don't orthonormalize because our basis is already
3112 # orthonormal with respect to our inner-product.
3113 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
3115 # But also don't pass check_field=False here, because the user
3116 # can pass in a field!
3117 super().__init
__(B
, *args
, **kwargs
)
3120 def _max_random_instance_size():
3122 The maximum dimension of a random JordanSpinEJA.
3127 def random_instance(cls
, **kwargs
):
3129 Return a random instance of this type of algebra.
3131 Needed here to override the implementation for ``BilinearFormEJA``.
3133 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
3134 return cls(n
, **kwargs
)
3137 class TrivialEJA(ConcreteEJA
, RationalBasisEJA
):
3139 The trivial Euclidean Jordan algebra consisting of only a zero element.
3143 sage: from mjo.eja.eja_algebra import TrivialEJA
3147 sage: J = TrivialEJA()
3154 sage: 7*J.one()*12*J.one()
3156 sage: J.one().inner_product(J.one())
3158 sage: J.one().norm()
3160 sage: J.one().subalgebra_generated_by()
3161 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
3166 def __init__(self
, **kwargs
):
3167 jordan_product
= lambda x
,y
: x
3168 inner_product
= lambda x
,y
: 0
3171 # New defaults for keyword arguments
3172 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
3173 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
3175 super().__init
__(basis
,
3181 # The rank is zero using my definition, namely the dimension of the
3182 # largest subalgebra generated by any element.
3183 self
.rank
.set_cache(0)
3184 self
.one
.set_cache( self
.zero() )
3187 def random_instance(cls
, **kwargs
):
3188 # We don't take a "size" argument so the superclass method is
3189 # inappropriate for us.
3190 return cls(**kwargs
)
3193 class CartesianProductEJA(FiniteDimensionalEJA
):
3195 The external (orthogonal) direct sum of two or more Euclidean
3196 Jordan algebras. Every Euclidean Jordan algebra decomposes into an
3197 orthogonal direct sum of simple Euclidean Jordan algebras which is
3198 then isometric to a Cartesian product, so no generality is lost by
3199 providing only this construction.
3203 sage: from mjo.eja.eja_algebra import (random_eja,
3204 ....: CartesianProductEJA,
3206 ....: JordanSpinEJA,
3207 ....: RealSymmetricEJA)
3211 The Jordan product is inherited from our factors and implemented by
3212 our CombinatorialFreeModule Cartesian product superclass::
3214 sage: set_random_seed()
3215 sage: J1 = HadamardEJA(2)
3216 sage: J2 = RealSymmetricEJA(2)
3217 sage: J = cartesian_product([J1,J2])
3218 sage: x,y = J.random_elements(2)
3222 The ability to retrieve the original factors is implemented by our
3223 CombinatorialFreeModule Cartesian product superclass::
3225 sage: J1 = HadamardEJA(2, field=QQ)
3226 sage: J2 = JordanSpinEJA(3, field=QQ)
3227 sage: J = cartesian_product([J1,J2])
3228 sage: J.cartesian_factors()
3229 (Euclidean Jordan algebra of dimension 2 over Rational Field,
3230 Euclidean Jordan algebra of dimension 3 over Rational Field)
3232 You can provide more than two factors::
3234 sage: J1 = HadamardEJA(2)
3235 sage: J2 = JordanSpinEJA(3)
3236 sage: J3 = RealSymmetricEJA(3)
3237 sage: cartesian_product([J1,J2,J3])
3238 Euclidean Jordan algebra of dimension 2 over Algebraic Real
3239 Field (+) Euclidean Jordan algebra of dimension 3 over Algebraic
3240 Real Field (+) Euclidean Jordan algebra of dimension 6 over
3241 Algebraic Real Field
3243 Rank is additive on a Cartesian product::
3245 sage: J1 = HadamardEJA(1)
3246 sage: J2 = RealSymmetricEJA(2)
3247 sage: J = cartesian_product([J1,J2])
3248 sage: J1.rank.clear_cache()
3249 sage: J2.rank.clear_cache()
3250 sage: J.rank.clear_cache()
3253 sage: J.rank() == J1.rank() + J2.rank()
3256 The same rank computation works over the rationals, with whatever
3259 sage: J1 = HadamardEJA(1, field=QQ, orthonormalize=False)
3260 sage: J2 = RealSymmetricEJA(2, field=QQ, orthonormalize=False)
3261 sage: J = cartesian_product([J1,J2])
3262 sage: J1.rank.clear_cache()
3263 sage: J2.rank.clear_cache()
3264 sage: J.rank.clear_cache()
3267 sage: J.rank() == J1.rank() + J2.rank()
3270 The product algebra will be associative if and only if all of its
3271 components are associative::
3273 sage: J1 = HadamardEJA(2)
3274 sage: J1.is_associative()
3276 sage: J2 = HadamardEJA(3)
3277 sage: J2.is_associative()
3279 sage: J3 = RealSymmetricEJA(3)
3280 sage: J3.is_associative()
3282 sage: CP1 = cartesian_product([J1,J2])
3283 sage: CP1.is_associative()
3285 sage: CP2 = cartesian_product([J1,J3])
3286 sage: CP2.is_associative()
3289 Cartesian products of Cartesian products work::
3291 sage: J1 = JordanSpinEJA(1)
3292 sage: J2 = JordanSpinEJA(1)
3293 sage: J3 = JordanSpinEJA(1)
3294 sage: J = cartesian_product([J1,cartesian_product([J2,J3])])
3295 sage: J.multiplication_table()
3296 +----++----+----+----+
3297 | * || b0 | b1 | b2 |
3298 +====++====+====+====+
3299 | b0 || b0 | 0 | 0 |
3300 +----++----+----+----+
3301 | b1 || 0 | b1 | 0 |
3302 +----++----+----+----+
3303 | b2 || 0 | 0 | b2 |
3304 +----++----+----+----+
3305 sage: HadamardEJA(3).multiplication_table()
3306 +----++----+----+----+
3307 | * || b0 | b1 | b2 |
3308 +====++====+====+====+
3309 | b0 || b0 | 0 | 0 |
3310 +----++----+----+----+
3311 | b1 || 0 | b1 | 0 |
3312 +----++----+----+----+
3313 | b2 || 0 | 0 | b2 |
3314 +----++----+----+----+
3318 All factors must share the same base field::
3320 sage: J1 = HadamardEJA(2, field=QQ)
3321 sage: J2 = RealSymmetricEJA(2)
3322 sage: CartesianProductEJA((J1,J2))
3323 Traceback (most recent call last):
3325 ValueError: all factors must share the same base field
3327 The cached unit element is the same one that would be computed::
3329 sage: set_random_seed() # long time
3330 sage: J1 = random_eja() # long time
3331 sage: J2 = random_eja() # long time
3332 sage: J = cartesian_product([J1,J2]) # long time
3333 sage: actual = J.one() # long time
3334 sage: J.one.clear_cache() # long time
3335 sage: expected = J.one() # long time
3336 sage: actual == expected # long time
3340 Element
= FiniteDimensionalEJAElement
3343 def __init__(self
, factors
, **kwargs
):
3348 self
._sets
= factors
3350 field
= factors
[0].base_ring()
3351 if not all( J
.base_ring() == field
for J
in factors
):
3352 raise ValueError("all factors must share the same base field")
3354 associative
= all( f
.is_associative() for f
in factors
)
3356 MS
= self
.matrix_space()
3360 for b
in factors
[i
].matrix_basis():
3365 basis
= tuple( MS(b
) for b
in basis
)
3367 # Define jordan/inner products that operate on that matrix_basis.
3368 def jordan_product(x
,y
):
3370 (factors
[i
](x
[i
])*factors
[i
](y
[i
])).to_matrix()
3374 def inner_product(x
, y
):
3376 factors
[i
](x
[i
]).inner_product(factors
[i
](y
[i
]))
3380 # There's no need to check the field since it already came
3381 # from an EJA. Likewise the axioms are guaranteed to be
3382 # satisfied, unless the guy writing this class sucks.
3384 # If you want the basis to be orthonormalized, orthonormalize
3386 FiniteDimensionalEJA
.__init
__(self
,
3391 orthonormalize
=False,
3392 associative
=associative
,
3393 cartesian_product
=True,
3397 ones
= tuple(J
.one().to_matrix() for J
in factors
)
3398 self
.one
.set_cache(self(ones
))
3399 self
.rank
.set_cache(sum(J
.rank() for J
in factors
))
3401 def cartesian_factors(self
):
3402 # Copy/pasted from CombinatorialFreeModule_CartesianProduct.
3405 def cartesian_factor(self
, i
):
3407 Return the ``i``th factor of this algebra.
3409 return self
._sets
[i
]
3412 # Copy/pasted from CombinatorialFreeModule_CartesianProduct.
3413 from sage
.categories
.cartesian_product
import cartesian_product
3414 return cartesian_product
.symbol
.join("%s" % factor
3415 for factor
in self
._sets
)
3417 def matrix_space(self
):
3419 Return the space that our matrix basis lives in as a Cartesian
3422 We don't simply use the ``cartesian_product()`` functor here
3423 because it acts differently on SageMath MatrixSpaces and our
3424 custom MatrixAlgebras, which are CombinatorialFreeModules. We
3425 always want the result to be represented (and indexed) as
3430 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
3432 ....: OctonionHermitianEJA,
3433 ....: RealSymmetricEJA)
3437 sage: J1 = HadamardEJA(1)
3438 sage: J2 = RealSymmetricEJA(2)
3439 sage: J = cartesian_product([J1,J2])
3440 sage: J.matrix_space()
3441 The Cartesian product of (Full MatrixSpace of 1 by 1 dense
3442 matrices over Algebraic Real Field, Full MatrixSpace of 2
3443 by 2 dense matrices over Algebraic Real Field)
3447 sage: J1 = ComplexHermitianEJA(1)
3448 sage: J2 = ComplexHermitianEJA(1)
3449 sage: J = cartesian_product([J1,J2])
3450 sage: J.one().to_matrix()[0]
3453 sage: J.one().to_matrix()[1]
3459 sage: J1 = OctonionHermitianEJA(1)
3460 sage: J2 = OctonionHermitianEJA(1)
3461 sage: J = cartesian_product([J1,J2])
3462 sage: J.one().to_matrix()[0]
3466 sage: J.one().to_matrix()[1]
3472 scalars
= self
.cartesian_factor(0).base_ring()
3474 # This category isn't perfect, but is good enough for what we
3476 cat
= MagmaticAlgebras(scalars
).FiniteDimensional().WithBasis()
3477 cat
= cat
.Unital().CartesianProducts()
3478 factors
= tuple( J
.matrix_space() for J
in self
.cartesian_factors() )
3480 from sage
.sets
.cartesian_product
import CartesianProduct
3481 return CartesianProduct(factors
, cat
)
3485 def cartesian_projection(self
, i
):
3489 sage: from mjo.eja.eja_algebra import (random_eja,
3490 ....: JordanSpinEJA,
3492 ....: RealSymmetricEJA,
3493 ....: ComplexHermitianEJA)
3497 The projection morphisms are Euclidean Jordan algebra
3500 sage: J1 = HadamardEJA(2)
3501 sage: J2 = RealSymmetricEJA(2)
3502 sage: J = cartesian_product([J1,J2])
3503 sage: J.cartesian_projection(0)
3504 Linear operator between finite-dimensional Euclidean Jordan
3505 algebras represented by the matrix:
3508 Domain: Euclidean Jordan algebra of dimension 2 over Algebraic
3509 Real Field (+) Euclidean Jordan algebra of dimension 3 over
3510 Algebraic Real Field
3511 Codomain: Euclidean Jordan algebra of dimension 2 over Algebraic
3513 sage: J.cartesian_projection(1)
3514 Linear operator between finite-dimensional Euclidean Jordan
3515 algebras represented by the matrix:
3519 Domain: Euclidean Jordan algebra of dimension 2 over Algebraic
3520 Real Field (+) Euclidean Jordan algebra of dimension 3 over
3521 Algebraic Real Field
3522 Codomain: Euclidean Jordan algebra of dimension 3 over Algebraic
3525 The projections work the way you'd expect on the vector
3526 representation of an element::
3528 sage: J1 = JordanSpinEJA(2)
3529 sage: J2 = ComplexHermitianEJA(2)
3530 sage: J = cartesian_product([J1,J2])
3531 sage: pi_left = J.cartesian_projection(0)
3532 sage: pi_right = J.cartesian_projection(1)
3533 sage: pi_left(J.one()).to_vector()
3535 sage: pi_right(J.one()).to_vector()
3537 sage: J.one().to_vector()
3542 The answer never changes::
3544 sage: set_random_seed()
3545 sage: J1 = random_eja()
3546 sage: J2 = random_eja()
3547 sage: J = cartesian_product([J1,J2])
3548 sage: P0 = J.cartesian_projection(0)
3549 sage: P1 = J.cartesian_projection(0)
3554 offset
= sum( self
.cartesian_factor(k
).dimension()
3556 Ji
= self
.cartesian_factor(i
)
3557 Pi
= self
._module
_morphism
(lambda j
: Ji
.monomial(j
- offset
),
3560 return FiniteDimensionalEJAOperator(self
,Ji
,Pi
.matrix())
3563 def cartesian_embedding(self
, i
):
3567 sage: from mjo.eja.eja_algebra import (random_eja,
3568 ....: JordanSpinEJA,
3570 ....: RealSymmetricEJA)
3574 The embedding morphisms are Euclidean Jordan algebra
3577 sage: J1 = HadamardEJA(2)
3578 sage: J2 = RealSymmetricEJA(2)
3579 sage: J = cartesian_product([J1,J2])
3580 sage: J.cartesian_embedding(0)
3581 Linear operator between finite-dimensional Euclidean Jordan
3582 algebras represented by the matrix:
3588 Domain: Euclidean Jordan algebra of dimension 2 over
3589 Algebraic Real Field
3590 Codomain: Euclidean Jordan algebra of dimension 2 over
3591 Algebraic Real Field (+) Euclidean Jordan algebra of
3592 dimension 3 over Algebraic Real Field
3593 sage: J.cartesian_embedding(1)
3594 Linear operator between finite-dimensional Euclidean Jordan
3595 algebras represented by the matrix:
3601 Domain: Euclidean Jordan algebra of dimension 3 over
3602 Algebraic Real Field
3603 Codomain: Euclidean Jordan algebra of dimension 2 over
3604 Algebraic Real Field (+) Euclidean Jordan algebra of
3605 dimension 3 over Algebraic Real Field
3607 The embeddings work the way you'd expect on the vector
3608 representation of an element::
3610 sage: J1 = JordanSpinEJA(3)
3611 sage: J2 = RealSymmetricEJA(2)
3612 sage: J = cartesian_product([J1,J2])
3613 sage: iota_left = J.cartesian_embedding(0)
3614 sage: iota_right = J.cartesian_embedding(1)
3615 sage: iota_left(J1.zero()) == J.zero()
3617 sage: iota_right(J2.zero()) == J.zero()
3619 sage: J1.one().to_vector()
3621 sage: iota_left(J1.one()).to_vector()
3623 sage: J2.one().to_vector()
3625 sage: iota_right(J2.one()).to_vector()
3627 sage: J.one().to_vector()
3632 The answer never changes::
3634 sage: set_random_seed()
3635 sage: J1 = random_eja()
3636 sage: J2 = random_eja()
3637 sage: J = cartesian_product([J1,J2])
3638 sage: E0 = J.cartesian_embedding(0)
3639 sage: E1 = J.cartesian_embedding(0)
3643 Composing a projection with the corresponding inclusion should
3644 produce the identity map, and mismatching them should produce
3647 sage: set_random_seed()
3648 sage: J1 = random_eja()
3649 sage: J2 = random_eja()
3650 sage: J = cartesian_product([J1,J2])
3651 sage: iota_left = J.cartesian_embedding(0)
3652 sage: iota_right = J.cartesian_embedding(1)
3653 sage: pi_left = J.cartesian_projection(0)
3654 sage: pi_right = J.cartesian_projection(1)
3655 sage: pi_left*iota_left == J1.one().operator()
3657 sage: pi_right*iota_right == J2.one().operator()
3659 sage: (pi_left*iota_right).is_zero()
3661 sage: (pi_right*iota_left).is_zero()
3665 offset
= sum( self
.cartesian_factor(k
).dimension()
3667 Ji
= self
.cartesian_factor(i
)
3668 Ei
= Ji
._module
_morphism
(lambda j
: self
.monomial(j
+ offset
),
3670 return FiniteDimensionalEJAOperator(Ji
,self
,Ei
.matrix())
3674 FiniteDimensionalEJA
.CartesianProduct
= CartesianProductEJA
3676 class RationalBasisCartesianProductEJA(CartesianProductEJA
,
3679 A separate class for products of algebras for which we know a
3684 sage: from mjo.eja.eja_algebra import (HadamardEJA,
3685 ....: JordanSpinEJA,
3686 ....: OctonionHermitianEJA,
3687 ....: RealSymmetricEJA)
3691 This gives us fast characteristic polynomial computations in
3692 product algebras, too::
3695 sage: J1 = JordanSpinEJA(2)
3696 sage: J2 = RealSymmetricEJA(3)
3697 sage: J = cartesian_product([J1,J2])
3698 sage: J.characteristic_polynomial_of().degree()
3705 The ``cartesian_product()`` function only uses the first factor to
3706 decide where the result will live; thus we have to be careful to
3707 check that all factors do indeed have a `_rational_algebra` member
3708 before we try to access it::
3710 sage: J1 = OctonionHermitianEJA(1) # no rational basis
3711 sage: J2 = HadamardEJA(2)
3712 sage: cartesian_product([J1,J2])
3713 Euclidean Jordan algebra of dimension 1 over Algebraic Real Field
3714 (+) Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
3715 sage: cartesian_product([J2,J1])
3716 Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
3717 (+) Euclidean Jordan algebra of dimension 1 over Algebraic Real Field
3720 def __init__(self
, algebras
, **kwargs
):
3721 CartesianProductEJA
.__init
__(self
, algebras
, **kwargs
)
3723 self
._rational
_algebra
= None
3724 if self
.vector_space().base_field() is not QQ
:
3725 if all( hasattr(r
, "_rational_algebra") for r
in algebras
):
3726 self
._rational
_algebra
= cartesian_product([
3727 r
._rational
_algebra
for r
in algebras
3731 RationalBasisEJA
.CartesianProduct
= RationalBasisCartesianProductEJA
3733 def random_eja(*args
, **kwargs
):
3734 J1
= ConcreteEJA
.random_instance(*args
, **kwargs
)
3736 # This might make Cartesian products appear roughly as often as
3737 # any other ConcreteEJA.
3738 if ZZ
.random_element(len(ConcreteEJA
.__subclasses
__()) + 1) == 0:
3739 # Use random_eja() again so we can get more than two factors.
3740 J2
= random_eja(*args
, **kwargs
)
3741 J
= cartesian_product([J1
,J2
])