2 Representations and constructions for Euclidean Jordan algebras.
4 A Euclidean Jordan algebra is a Jordan algebra that has some
7 1. It is finite-dimensional.
8 2. Its scalar field is the real numbers.
9 3a. An inner product is defined on it, and...
10 3b. That inner product is compatible with the Jordan product
11 in the sense that `<x*y,z> = <y,x*z>` for all elements
12 `x,y,z` in the algebra.
14 Every Euclidean Jordan algebra is formally-real: for any two elements
15 `x` and `y` in the algebra, `x^{2} + y^{2} = 0` implies that `x = y =
16 0`. Conversely, every finite-dimensional formally-real Jordan algebra
17 can be made into a Euclidean Jordan algebra with an appropriate choice
20 Formally-real Jordan algebras were originally studied as a framework
21 for quantum mechanics. Today, Euclidean Jordan algebras are crucial in
22 symmetric cone optimization, since every symmetric cone arises as the
23 cone of squares in some Euclidean Jordan algebra.
25 It is known that every Euclidean Jordan algebra decomposes into an
26 orthogonal direct sum (essentially, a Cartesian product) of simple
27 algebras, and that moreover, up to Jordan-algebra isomorphism, there
28 are only five families of simple algebras. We provide constructions
29 for these simple algebras:
31 * :class:`BilinearFormEJA`
32 * :class:`RealSymmetricEJA`
33 * :class:`ComplexHermitianEJA`
34 * :class:`QuaternionHermitianEJA`
36 Missing from this list is the algebra of three-by-three octononion
37 Hermitian matrices, as there is (as of yet) no implementation of the
38 octonions in SageMath. In addition to these, we provide two other
39 example constructions,
41 * :class:`HadamardEJA`
44 The Jordan spin algebra is a bilinear form algebra where the bilinear
45 form is the identity. The Hadamard EJA is simply a Cartesian product
46 of one-dimensional spin algebras. And last but not least, the trivial
47 EJA is exactly what you think. Cartesian products of these are also
48 supported using the usual ``cartesian_product()`` function; as a
49 result, we support (up to isomorphism) all Euclidean Jordan algebras
50 that don't involve octonions.
54 sage: from mjo.eja.eja_algebra import random_eja
59 Euclidean Jordan algebra of dimension...
62 from itertools
import repeat
64 from sage
.algebras
.quatalg
.quaternion_algebra
import QuaternionAlgebra
65 from sage
.categories
.magmatic_algebras
import MagmaticAlgebras
66 from sage
.categories
.sets_cat
import cartesian_product
67 from sage
.combinat
.free_module
import CombinatorialFreeModule
68 from sage
.matrix
.constructor
import matrix
69 from sage
.matrix
.matrix_space
import MatrixSpace
70 from sage
.misc
.cachefunc
import cached_method
71 from sage
.misc
.table
import table
72 from sage
.modules
.free_module
import FreeModule
, VectorSpace
73 from sage
.rings
.all
import (ZZ
, QQ
, AA
, QQbar
, RR
, RLF
, CLF
,
76 from mjo
.eja
.eja_element
import FiniteDimensionalEJAElement
77 from mjo
.eja
.eja_operator
import FiniteDimensionalEJAOperator
78 from mjo
.eja
.eja_utils
import _all2list
, _mat2vec
80 class FiniteDimensionalEJA(CombinatorialFreeModule
):
82 A finite-dimensional Euclidean Jordan algebra.
86 - ``basis`` -- a tuple; a tuple of basis elements in "matrix
87 form," which must be the same form as the arguments to
88 ``jordan_product`` and ``inner_product``. In reality, "matrix
89 form" can be either vectors, matrices, or a Cartesian product
90 (ordered tuple) of vectors or matrices. All of these would
91 ideally be vector spaces in sage with no special-casing
92 needed; but in reality we turn vectors into column-matrices
93 and Cartesian products `(a,b)` into column matrices
94 `(a,b)^{T}` after converting `a` and `b` themselves.
96 - ``jordan_product`` -- a function; afunction of two ``basis``
97 elements (in matrix form) that returns their jordan product,
98 also in matrix form; this will be applied to ``basis`` to
99 compute a multiplication table for the algebra.
101 - ``inner_product`` -- a function; a function of two ``basis``
102 elements (in matrix form) that returns their inner
103 product. This will be applied to ``basis`` to compute an
104 inner-product table (basically a matrix) for this algebra.
106 - ``field`` -- a subfield of the reals (default: ``AA``); the scalar
107 field for the algebra.
109 - ``orthonormalize`` -- boolean (default: ``True``); whether or
110 not to orthonormalize the basis. Doing so is expensive and
111 generally rules out using the rationals as your ``field``, but
112 is required for spectral decompositions.
116 sage: from mjo.eja.eja_algebra import random_eja
120 We should compute that an element subalgebra is associative even
121 if we circumvent the element method::
123 sage: set_random_seed()
124 sage: J = random_eja(field=QQ,orthonormalize=False)
125 sage: x = J.random_element()
126 sage: A = x.subalgebra_generated_by(orthonormalize=False)
127 sage: basis = tuple(b.superalgebra_element() for b in A.basis())
128 sage: J.subalgebra(basis, orthonormalize=False).is_associative()
132 Element
= FiniteDimensionalEJAElement
141 cartesian_product
=False,
149 if not field
.is_subring(RR
):
150 # Note: this does return true for the real algebraic
151 # field, the rationals, and any quadratic field where
152 # we've specified a real embedding.
153 raise ValueError("scalar field is not real")
156 # Check commutativity of the Jordan and inner-products.
157 # This has to be done before we build the multiplication
158 # and inner-product tables/matrices, because we take
159 # advantage of symmetry in the process.
160 if not all( jordan_product(bi
,bj
) == jordan_product(bj
,bi
)
163 raise ValueError("Jordan product is not commutative")
165 if not all( inner_product(bi
,bj
) == inner_product(bj
,bi
)
168 raise ValueError("inner-product is not commutative")
171 category
= MagmaticAlgebras(field
).FiniteDimensional()
172 category
= category
.WithBasis().Unital().Commutative()
174 if associative
is None:
175 # We should figure it out. As with check_axioms, we have to do
176 # this without the help of the _jordan_product_is_associative()
177 # method because we need to know the category before we
178 # initialize the algebra.
179 associative
= all( jordan_product(jordan_product(bi
,bj
),bk
)
181 jordan_product(bi
,jordan_product(bj
,bk
))
187 # Element subalgebras can take advantage of this.
188 category
= category
.Associative()
189 if cartesian_product
:
190 # Use join() here because otherwise we only get the
191 # "Cartesian product of..." and not the things themselves.
192 category
= category
.join([category
,
193 category
.CartesianProducts()])
195 # Call the superclass constructor so that we can use its from_vector()
196 # method to build our multiplication table.
197 CombinatorialFreeModule
.__init
__(self
,
204 # Now comes all of the hard work. We'll be constructing an
205 # ambient vector space V that our (vectorized) basis lives in,
206 # as well as a subspace W of V spanned by those (vectorized)
207 # basis elements. The W-coordinates are the coefficients that
208 # we see in things like x = 1*b1 + 2*b2.
213 degree
= len(_all2list(basis
[0]))
215 # Build an ambient space that fits our matrix basis when
216 # written out as "long vectors."
217 V
= VectorSpace(field
, degree
)
219 # The matrix that will hole the orthonormal -> unorthonormal
220 # coordinate transformation.
221 self
._deortho
_matrix
= None
224 # Save a copy of the un-orthonormalized basis for later.
225 # Convert it to ambient V (vector) coordinates while we're
226 # at it, because we'd have to do it later anyway.
227 deortho_vector_basis
= tuple( V(_all2list(b
)) for b
in basis
)
229 from mjo
.eja
.eja_utils
import gram_schmidt
230 basis
= tuple(gram_schmidt(basis
, inner_product
))
232 # Save the (possibly orthonormalized) matrix basis for
234 self
._matrix
_basis
= basis
236 # Now create the vector space for the algebra, which will have
237 # its own set of non-ambient coordinates (in terms of the
239 vector_basis
= tuple( V(_all2list(b
)) for b
in basis
)
240 W
= V
.span_of_basis( vector_basis
, check
=check_axioms
)
243 # Now "W" is the vector space of our algebra coordinates. The
244 # variables "X1", "X2",... refer to the entries of vectors in
245 # W. Thus to convert back and forth between the orthonormal
246 # coordinates and the given ones, we need to stick the original
248 U
= V
.span_of_basis( deortho_vector_basis
, check
=check_axioms
)
249 self
._deortho
_matrix
= matrix( U
.coordinate_vector(q
)
250 for q
in vector_basis
)
253 # Now we actually compute the multiplication and inner-product
254 # tables/matrices using the possibly-orthonormalized basis.
255 self
._inner
_product
_matrix
= matrix
.identity(field
, n
)
256 self
._multiplication
_table
= [ [0 for j
in range(i
+1)]
259 # Note: the Jordan and inner-products are defined in terms
260 # of the ambient basis. It's important that their arguments
261 # are in ambient coordinates as well.
264 # ortho basis w.r.t. ambient coords
268 # The jordan product returns a matrixy answer, so we
269 # have to convert it to the algebra coordinates.
270 elt
= jordan_product(q_i
, q_j
)
271 elt
= W
.coordinate_vector(V(_all2list(elt
)))
272 self
._multiplication
_table
[i
][j
] = self
.from_vector(elt
)
274 if not orthonormalize
:
275 # If we're orthonormalizing the basis with respect
276 # to an inner-product, then the inner-product
277 # matrix with respect to the resulting basis is
278 # just going to be the identity.
279 ip
= inner_product(q_i
, q_j
)
280 self
._inner
_product
_matrix
[i
,j
] = ip
281 self
._inner
_product
_matrix
[j
,i
] = ip
283 self
._inner
_product
_matrix
._cache
= {'hermitian': True}
284 self
._inner
_product
_matrix
.set_immutable()
287 if not self
._is
_jordanian
():
288 raise ValueError("Jordan identity does not hold")
289 if not self
._inner
_product
_is
_associative
():
290 raise ValueError("inner product is not associative")
293 def _coerce_map_from_base_ring(self
):
295 Disable the map from the base ring into the algebra.
297 Performing a nonsense conversion like this automatically
298 is counterpedagogical. The fallback is to try the usual
299 element constructor, which should also fail.
303 sage: from mjo.eja.eja_algebra import random_eja
307 sage: set_random_seed()
308 sage: J = random_eja()
310 Traceback (most recent call last):
312 ValueError: not an element of this algebra
318 def product_on_basis(self
, i
, j
):
320 Returns the Jordan product of the `i` and `j`th basis elements.
322 This completely defines the Jordan product on the algebra, and
323 is used direclty by our superclass machinery to implement
328 sage: from mjo.eja.eja_algebra import random_eja
332 sage: set_random_seed()
333 sage: J = random_eja()
334 sage: n = J.dimension()
337 sage: bi_bj = J.zero()*J.zero()
339 ....: i = ZZ.random_element(n)
340 ....: j = ZZ.random_element(n)
341 ....: bi = J.monomial(i)
342 ....: bj = J.monomial(j)
343 ....: bi_bj = J.product_on_basis(i,j)
348 # We only stored the lower-triangular portion of the
349 # multiplication table.
351 return self
._multiplication
_table
[i
][j
]
353 return self
._multiplication
_table
[j
][i
]
355 def inner_product(self
, x
, y
):
357 The inner product associated with this Euclidean Jordan algebra.
359 Defaults to the trace inner product, but can be overridden by
360 subclasses if they are sure that the necessary properties are
365 sage: from mjo.eja.eja_algebra import (random_eja,
367 ....: BilinearFormEJA)
371 Our inner product is "associative," which means the following for
372 a symmetric bilinear form::
374 sage: set_random_seed()
375 sage: J = random_eja()
376 sage: x,y,z = J.random_elements(3)
377 sage: (x*y).inner_product(z) == y.inner_product(x*z)
382 Ensure that this is the usual inner product for the algebras
385 sage: set_random_seed()
386 sage: J = HadamardEJA.random_instance()
387 sage: x,y = J.random_elements(2)
388 sage: actual = x.inner_product(y)
389 sage: expected = x.to_vector().inner_product(y.to_vector())
390 sage: actual == expected
393 Ensure that this is one-half of the trace inner-product in a
394 BilinearFormEJA that isn't just the reals (when ``n`` isn't
395 one). This is in Faraut and Koranyi, and also my "On the
398 sage: set_random_seed()
399 sage: J = BilinearFormEJA.random_instance()
400 sage: n = J.dimension()
401 sage: x = J.random_element()
402 sage: y = J.random_element()
403 sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
407 B
= self
._inner
_product
_matrix
408 return (B
*x
.to_vector()).inner_product(y
.to_vector())
411 def is_associative(self
):
413 Return whether or not this algebra's Jordan product is associative.
417 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
421 sage: J = ComplexHermitianEJA(3, field=QQ, orthonormalize=False)
422 sage: J.is_associative()
424 sage: x = sum(J.gens())
425 sage: A = x.subalgebra_generated_by(orthonormalize=False)
426 sage: A.is_associative()
430 return "Associative" in self
.category().axioms()
432 def _is_commutative(self
):
434 Whether or not this algebra's multiplication table is commutative.
436 This method should of course always return ``True``, unless
437 this algebra was constructed with ``check_axioms=False`` and
438 passed an invalid multiplication table.
440 return all( x
*y
== y
*x
for x
in self
.gens() for y
in self
.gens() )
442 def _is_jordanian(self
):
444 Whether or not this algebra's multiplication table respects the
445 Jordan identity `(x^{2})(xy) = x(x^{2}y)`.
447 We only check one arrangement of `x` and `y`, so for a
448 ``True`` result to be truly true, you should also check
449 :meth:`_is_commutative`. This method should of course always
450 return ``True``, unless this algebra was constructed with
451 ``check_axioms=False`` and passed an invalid multiplication table.
453 return all( (self
.monomial(i
)**2)*(self
.monomial(i
)*self
.monomial(j
))
455 (self
.monomial(i
))*((self
.monomial(i
)**2)*self
.monomial(j
))
456 for i
in range(self
.dimension())
457 for j
in range(self
.dimension()) )
459 def _jordan_product_is_associative(self
):
461 Return whether or not this algebra's Jordan product is
462 associative; that is, whether or not `x*(y*z) = (x*y)*z`
465 This method should agree with :meth:`is_associative` unless
466 you lied about the value of the ``associative`` parameter
467 when you constructed the algebra.
471 sage: from mjo.eja.eja_algebra import (random_eja,
472 ....: RealSymmetricEJA,
473 ....: ComplexHermitianEJA,
474 ....: QuaternionHermitianEJA)
478 sage: J = RealSymmetricEJA(4, orthonormalize=False)
479 sage: J._jordan_product_is_associative()
481 sage: x = sum(J.gens())
482 sage: A = x.subalgebra_generated_by()
483 sage: A._jordan_product_is_associative()
488 sage: J = ComplexHermitianEJA(2,field=QQ,orthonormalize=False)
489 sage: J._jordan_product_is_associative()
491 sage: x = sum(J.gens())
492 sage: A = x.subalgebra_generated_by(orthonormalize=False)
493 sage: A._jordan_product_is_associative()
498 sage: J = QuaternionHermitianEJA(2)
499 sage: J._jordan_product_is_associative()
501 sage: x = sum(J.gens())
502 sage: A = x.subalgebra_generated_by()
503 sage: A._jordan_product_is_associative()
508 The values we've presupplied to the constructors agree with
511 sage: set_random_seed()
512 sage: J = random_eja()
513 sage: J.is_associative() == J._jordan_product_is_associative()
519 # Used to check whether or not something is zero.
522 # I don't know of any examples that make this magnitude
523 # necessary because I don't know how to make an
524 # associative algebra when the element subalgebra
525 # construction is unreliable (as it is over RDF; we can't
526 # find the degree of an element because we can't compute
527 # the rank of a matrix). But even multiplication of floats
528 # is non-associative, so *some* epsilon is needed... let's
529 # just take the one from _inner_product_is_associative?
532 for i
in range(self
.dimension()):
533 for j
in range(self
.dimension()):
534 for k
in range(self
.dimension()):
538 diff
= (x
*y
)*z
- x
*(y
*z
)
540 if diff
.norm() > epsilon
:
545 def _inner_product_is_associative(self
):
547 Return whether or not this algebra's inner product `B` is
548 associative; that is, whether or not `B(xy,z) = B(x,yz)`.
550 This method should of course always return ``True``, unless
551 this algebra was constructed with ``check_axioms=False`` and
552 passed an invalid Jordan or inner-product.
556 # Used to check whether or not something is zero.
559 # This choice is sufficient to allow the construction of
560 # QuaternionHermitianEJA(2, field=RDF) with check_axioms=True.
563 for i
in range(self
.dimension()):
564 for j
in range(self
.dimension()):
565 for k
in range(self
.dimension()):
569 diff
= (x
*y
).inner_product(z
) - x
.inner_product(y
*z
)
571 if diff
.abs() > epsilon
:
576 def _element_constructor_(self
, elt
):
578 Construct an element of this algebra from its vector or matrix
581 This gets called only after the parent element _call_ method
582 fails to find a coercion for the argument.
586 sage: from mjo.eja.eja_algebra import (random_eja,
589 ....: RealSymmetricEJA)
593 The identity in `S^n` is converted to the identity in the EJA::
595 sage: J = RealSymmetricEJA(3)
596 sage: I = matrix.identity(QQ,3)
597 sage: J(I) == J.one()
600 This skew-symmetric matrix can't be represented in the EJA::
602 sage: J = RealSymmetricEJA(3)
603 sage: A = matrix(QQ,3, lambda i,j: i-j)
605 Traceback (most recent call last):
607 ValueError: not an element of this algebra
609 Tuples work as well, provided that the matrix basis for the
610 algebra consists of them::
612 sage: J1 = HadamardEJA(3)
613 sage: J2 = RealSymmetricEJA(2)
614 sage: J = cartesian_product([J1,J2])
615 sage: J( (J1.matrix_basis()[1], J2.matrix_basis()[2]) )
620 Ensure that we can convert any element back and forth
621 faithfully between its matrix and algebra representations::
623 sage: set_random_seed()
624 sage: J = random_eja()
625 sage: x = J.random_element()
626 sage: J(x.to_matrix()) == x
629 We cannot coerce elements between algebras just because their
630 matrix representations are compatible::
632 sage: J1 = HadamardEJA(3)
633 sage: J2 = JordanSpinEJA(3)
635 Traceback (most recent call last):
637 ValueError: not an element of this algebra
639 Traceback (most recent call last):
641 ValueError: not an element of this algebra
643 msg
= "not an element of this algebra"
644 if elt
in self
.base_ring():
645 # Ensure that no base ring -> algebra coercion is performed
646 # by this method. There's some stupidity in sage that would
647 # otherwise propagate to this method; for example, sage thinks
648 # that the integer 3 belongs to the space of 2-by-2 matrices.
649 raise ValueError(msg
)
652 # Try to convert a vector into a column-matrix...
654 except (AttributeError, TypeError):
655 # and ignore failure, because we weren't really expecting
656 # a vector as an argument anyway.
659 if elt
not in self
.matrix_space():
660 raise ValueError(msg
)
662 # Thanks for nothing! Matrix spaces aren't vector spaces in
663 # Sage, so we have to figure out its matrix-basis coordinates
664 # ourselves. We use the basis space's ring instead of the
665 # element's ring because the basis space might be an algebraic
666 # closure whereas the base ring of the 3-by-3 identity matrix
667 # could be QQ instead of QQbar.
669 # And, we also have to handle Cartesian product bases (when
670 # the matrix basis consists of tuples) here. The "good news"
671 # is that we're already converting everything to long vectors,
672 # and that strategy works for tuples as well.
674 # We pass check=False because the matrix basis is "guaranteed"
675 # to be linearly independent... right? Ha ha.
677 V
= VectorSpace(self
.base_ring(), len(elt
))
678 W
= V
.span_of_basis( (V(_all2list(s
)) for s
in self
.matrix_basis()),
682 coords
= W
.coordinate_vector(V(elt
))
683 except ArithmeticError: # vector is not in free module
684 raise ValueError(msg
)
686 return self
.from_vector(coords
)
690 Return a string representation of ``self``.
694 sage: from mjo.eja.eja_algebra import JordanSpinEJA
698 Ensure that it says what we think it says::
700 sage: JordanSpinEJA(2, field=AA)
701 Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
702 sage: JordanSpinEJA(3, field=RDF)
703 Euclidean Jordan algebra of dimension 3 over Real Double Field
706 fmt
= "Euclidean Jordan algebra of dimension {} over {}"
707 return fmt
.format(self
.dimension(), self
.base_ring())
711 def characteristic_polynomial_of(self
):
713 Return the algebra's "characteristic polynomial of" function,
714 which is itself a multivariate polynomial that, when evaluated
715 at the coordinates of some algebra element, returns that
716 element's characteristic polynomial.
718 The resulting polynomial has `n+1` variables, where `n` is the
719 dimension of this algebra. The first `n` variables correspond to
720 the coordinates of an algebra element: when evaluated at the
721 coordinates of an algebra element with respect to a certain
722 basis, the result is a univariate polynomial (in the one
723 remaining variable ``t``), namely the characteristic polynomial
728 sage: from mjo.eja.eja_algebra import JordanSpinEJA, TrivialEJA
732 The characteristic polynomial in the spin algebra is given in
733 Alizadeh, Example 11.11::
735 sage: J = JordanSpinEJA(3)
736 sage: p = J.characteristic_polynomial_of(); p
737 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
738 sage: xvec = J.one().to_vector()
742 By definition, the characteristic polynomial is a monic
743 degree-zero polynomial in a rank-zero algebra. Note that
744 Cayley-Hamilton is indeed satisfied since the polynomial
745 ``1`` evaluates to the identity element of the algebra on
748 sage: J = TrivialEJA()
749 sage: J.characteristic_polynomial_of()
756 # The list of coefficient polynomials a_0, a_1, a_2, ..., a_(r-1).
757 a
= self
._charpoly
_coefficients
()
759 # We go to a bit of trouble here to reorder the
760 # indeterminates, so that it's easier to evaluate the
761 # characteristic polynomial at x's coordinates and get back
762 # something in terms of t, which is what we want.
763 S
= PolynomialRing(self
.base_ring(),'t')
767 S
= PolynomialRing(S
, R
.variable_names())
770 return (t
**r
+ sum( a
[k
]*(t
**k
) for k
in range(r
) ))
772 def coordinate_polynomial_ring(self
):
774 The multivariate polynomial ring in which this algebra's
775 :meth:`characteristic_polynomial_of` lives.
779 sage: from mjo.eja.eja_algebra import (HadamardEJA,
780 ....: RealSymmetricEJA)
784 sage: J = HadamardEJA(2)
785 sage: J.coordinate_polynomial_ring()
786 Multivariate Polynomial Ring in X1, X2...
787 sage: J = RealSymmetricEJA(3,field=QQ,orthonormalize=False)
788 sage: J.coordinate_polynomial_ring()
789 Multivariate Polynomial Ring in X1, X2, X3, X4, X5, X6...
792 var_names
= tuple( "X%d" % z
for z
in range(1, self
.dimension()+1) )
793 return PolynomialRing(self
.base_ring(), var_names
)
795 def inner_product(self
, x
, y
):
797 The inner product associated with this Euclidean Jordan algebra.
799 Defaults to the trace inner product, but can be overridden by
800 subclasses if they are sure that the necessary properties are
805 sage: from mjo.eja.eja_algebra import (random_eja,
807 ....: BilinearFormEJA)
811 Our inner product is "associative," which means the following for
812 a symmetric bilinear form::
814 sage: set_random_seed()
815 sage: J = random_eja()
816 sage: x,y,z = J.random_elements(3)
817 sage: (x*y).inner_product(z) == y.inner_product(x*z)
822 Ensure that this is the usual inner product for the algebras
825 sage: set_random_seed()
826 sage: J = HadamardEJA.random_instance()
827 sage: x,y = J.random_elements(2)
828 sage: actual = x.inner_product(y)
829 sage: expected = x.to_vector().inner_product(y.to_vector())
830 sage: actual == expected
833 Ensure that this is one-half of the trace inner-product in a
834 BilinearFormEJA that isn't just the reals (when ``n`` isn't
835 one). This is in Faraut and Koranyi, and also my "On the
838 sage: set_random_seed()
839 sage: J = BilinearFormEJA.random_instance()
840 sage: n = J.dimension()
841 sage: x = J.random_element()
842 sage: y = J.random_element()
843 sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
846 B
= self
._inner
_product
_matrix
847 return (B
*x
.to_vector()).inner_product(y
.to_vector())
850 def is_trivial(self
):
852 Return whether or not this algebra is trivial.
854 A trivial algebra contains only the zero element.
858 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
863 sage: J = ComplexHermitianEJA(3)
869 sage: J = TrivialEJA()
874 return self
.dimension() == 0
877 def multiplication_table(self
):
879 Return a visual representation of this algebra's multiplication
880 table (on basis elements).
884 sage: from mjo.eja.eja_algebra import JordanSpinEJA
888 sage: J = JordanSpinEJA(4)
889 sage: J.multiplication_table()
890 +----++----+----+----+----+
891 | * || b0 | b1 | b2 | b3 |
892 +====++====+====+====+====+
893 | b0 || b0 | b1 | b2 | b3 |
894 +----++----+----+----+----+
895 | b1 || b1 | b0 | 0 | 0 |
896 +----++----+----+----+----+
897 | b2 || b2 | 0 | b0 | 0 |
898 +----++----+----+----+----+
899 | b3 || b3 | 0 | 0 | b0 |
900 +----++----+----+----+----+
904 # Prepend the header row.
905 M
= [["*"] + list(self
.gens())]
907 # And to each subsequent row, prepend an entry that belongs to
908 # the left-side "header column."
909 M
+= [ [self
.monomial(i
)] + [ self
.monomial(i
)*self
.monomial(j
)
913 return table(M
, header_row
=True, header_column
=True, frame
=True)
916 def matrix_basis(self
):
918 Return an (often more natural) representation of this algebras
919 basis as an ordered tuple of matrices.
921 Every finite-dimensional Euclidean Jordan Algebra is a, up to
922 Jordan isomorphism, a direct sum of five simple
923 algebras---four of which comprise Hermitian matrices. And the
924 last type of algebra can of course be thought of as `n`-by-`1`
925 column matrices (ambiguusly called column vectors) to avoid
926 special cases. As a result, matrices (and column vectors) are
927 a natural representation format for Euclidean Jordan algebra
930 But, when we construct an algebra from a basis of matrices,
931 those matrix representations are lost in favor of coordinate
932 vectors *with respect to* that basis. We could eventually
933 convert back if we tried hard enough, but having the original
934 representations handy is valuable enough that we simply store
935 them and return them from this method.
937 Why implement this for non-matrix algebras? Avoiding special
938 cases for the :class:`BilinearFormEJA` pays with simplicity in
939 its own right. But mainly, we would like to be able to assume
940 that elements of a :class:`CartesianProductEJA` can be displayed
941 nicely, without having to have special classes for direct sums
942 one of whose components was a matrix algebra.
946 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
947 ....: RealSymmetricEJA)
951 sage: J = RealSymmetricEJA(2)
953 Finite family {0: b0, 1: b1, 2: b2}
954 sage: J.matrix_basis()
956 [1 0] [ 0 0.7071067811865475?] [0 0]
957 [0 0], [0.7071067811865475? 0], [0 1]
962 sage: J = JordanSpinEJA(2)
964 Finite family {0: b0, 1: b1}
965 sage: J.matrix_basis()
971 return self
._matrix
_basis
974 def matrix_space(self
):
976 Return the matrix space in which this algebra's elements live, if
977 we think of them as matrices (including column vectors of the
980 "By default" this will be an `n`-by-`1` column-matrix space,
981 except when the algebra is trivial. There it's `n`-by-`n`
982 (where `n` is zero), to ensure that two elements of the matrix
983 space (empty matrices) can be multiplied. For algebras of
984 matrices, this returns the space in which their
985 real embeddings live.
989 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
991 ....: QuaternionHermitianEJA,
996 By default, the matrix representation is just a column-matrix
997 equivalent to the vector representation::
999 sage: J = JordanSpinEJA(3)
1000 sage: J.matrix_space()
1001 Full MatrixSpace of 3 by 1 dense matrices over Algebraic
1004 The matrix representation in the trivial algebra is
1005 zero-by-zero instead of the usual `n`-by-one::
1007 sage: J = TrivialEJA()
1008 sage: J.matrix_space()
1009 Full MatrixSpace of 0 by 0 dense matrices over Algebraic
1012 The matrix space for complex/quaternion Hermitian matrix EJA
1013 is the space in which their real-embeddings live, not the
1014 original complex/quaternion matrix space::
1016 sage: J = ComplexHermitianEJA(2,field=QQ,orthonormalize=False)
1017 sage: J.matrix_space()
1018 Full MatrixSpace of 4 by 4 dense matrices over Rational Field
1019 sage: J = QuaternionHermitianEJA(1,field=QQ,orthonormalize=False)
1020 sage: J.matrix_space()
1021 Full MatrixSpace of 4 by 4 dense matrices over Rational Field
1024 if self
.is_trivial():
1025 return MatrixSpace(self
.base_ring(), 0)
1027 return self
.matrix_basis()[0].parent()
1033 Return the unit element of this algebra.
1037 sage: from mjo.eja.eja_algebra import (HadamardEJA,
1042 We can compute unit element in the Hadamard EJA::
1044 sage: J = HadamardEJA(5)
1046 b0 + b1 + b2 + b3 + b4
1048 The unit element in the Hadamard EJA is inherited in the
1049 subalgebras generated by its elements::
1051 sage: J = HadamardEJA(5)
1053 b0 + b1 + b2 + b3 + b4
1054 sage: x = sum(J.gens())
1055 sage: A = x.subalgebra_generated_by(orthonormalize=False)
1058 sage: A.one().superalgebra_element()
1059 b0 + b1 + b2 + b3 + b4
1063 The identity element acts like the identity, regardless of
1064 whether or not we orthonormalize::
1066 sage: set_random_seed()
1067 sage: J = random_eja()
1068 sage: x = J.random_element()
1069 sage: J.one()*x == x and x*J.one() == x
1071 sage: A = x.subalgebra_generated_by()
1072 sage: y = A.random_element()
1073 sage: A.one()*y == y and y*A.one() == y
1078 sage: set_random_seed()
1079 sage: J = random_eja(field=QQ, orthonormalize=False)
1080 sage: x = J.random_element()
1081 sage: J.one()*x == x and x*J.one() == x
1083 sage: A = x.subalgebra_generated_by(orthonormalize=False)
1084 sage: y = A.random_element()
1085 sage: A.one()*y == y and y*A.one() == y
1088 The matrix of the unit element's operator is the identity,
1089 regardless of the base field and whether or not we
1092 sage: set_random_seed()
1093 sage: J = random_eja()
1094 sage: actual = J.one().operator().matrix()
1095 sage: expected = matrix.identity(J.base_ring(), J.dimension())
1096 sage: actual == expected
1098 sage: x = J.random_element()
1099 sage: A = x.subalgebra_generated_by()
1100 sage: actual = A.one().operator().matrix()
1101 sage: expected = matrix.identity(A.base_ring(), A.dimension())
1102 sage: actual == expected
1107 sage: set_random_seed()
1108 sage: J = random_eja(field=QQ, orthonormalize=False)
1109 sage: actual = J.one().operator().matrix()
1110 sage: expected = matrix.identity(J.base_ring(), J.dimension())
1111 sage: actual == expected
1113 sage: x = J.random_element()
1114 sage: A = x.subalgebra_generated_by(orthonormalize=False)
1115 sage: actual = A.one().operator().matrix()
1116 sage: expected = matrix.identity(A.base_ring(), A.dimension())
1117 sage: actual == expected
1120 Ensure that the cached unit element (often precomputed by
1121 hand) agrees with the computed one::
1123 sage: set_random_seed()
1124 sage: J = random_eja()
1125 sage: cached = J.one()
1126 sage: J.one.clear_cache()
1127 sage: J.one() == cached
1132 sage: set_random_seed()
1133 sage: J = random_eja(field=QQ, orthonormalize=False)
1134 sage: cached = J.one()
1135 sage: J.one.clear_cache()
1136 sage: J.one() == cached
1140 # We can brute-force compute the matrices of the operators
1141 # that correspond to the basis elements of this algebra.
1142 # If some linear combination of those basis elements is the
1143 # algebra identity, then the same linear combination of
1144 # their matrices has to be the identity matrix.
1146 # Of course, matrices aren't vectors in sage, so we have to
1147 # appeal to the "long vectors" isometry.
1148 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
1150 # Now we use basic linear algebra to find the coefficients,
1151 # of the matrices-as-vectors-linear-combination, which should
1152 # work for the original algebra basis too.
1153 A
= matrix(self
.base_ring(), oper_vecs
)
1155 # We used the isometry on the left-hand side already, but we
1156 # still need to do it for the right-hand side. Recall that we
1157 # wanted something that summed to the identity matrix.
1158 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
1160 # Now if there's an identity element in the algebra, this
1161 # should work. We solve on the left to avoid having to
1162 # transpose the matrix "A".
1163 return self
.from_vector(A
.solve_left(b
))
1166 def peirce_decomposition(self
, c
):
1168 The Peirce decomposition of this algebra relative to the
1171 In the future, this can be extended to a complete system of
1172 orthogonal idempotents.
1176 - ``c`` -- an idempotent of this algebra.
1180 A triple (J0, J5, J1) containing two subalgebras and one subspace
1183 - ``J0`` -- the algebra on the eigenspace of ``c.operator()``
1184 corresponding to the eigenvalue zero.
1186 - ``J5`` -- the eigenspace (NOT a subalgebra) of ``c.operator()``
1187 corresponding to the eigenvalue one-half.
1189 - ``J1`` -- the algebra on the eigenspace of ``c.operator()``
1190 corresponding to the eigenvalue one.
1192 These are the only possible eigenspaces for that operator, and this
1193 algebra is a direct sum of them. The spaces ``J0`` and ``J1`` are
1194 orthogonal, and are subalgebras of this algebra with the appropriate
1199 sage: from mjo.eja.eja_algebra import random_eja, RealSymmetricEJA
1203 The canonical example comes from the symmetric matrices, which
1204 decompose into diagonal and off-diagonal parts::
1206 sage: J = RealSymmetricEJA(3)
1207 sage: C = matrix(QQ, [ [1,0,0],
1211 sage: J0,J5,J1 = J.peirce_decomposition(c)
1213 Euclidean Jordan algebra of dimension 1...
1215 Vector space of degree 6 and dimension 2...
1217 Euclidean Jordan algebra of dimension 3...
1218 sage: J0.one().to_matrix()
1222 sage: orig_df = AA.options.display_format
1223 sage: AA.options.display_format = 'radical'
1224 sage: J.from_vector(J5.basis()[0]).to_matrix()
1228 sage: J.from_vector(J5.basis()[1]).to_matrix()
1232 sage: AA.options.display_format = orig_df
1233 sage: J1.one().to_matrix()
1240 Every algebra decomposes trivially with respect to its identity
1243 sage: set_random_seed()
1244 sage: J = random_eja()
1245 sage: J0,J5,J1 = J.peirce_decomposition(J.one())
1246 sage: J0.dimension() == 0 and J5.dimension() == 0
1248 sage: J1.superalgebra() == J and J1.dimension() == J.dimension()
1251 The decomposition is into eigenspaces, and its components are
1252 therefore necessarily orthogonal. Moreover, the identity
1253 elements in the two subalgebras are the projections onto their
1254 respective subspaces of the superalgebra's identity element::
1256 sage: set_random_seed()
1257 sage: J = random_eja()
1258 sage: x = J.random_element()
1259 sage: if not J.is_trivial():
1260 ....: while x.is_nilpotent():
1261 ....: x = J.random_element()
1262 sage: c = x.subalgebra_idempotent()
1263 sage: J0,J5,J1 = J.peirce_decomposition(c)
1265 sage: for (w,y,z) in zip(J0.basis(), J5.basis(), J1.basis()):
1266 ....: w = w.superalgebra_element()
1267 ....: y = J.from_vector(y)
1268 ....: z = z.superalgebra_element()
1269 ....: ipsum += w.inner_product(y).abs()
1270 ....: ipsum += w.inner_product(z).abs()
1271 ....: ipsum += y.inner_product(z).abs()
1274 sage: J1(c) == J1.one()
1276 sage: J0(J.one() - c) == J0.one()
1280 if not c
.is_idempotent():
1281 raise ValueError("element is not idempotent: %s" % c
)
1283 # Default these to what they should be if they turn out to be
1284 # trivial, because eigenspaces_left() won't return eigenvalues
1285 # corresponding to trivial spaces (e.g. it returns only the
1286 # eigenspace corresponding to lambda=1 if you take the
1287 # decomposition relative to the identity element).
1288 trivial
= self
.subalgebra(())
1289 J0
= trivial
# eigenvalue zero
1290 J5
= VectorSpace(self
.base_ring(), 0) # eigenvalue one-half
1291 J1
= trivial
# eigenvalue one
1293 for (eigval
, eigspace
) in c
.operator().matrix().right_eigenspaces():
1294 if eigval
== ~
(self
.base_ring()(2)):
1297 gens
= tuple( self
.from_vector(b
) for b
in eigspace
.basis() )
1298 subalg
= self
.subalgebra(gens
, check_axioms
=False)
1304 raise ValueError("unexpected eigenvalue: %s" % eigval
)
1309 def random_element(self
, thorough
=False):
1311 Return a random element of this algebra.
1313 Our algebra superclass method only returns a linear
1314 combination of at most two basis elements. We instead
1315 want the vector space "random element" method that
1316 returns a more diverse selection.
1320 - ``thorough`` -- (boolean; default False) whether or not we
1321 should generate irrational coefficients for the random
1322 element when our base ring is irrational; this slows the
1323 algebra operations to a crawl, but any truly random method
1327 # For a general base ring... maybe we can trust this to do the
1328 # right thing? Unlikely, but.
1329 V
= self
.vector_space()
1330 v
= V
.random_element()
1332 if self
.base_ring() is AA
:
1333 # The "random element" method of the algebraic reals is
1334 # stupid at the moment, and only returns integers between
1335 # -2 and 2, inclusive:
1337 # https://trac.sagemath.org/ticket/30875
1339 # Instead, we implement our own "random vector" method,
1340 # and then coerce that into the algebra. We use the vector
1341 # space degree here instead of the dimension because a
1342 # subalgebra could (for example) be spanned by only two
1343 # vectors, each with five coordinates. We need to
1344 # generate all five coordinates.
1346 v
*= QQbar
.random_element().real()
1348 v
*= QQ
.random_element()
1350 return self
.from_vector(V
.coordinate_vector(v
))
1352 def random_elements(self
, count
, thorough
=False):
1354 Return ``count`` random elements as a tuple.
1358 - ``thorough`` -- (boolean; default False) whether or not we
1359 should generate irrational coefficients for the random
1360 elements when our base ring is irrational; this slows the
1361 algebra operations to a crawl, but any truly random method
1366 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1370 sage: J = JordanSpinEJA(3)
1371 sage: x,y,z = J.random_elements(3)
1372 sage: all( [ x in J, y in J, z in J ])
1374 sage: len( J.random_elements(10) ) == 10
1378 return tuple( self
.random_element(thorough
)
1379 for idx
in range(count
) )
1383 def _charpoly_coefficients(self
):
1385 The `r` polynomial coefficients of the "characteristic polynomial
1390 sage: from mjo.eja.eja_algebra import random_eja
1394 The theory shows that these are all homogeneous polynomials of
1397 sage: set_random_seed()
1398 sage: J = random_eja()
1399 sage: all(p.is_homogeneous() for p in J._charpoly_coefficients())
1403 n
= self
.dimension()
1404 R
= self
.coordinate_polynomial_ring()
1406 F
= R
.fraction_field()
1409 # From a result in my book, these are the entries of the
1410 # basis representation of L_x.
1411 return sum( vars[k
]*self
.monomial(k
).operator().matrix()[i
,j
]
1414 L_x
= matrix(F
, n
, n
, L_x_i_j
)
1417 if self
.rank
.is_in_cache():
1419 # There's no need to pad the system with redundant
1420 # columns if we *know* they'll be redundant.
1423 # Compute an extra power in case the rank is equal to
1424 # the dimension (otherwise, we would stop at x^(r-1)).
1425 x_powers
= [ (L_x
**k
)*self
.one().to_vector()
1426 for k
in range(n
+1) ]
1427 A
= matrix
.column(F
, x_powers
[:n
])
1428 AE
= A
.extended_echelon_form()
1435 # The theory says that only the first "r" coefficients are
1436 # nonzero, and they actually live in the original polynomial
1437 # ring and not the fraction field. We negate them because in
1438 # the actual characteristic polynomial, they get moved to the
1439 # other side where x^r lives. We don't bother to trim A_rref
1440 # down to a square matrix and solve the resulting system,
1441 # because the upper-left r-by-r portion of A_rref is
1442 # guaranteed to be the identity matrix, so e.g.
1444 # A_rref.solve_right(Y)
1446 # would just be returning Y.
1447 return (-E
*b
)[:r
].change_ring(R
)
1452 Return the rank of this EJA.
1454 This is a cached method because we know the rank a priori for
1455 all of the algebras we can construct. Thus we can avoid the
1456 expensive ``_charpoly_coefficients()`` call unless we truly
1457 need to compute the whole characteristic polynomial.
1461 sage: from mjo.eja.eja_algebra import (HadamardEJA,
1462 ....: JordanSpinEJA,
1463 ....: RealSymmetricEJA,
1464 ....: ComplexHermitianEJA,
1465 ....: QuaternionHermitianEJA,
1470 The rank of the Jordan spin algebra is always two::
1472 sage: JordanSpinEJA(2).rank()
1474 sage: JordanSpinEJA(3).rank()
1476 sage: JordanSpinEJA(4).rank()
1479 The rank of the `n`-by-`n` Hermitian real, complex, or
1480 quaternion matrices is `n`::
1482 sage: RealSymmetricEJA(4).rank()
1484 sage: ComplexHermitianEJA(3).rank()
1486 sage: QuaternionHermitianEJA(2).rank()
1491 Ensure that every EJA that we know how to construct has a
1492 positive integer rank, unless the algebra is trivial in
1493 which case its rank will be zero::
1495 sage: set_random_seed()
1496 sage: J = random_eja()
1500 sage: r > 0 or (r == 0 and J.is_trivial())
1503 Ensure that computing the rank actually works, since the ranks
1504 of all simple algebras are known and will be cached by default::
1506 sage: set_random_seed() # long time
1507 sage: J = random_eja() # long time
1508 sage: cached = J.rank() # long time
1509 sage: J.rank.clear_cache() # long time
1510 sage: J.rank() == cached # long time
1514 return len(self
._charpoly
_coefficients
())
1517 def subalgebra(self
, basis
, **kwargs
):
1519 Create a subalgebra of this algebra from the given basis.
1521 from mjo
.eja
.eja_subalgebra
import FiniteDimensionalEJASubalgebra
1522 return FiniteDimensionalEJASubalgebra(self
, basis
, **kwargs
)
1525 def vector_space(self
):
1527 Return the vector space that underlies this algebra.
1531 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1535 sage: J = RealSymmetricEJA(2)
1536 sage: J.vector_space()
1537 Vector space of dimension 3 over...
1540 return self
.zero().to_vector().parent().ambient_vector_space()
1544 class RationalBasisEJA(FiniteDimensionalEJA
):
1546 New class for algebras whose supplied basis elements have all rational entries.
1550 sage: from mjo.eja.eja_algebra import BilinearFormEJA
1554 The supplied basis is orthonormalized by default::
1556 sage: B = matrix(QQ, [[1, 0, 0], [0, 25, -32], [0, -32, 41]])
1557 sage: J = BilinearFormEJA(B)
1558 sage: J.matrix_basis()
1575 # Abuse the check_field parameter to check that the entries of
1576 # out basis (in ambient coordinates) are in the field QQ.
1577 if not all( all(b_i
in QQ
for b_i
in b
.list()) for b
in basis
):
1578 raise TypeError("basis not rational")
1580 super().__init
__(basis
,
1584 check_field
=check_field
,
1587 self
._rational
_algebra
= None
1589 # There's no point in constructing the extra algebra if this
1590 # one is already rational.
1592 # Note: the same Jordan and inner-products work here,
1593 # because they are necessarily defined with respect to
1594 # ambient coordinates and not any particular basis.
1595 self
._rational
_algebra
= FiniteDimensionalEJA(
1600 associative
=self
.is_associative(),
1601 orthonormalize
=False,
1606 def _charpoly_coefficients(self
):
1610 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
1611 ....: JordanSpinEJA)
1615 The base ring of the resulting polynomial coefficients is what
1616 it should be, and not the rationals (unless the algebra was
1617 already over the rationals)::
1619 sage: J = JordanSpinEJA(3)
1620 sage: J._charpoly_coefficients()
1621 (X1^2 - X2^2 - X3^2, -2*X1)
1622 sage: a0 = J._charpoly_coefficients()[0]
1624 Algebraic Real Field
1625 sage: a0.base_ring()
1626 Algebraic Real Field
1629 if self
._rational
_algebra
is None:
1630 # There's no need to construct *another* algebra over the
1631 # rationals if this one is already over the
1632 # rationals. Likewise, if we never orthonormalized our
1633 # basis, we might as well just use the given one.
1634 return super()._charpoly
_coefficients
()
1636 # Do the computation over the rationals. The answer will be
1637 # the same, because all we've done is a change of basis.
1638 # Then, change back from QQ to our real base ring
1639 a
= ( a_i
.change_ring(self
.base_ring())
1640 for a_i
in self
._rational
_algebra
._charpoly
_coefficients
() )
1642 if self
._deortho
_matrix
is None:
1643 # This can happen if our base ring was, say, AA and we
1644 # chose not to (or didn't need to) orthonormalize. It's
1645 # still faster to do the computations over QQ even if
1646 # the numbers in the boxes stay the same.
1649 # Otherwise, convert the coordinate variables back to the
1650 # deorthonormalized ones.
1651 R
= self
.coordinate_polynomial_ring()
1652 from sage
.modules
.free_module_element
import vector
1653 X
= vector(R
, R
.gens())
1654 BX
= self
._deortho
_matrix
*X
1656 subs_dict
= { X[i]: BX[i] for i in range(len(X)) }
1657 return tuple( a_i
.subs(subs_dict
) for a_i
in a
)
1659 class ConcreteEJA(FiniteDimensionalEJA
):
1661 A class for the Euclidean Jordan algebras that we know by name.
1663 These are the Jordan algebras whose basis, multiplication table,
1664 rank, and so on are known a priori. More to the point, they are
1665 the Euclidean Jordan algebras for which we are able to conjure up
1666 a "random instance."
1670 sage: from mjo.eja.eja_algebra import ConcreteEJA
1674 Our basis is normalized with respect to the algebra's inner
1675 product, unless we specify otherwise::
1677 sage: set_random_seed()
1678 sage: J = ConcreteEJA.random_instance()
1679 sage: all( b.norm() == 1 for b in J.gens() )
1682 Since our basis is orthonormal with respect to the algebra's inner
1683 product, and since we know that this algebra is an EJA, any
1684 left-multiplication operator's matrix will be symmetric because
1685 natural->EJA basis representation is an isometry and within the
1686 EJA the operator is self-adjoint by the Jordan axiom::
1688 sage: set_random_seed()
1689 sage: J = ConcreteEJA.random_instance()
1690 sage: x = J.random_element()
1691 sage: x.operator().is_self_adjoint()
1696 def _max_random_instance_size():
1698 Return an integer "size" that is an upper bound on the size of
1699 this algebra when it is used in a random test
1700 case. Unfortunately, the term "size" is ambiguous -- when
1701 dealing with `R^n` under either the Hadamard or Jordan spin
1702 product, the "size" refers to the dimension `n`. When dealing
1703 with a matrix algebra (real symmetric or complex/quaternion
1704 Hermitian), it refers to the size of the matrix, which is far
1705 less than the dimension of the underlying vector space.
1707 This method must be implemented in each subclass.
1709 raise NotImplementedError
1712 def random_instance(cls
, *args
, **kwargs
):
1714 Return a random instance of this type of algebra.
1716 This method should be implemented in each subclass.
1718 from sage
.misc
.prandom
import choice
1719 eja_class
= choice(cls
.__subclasses
__())
1721 # These all bubble up to the RationalBasisEJA superclass
1722 # constructor, so any (kw)args valid there are also valid
1724 return eja_class
.random_instance(*args
, **kwargs
)
1729 def jordan_product(X
,Y
):
1730 return (X
*Y
+ Y
*X
)/2
1733 def trace_inner_product(X
,Y
):
1735 A trace inner-product for matrices that aren't embedded in the
1736 reals. It takes MATRICES as arguments, not EJA elements.
1738 return (X
*Y
).trace().real()
1740 class RealEmbeddedMatrixEJA(MatrixEJA
):
1742 def dimension_over_reals():
1744 The dimension of this matrix's base ring over the reals.
1746 The reals are dimension one over themselves, obviously; that's
1747 just `\mathbb{R}^{1}`. Likewise, the complex numbers `a + bi`
1748 have dimension two. Finally, the quaternions have dimension
1749 four over the reals.
1751 This is used to determine the size of the matrix returned from
1752 :meth:`real_embed`, among other things.
1754 raise NotImplementedError
1757 def real_embed(cls
,M
):
1759 Embed the matrix ``M`` into a space of real matrices.
1761 The matrix ``M`` can have entries in any field at the moment:
1762 the real numbers, complex numbers, or quaternions. And although
1763 they are not a field, we can probably support octonions at some
1764 point, too. This function returns a real matrix that "acts like"
1765 the original with respect to matrix multiplication; i.e.
1767 real_embed(M*N) = real_embed(M)*real_embed(N)
1770 if M
.ncols() != M
.nrows():
1771 raise ValueError("the matrix 'M' must be square")
1776 def real_unembed(cls
,M
):
1778 The inverse of :meth:`real_embed`.
1780 if M
.ncols() != M
.nrows():
1781 raise ValueError("the matrix 'M' must be square")
1782 if not ZZ(M
.nrows()).mod(cls
.dimension_over_reals()).is_zero():
1783 raise ValueError("the matrix 'M' must be a real embedding")
1788 def trace_inner_product(cls
,X
,Y
):
1790 Compute the trace inner-product of two real-embeddings.
1794 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
1795 ....: QuaternionHermitianEJA)
1799 sage: set_random_seed()
1800 sage: J = ComplexHermitianEJA.random_instance()
1801 sage: x,y = J.random_elements(2)
1802 sage: Xe = x.to_matrix()
1803 sage: Ye = y.to_matrix()
1804 sage: X = J.real_unembed(Xe)
1805 sage: Y = J.real_unembed(Ye)
1806 sage: expected = (X*Y).trace().real()
1807 sage: actual = J.trace_inner_product(Xe,Ye)
1808 sage: actual == expected
1813 sage: set_random_seed()
1814 sage: J = QuaternionHermitianEJA.random_instance()
1815 sage: x,y = J.random_elements(2)
1816 sage: Xe = x.to_matrix()
1817 sage: Ye = y.to_matrix()
1818 sage: X = J.real_unembed(Xe)
1819 sage: Y = J.real_unembed(Ye)
1820 sage: expected = (X*Y).trace().coefficient_tuple()[0]
1821 sage: actual = J.trace_inner_product(Xe,Ye)
1822 sage: actual == expected
1826 # This does in fact compute the real part of the trace.
1827 # If we compute the trace of e.g. a complex matrix M,
1828 # then we do so by adding up its diagonal entries --
1829 # call them z_1 through z_n. The real embedding of z_1
1830 # will be a 2-by-2 REAL matrix [a, b; -b, a] whose trace
1831 # as a REAL matrix will be 2*a = 2*Re(z_1). And so forth.
1832 return (X
*Y
).trace()/cls
.dimension_over_reals()
1834 class RealSymmetricEJA(ConcreteEJA
, RationalBasisEJA
, MatrixEJA
):
1836 The rank-n simple EJA consisting of real symmetric n-by-n
1837 matrices, the usual symmetric Jordan product, and the trace inner
1838 product. It has dimension `(n^2 + n)/2` over the reals.
1842 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1846 sage: J = RealSymmetricEJA(2)
1847 sage: b0, b1, b2 = J.gens()
1855 In theory, our "field" can be any subfield of the reals::
1857 sage: RealSymmetricEJA(2, field=RDF, check_axioms=True)
1858 Euclidean Jordan algebra of dimension 3 over Real Double Field
1859 sage: RealSymmetricEJA(2, field=RR, check_axioms=True)
1860 Euclidean Jordan algebra of dimension 3 over Real Field with
1861 53 bits of precision
1865 The dimension of this algebra is `(n^2 + n) / 2`::
1867 sage: set_random_seed()
1868 sage: n_max = RealSymmetricEJA._max_random_instance_size()
1869 sage: n = ZZ.random_element(1, n_max)
1870 sage: J = RealSymmetricEJA(n)
1871 sage: J.dimension() == (n^2 + n)/2
1874 The Jordan multiplication is what we think it is::
1876 sage: set_random_seed()
1877 sage: J = RealSymmetricEJA.random_instance()
1878 sage: x,y = J.random_elements(2)
1879 sage: actual = (x*y).to_matrix()
1880 sage: X = x.to_matrix()
1881 sage: Y = y.to_matrix()
1882 sage: expected = (X*Y + Y*X)/2
1883 sage: actual == expected
1885 sage: J(expected) == x*y
1888 We can change the generator prefix::
1890 sage: RealSymmetricEJA(3, prefix='q').gens()
1891 (q0, q1, q2, q3, q4, q5)
1893 We can construct the (trivial) algebra of rank zero::
1895 sage: RealSymmetricEJA(0)
1896 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1900 def _denormalized_basis(cls
, n
, field
):
1902 Return a basis for the space of real symmetric n-by-n matrices.
1906 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1910 sage: set_random_seed()
1911 sage: n = ZZ.random_element(1,5)
1912 sage: B = RealSymmetricEJA._denormalized_basis(n,ZZ)
1913 sage: all( M.is_symmetric() for M in B)
1917 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1921 for j
in range(i
+1):
1922 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
1926 Sij
= Eij
+ Eij
.transpose()
1932 def _max_random_instance_size():
1933 return 4 # Dimension 10
1936 def random_instance(cls
, **kwargs
):
1938 Return a random instance of this type of algebra.
1940 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1941 return cls(n
, **kwargs
)
1943 def __init__(self
, n
, field
=AA
, **kwargs
):
1944 # We know this is a valid EJA, but will double-check
1945 # if the user passes check_axioms=True.
1946 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
1952 super().__init
__(self
._denormalized
_basis
(n
,field
),
1953 self
.jordan_product
,
1954 self
.trace_inner_product
,
1956 associative
=associative
,
1959 # TODO: this could be factored out somehow, but is left here
1960 # because the MatrixEJA is not presently a subclass of the
1961 # FDEJA class that defines rank() and one().
1962 self
.rank
.set_cache(n
)
1963 idV
= self
.matrix_space().one()
1964 self
.one
.set_cache(self(idV
))
1968 class ComplexMatrixEJA(RealEmbeddedMatrixEJA
):
1969 # A manual dictionary-cache for the complex_extension() method,
1970 # since apparently @classmethods can't also be @cached_methods.
1971 _complex_extension
= {}
1974 def complex_extension(cls
,field
):
1976 The complex field that we embed/unembed, as an extension
1977 of the given ``field``.
1979 if field
in cls
._complex
_extension
:
1980 return cls
._complex
_extension
[field
]
1982 # Sage doesn't know how to adjoin the complex "i" (the root of
1983 # x^2 + 1) to a field in a general way. Here, we just enumerate
1984 # all of the cases that I have cared to support so far.
1986 # Sage doesn't know how to embed AA into QQbar, i.e. how
1987 # to adjoin sqrt(-1) to AA.
1989 elif not field
.is_exact():
1991 F
= field
.complex_field()
1993 # Works for QQ and... maybe some other fields.
1994 R
= PolynomialRing(field
, 'z')
1996 F
= field
.extension(z
**2 + 1, 'I', embedding
=CLF(-1).sqrt())
1998 cls
._complex
_extension
[field
] = F
2002 def dimension_over_reals():
2006 def real_embed(cls
,M
):
2008 Embed the n-by-n complex matrix ``M`` into the space of real
2009 matrices of size 2n-by-2n via the map the sends each entry `z = a +
2010 bi` to the block matrix ``[[a,b],[-b,a]]``.
2014 sage: from mjo.eja.eja_algebra import ComplexMatrixEJA
2018 sage: F = QuadraticField(-1, 'I')
2019 sage: x1 = F(4 - 2*i)
2020 sage: x2 = F(1 + 2*i)
2023 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
2024 sage: ComplexMatrixEJA.real_embed(M)
2033 Embedding is a homomorphism (isomorphism, in fact)::
2035 sage: set_random_seed()
2036 sage: n = ZZ.random_element(3)
2037 sage: F = QuadraticField(-1, 'I')
2038 sage: X = random_matrix(F, n)
2039 sage: Y = random_matrix(F, n)
2040 sage: Xe = ComplexMatrixEJA.real_embed(X)
2041 sage: Ye = ComplexMatrixEJA.real_embed(Y)
2042 sage: XYe = ComplexMatrixEJA.real_embed(X*Y)
2047 super().real_embed(M
)
2050 # We don't need any adjoined elements...
2051 field
= M
.base_ring().base_ring()
2057 blocks
.append(matrix(field
, 2, [ [ a
, b
],
2060 return matrix
.block(field
, n
, blocks
)
2064 def real_unembed(cls
,M
):
2066 The inverse of _embed_complex_matrix().
2070 sage: from mjo.eja.eja_algebra import ComplexMatrixEJA
2074 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
2075 ....: [-2, 1, -4, 3],
2076 ....: [ 9, 10, 11, 12],
2077 ....: [-10, 9, -12, 11] ])
2078 sage: ComplexMatrixEJA.real_unembed(A)
2080 [ 10*I + 9 12*I + 11]
2084 Unembedding is the inverse of embedding::
2086 sage: set_random_seed()
2087 sage: F = QuadraticField(-1, 'I')
2088 sage: M = random_matrix(F, 3)
2089 sage: Me = ComplexMatrixEJA.real_embed(M)
2090 sage: ComplexMatrixEJA.real_unembed(Me) == M
2094 super().real_unembed(M
)
2096 d
= cls
.dimension_over_reals()
2097 F
= cls
.complex_extension(M
.base_ring())
2100 # Go top-left to bottom-right (reading order), converting every
2101 # 2-by-2 block we see to a single complex element.
2103 for k
in range(n
/d
):
2104 for j
in range(n
/d
):
2105 submat
= M
[d
*k
:d
*k
+d
,d
*j
:d
*j
+d
]
2106 if submat
[0,0] != submat
[1,1]:
2107 raise ValueError('bad on-diagonal submatrix')
2108 if submat
[0,1] != -submat
[1,0]:
2109 raise ValueError('bad off-diagonal submatrix')
2110 z
= submat
[0,0] + submat
[0,1]*i
2113 return matrix(F
, n
/d
, elements
)
2116 class ComplexHermitianEJA(ConcreteEJA
, RationalBasisEJA
, ComplexMatrixEJA
):
2118 The rank-n simple EJA consisting of complex Hermitian n-by-n
2119 matrices over the real numbers, the usual symmetric Jordan product,
2120 and the real-part-of-trace inner product. It has dimension `n^2` over
2125 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
2129 In theory, our "field" can be any subfield of the reals::
2131 sage: ComplexHermitianEJA(2, field=RDF, check_axioms=True)
2132 Euclidean Jordan algebra of dimension 4 over Real Double Field
2133 sage: ComplexHermitianEJA(2, field=RR, check_axioms=True)
2134 Euclidean Jordan algebra of dimension 4 over Real Field with
2135 53 bits of precision
2139 The dimension of this algebra is `n^2`::
2141 sage: set_random_seed()
2142 sage: n_max = ComplexHermitianEJA._max_random_instance_size()
2143 sage: n = ZZ.random_element(1, n_max)
2144 sage: J = ComplexHermitianEJA(n)
2145 sage: J.dimension() == n^2
2148 The Jordan multiplication is what we think it is::
2150 sage: set_random_seed()
2151 sage: J = ComplexHermitianEJA.random_instance()
2152 sage: x,y = J.random_elements(2)
2153 sage: actual = (x*y).to_matrix()
2154 sage: X = x.to_matrix()
2155 sage: Y = y.to_matrix()
2156 sage: expected = (X*Y + Y*X)/2
2157 sage: actual == expected
2159 sage: J(expected) == x*y
2162 We can change the generator prefix::
2164 sage: ComplexHermitianEJA(2, prefix='z').gens()
2167 We can construct the (trivial) algebra of rank zero::
2169 sage: ComplexHermitianEJA(0)
2170 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2175 def _denormalized_basis(cls
, n
, field
):
2177 Returns a basis for the space of complex Hermitian n-by-n matrices.
2179 Why do we embed these? Basically, because all of numerical linear
2180 algebra assumes that you're working with vectors consisting of `n`
2181 entries from a field and scalars from the same field. There's no way
2182 to tell SageMath that (for example) the vectors contain complex
2183 numbers, while the scalar field is real.
2187 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
2191 sage: set_random_seed()
2192 sage: n = ZZ.random_element(1,5)
2193 sage: B = ComplexHermitianEJA._denormalized_basis(n,ZZ)
2194 sage: all( M.is_symmetric() for M in B)
2198 R
= PolynomialRing(ZZ
, 'z')
2200 F
= ZZ
.extension(z
**2 + 1, 'I')
2203 # This is like the symmetric case, but we need to be careful:
2205 # * We want conjugate-symmetry, not just symmetry.
2206 # * The diagonal will (as a result) be real.
2209 Eij
= matrix
.zero(F
,n
)
2211 for j
in range(i
+1):
2215 Sij
= cls
.real_embed(Eij
)
2218 # The second one has a minus because it's conjugated.
2219 Eij
[j
,i
] = 1 # Eij = Eij + Eij.transpose()
2220 Sij_real
= cls
.real_embed(Eij
)
2222 # Eij = I*Eij - I*Eij.transpose()
2225 Sij_imag
= cls
.real_embed(Eij
)
2231 # Since we embedded the entries, we can drop back to the
2232 # desired real "field" instead of the extension "F".
2233 return tuple( s
.change_ring(field
) for s
in S
)
2236 def __init__(self
, n
, field
=AA
, **kwargs
):
2237 # We know this is a valid EJA, but will double-check
2238 # if the user passes check_axioms=True.
2239 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2245 super().__init
__(self
._denormalized
_basis
(n
,field
),
2246 self
.jordan_product
,
2247 self
.trace_inner_product
,
2249 associative
=associative
,
2251 # TODO: this could be factored out somehow, but is left here
2252 # because the MatrixEJA is not presently a subclass of the
2253 # FDEJA class that defines rank() and one().
2254 self
.rank
.set_cache(n
)
2255 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
2256 self
.one
.set_cache(self(idV
))
2259 def _max_random_instance_size():
2260 return 3 # Dimension 9
2263 def random_instance(cls
, **kwargs
):
2265 Return a random instance of this type of algebra.
2267 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2268 return cls(n
, **kwargs
)
2270 class QuaternionMatrixEJA(RealEmbeddedMatrixEJA
):
2272 # A manual dictionary-cache for the quaternion_extension() method,
2273 # since apparently @classmethods can't also be @cached_methods.
2274 _quaternion_extension
= {}
2277 def quaternion_extension(cls
,field
):
2279 The quaternion field that we embed/unembed, as an extension
2280 of the given ``field``.
2282 if field
in cls
._quaternion
_extension
:
2283 return cls
._quaternion
_extension
[field
]
2285 Q
= QuaternionAlgebra(field
,-1,-1)
2287 cls
._quaternion
_extension
[field
] = Q
2291 def dimension_over_reals():
2295 def real_embed(cls
,M
):
2297 Embed the n-by-n quaternion matrix ``M`` into the space of real
2298 matrices of size 4n-by-4n by first sending each quaternion entry `z
2299 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
2300 c+di],[-c + di, a-bi]]`, and then embedding those into a real
2305 sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA
2309 sage: Q = QuaternionAlgebra(QQ,-1,-1)
2310 sage: i,j,k = Q.gens()
2311 sage: x = 1 + 2*i + 3*j + 4*k
2312 sage: M = matrix(Q, 1, [[x]])
2313 sage: QuaternionMatrixEJA.real_embed(M)
2319 Embedding is a homomorphism (isomorphism, in fact)::
2321 sage: set_random_seed()
2322 sage: n = ZZ.random_element(2)
2323 sage: Q = QuaternionAlgebra(QQ,-1,-1)
2324 sage: X = random_matrix(Q, n)
2325 sage: Y = random_matrix(Q, n)
2326 sage: Xe = QuaternionMatrixEJA.real_embed(X)
2327 sage: Ye = QuaternionMatrixEJA.real_embed(Y)
2328 sage: XYe = QuaternionMatrixEJA.real_embed(X*Y)
2333 super().real_embed(M
)
2334 quaternions
= M
.base_ring()
2337 F
= QuadraticField(-1, 'I')
2342 t
= z
.coefficient_tuple()
2347 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
2348 [-c
+ d
*i
, a
- b
*i
]])
2349 realM
= ComplexMatrixEJA
.real_embed(cplxM
)
2350 blocks
.append(realM
)
2352 # We should have real entries by now, so use the realest field
2353 # we've got for the return value.
2354 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
2359 def real_unembed(cls
,M
):
2361 The inverse of _embed_quaternion_matrix().
2365 sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA
2369 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
2370 ....: [-2, 1, -4, 3],
2371 ....: [-3, 4, 1, -2],
2372 ....: [-4, -3, 2, 1]])
2373 sage: QuaternionMatrixEJA.real_unembed(M)
2374 [1 + 2*i + 3*j + 4*k]
2378 Unembedding is the inverse of embedding::
2380 sage: set_random_seed()
2381 sage: Q = QuaternionAlgebra(QQ, -1, -1)
2382 sage: M = random_matrix(Q, 3)
2383 sage: Me = QuaternionMatrixEJA.real_embed(M)
2384 sage: QuaternionMatrixEJA.real_unembed(Me) == M
2388 super().real_unembed(M
)
2390 d
= cls
.dimension_over_reals()
2392 # Use the base ring of the matrix to ensure that its entries can be
2393 # multiplied by elements of the quaternion algebra.
2394 Q
= cls
.quaternion_extension(M
.base_ring())
2397 # Go top-left to bottom-right (reading order), converting every
2398 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
2401 for l
in range(n
/d
):
2402 for m
in range(n
/d
):
2403 submat
= ComplexMatrixEJA
.real_unembed(
2404 M
[d
*l
:d
*l
+d
,d
*m
:d
*m
+d
] )
2405 if submat
[0,0] != submat
[1,1].conjugate():
2406 raise ValueError('bad on-diagonal submatrix')
2407 if submat
[0,1] != -submat
[1,0].conjugate():
2408 raise ValueError('bad off-diagonal submatrix')
2409 z
= submat
[0,0].real()
2410 z
+= submat
[0,0].imag()*i
2411 z
+= submat
[0,1].real()*j
2412 z
+= submat
[0,1].imag()*k
2415 return matrix(Q
, n
/d
, elements
)
2418 class QuaternionHermitianEJA(ConcreteEJA
,
2420 QuaternionMatrixEJA
):
2422 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
2423 matrices, the usual symmetric Jordan product, and the
2424 real-part-of-trace inner product. It has dimension `2n^2 - n` over
2429 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2433 In theory, our "field" can be any subfield of the reals::
2435 sage: QuaternionHermitianEJA(2, field=RDF, check_axioms=True)
2436 Euclidean Jordan algebra of dimension 6 over Real Double Field
2437 sage: QuaternionHermitianEJA(2, field=RR, check_axioms=True)
2438 Euclidean Jordan algebra of dimension 6 over Real Field with
2439 53 bits of precision
2443 The dimension of this algebra is `2*n^2 - n`::
2445 sage: set_random_seed()
2446 sage: n_max = QuaternionHermitianEJA._max_random_instance_size()
2447 sage: n = ZZ.random_element(1, n_max)
2448 sage: J = QuaternionHermitianEJA(n)
2449 sage: J.dimension() == 2*(n^2) - n
2452 The Jordan multiplication is what we think it is::
2454 sage: set_random_seed()
2455 sage: J = QuaternionHermitianEJA.random_instance()
2456 sage: x,y = J.random_elements(2)
2457 sage: actual = (x*y).to_matrix()
2458 sage: X = x.to_matrix()
2459 sage: Y = y.to_matrix()
2460 sage: expected = (X*Y + Y*X)/2
2461 sage: actual == expected
2463 sage: J(expected) == x*y
2466 We can change the generator prefix::
2468 sage: QuaternionHermitianEJA(2, prefix='a').gens()
2469 (a0, a1, a2, a3, a4, a5)
2471 We can construct the (trivial) algebra of rank zero::
2473 sage: QuaternionHermitianEJA(0)
2474 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2478 def _denormalized_basis(cls
, n
, field
):
2480 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
2482 Why do we embed these? Basically, because all of numerical
2483 linear algebra assumes that you're working with vectors consisting
2484 of `n` entries from a field and scalars from the same field. There's
2485 no way to tell SageMath that (for example) the vectors contain
2486 complex numbers, while the scalar field is real.
2490 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2494 sage: set_random_seed()
2495 sage: n = ZZ.random_element(1,5)
2496 sage: B = QuaternionHermitianEJA._denormalized_basis(n,ZZ)
2497 sage: all( M.is_symmetric() for M in B )
2501 Q
= QuaternionAlgebra(QQ
,-1,-1)
2504 # This is like the symmetric case, but we need to be careful:
2506 # * We want conjugate-symmetry, not just symmetry.
2507 # * The diagonal will (as a result) be real.
2510 Eij
= matrix
.zero(Q
,n
)
2512 for j
in range(i
+1):
2516 Sij
= cls
.real_embed(Eij
)
2519 # The second, third, and fourth ones have a minus
2520 # because they're conjugated.
2521 # Eij = Eij + Eij.transpose()
2523 Sij_real
= cls
.real_embed(Eij
)
2525 # Eij = I*(Eij - Eij.transpose())
2528 Sij_I
= cls
.real_embed(Eij
)
2530 # Eij = J*(Eij - Eij.transpose())
2533 Sij_J
= cls
.real_embed(Eij
)
2535 # Eij = K*(Eij - Eij.transpose())
2538 Sij_K
= cls
.real_embed(Eij
)
2544 # Since we embedded the entries, we can drop back to the
2545 # desired real "field" instead of the quaternion algebra "Q".
2546 return tuple( s
.change_ring(field
) for s
in S
)
2549 def __init__(self
, n
, field
=AA
, **kwargs
):
2550 # We know this is a valid EJA, but will double-check
2551 # if the user passes check_axioms=True.
2552 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2558 super().__init
__(self
._denormalized
_basis
(n
,field
),
2559 self
.jordan_product
,
2560 self
.trace_inner_product
,
2562 associative
=associative
,
2565 # TODO: this could be factored out somehow, but is left here
2566 # because the MatrixEJA is not presently a subclass of the
2567 # FDEJA class that defines rank() and one().
2568 self
.rank
.set_cache(n
)
2569 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
2570 self
.one
.set_cache(self(idV
))
2574 def _max_random_instance_size():
2576 The maximum rank of a random QuaternionHermitianEJA.
2578 return 2 # Dimension 6
2581 def random_instance(cls
, **kwargs
):
2583 Return a random instance of this type of algebra.
2585 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2586 return cls(n
, **kwargs
)
2588 class OctonionHermitianEJA(FiniteDimensionalEJA
, MatrixEJA
):
2592 sage: from mjo.eja.eja_algebra import (FiniteDimensionalEJA,
2593 ....: OctonionHermitianEJA)
2597 The 3-by-3 algebra satisfies the axioms of an EJA::
2599 sage: OctonionHermitianEJA(3, # long time
2600 ....: field=QQ, # long time
2601 ....: orthonormalize=False, # long time
2602 ....: check_axioms=True) # long time
2603 Euclidean Jordan algebra of dimension 27 over Rational Field
2605 After a change-of-basis, the 2-by-2 algebra has the same
2606 multiplication table as the ten-dimensional Jordan spin algebra::
2608 sage: b = OctonionHermitianEJA._denormalized_basis(2,QQ)
2609 sage: basis = (b[0] + b[9],) + b[1:9] + (b[0] - b[9],)
2610 sage: jp = OctonionHermitianEJA.jordan_product
2611 sage: ip = OctonionHermitianEJA.trace_inner_product
2612 sage: J = FiniteDimensionalEJA(basis,
2616 ....: orthonormalize=False)
2617 sage: J.multiplication_table()
2618 +----++----+----+----+----+----+----+----+----+----+----+
2619 | * || b0 | b1 | b2 | b3 | b4 | b5 | b6 | b7 | b8 | b9 |
2620 +====++====+====+====+====+====+====+====+====+====+====+
2621 | b0 || b0 | b1 | b2 | b3 | b4 | b5 | b6 | b7 | b8 | b9 |
2622 +----++----+----+----+----+----+----+----+----+----+----+
2623 | b1 || b1 | b0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
2624 +----++----+----+----+----+----+----+----+----+----+----+
2625 | b2 || b2 | 0 | b0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
2626 +----++----+----+----+----+----+----+----+----+----+----+
2627 | b3 || b3 | 0 | 0 | b0 | 0 | 0 | 0 | 0 | 0 | 0 |
2628 +----++----+----+----+----+----+----+----+----+----+----+
2629 | b4 || b4 | 0 | 0 | 0 | b0 | 0 | 0 | 0 | 0 | 0 |
2630 +----++----+----+----+----+----+----+----+----+----+----+
2631 | b5 || b5 | 0 | 0 | 0 | 0 | b0 | 0 | 0 | 0 | 0 |
2632 +----++----+----+----+----+----+----+----+----+----+----+
2633 | b6 || b6 | 0 | 0 | 0 | 0 | 0 | b0 | 0 | 0 | 0 |
2634 +----++----+----+----+----+----+----+----+----+----+----+
2635 | b7 || b7 | 0 | 0 | 0 | 0 | 0 | 0 | b0 | 0 | 0 |
2636 +----++----+----+----+----+----+----+----+----+----+----+
2637 | b8 || b8 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | b0 | 0 |
2638 +----++----+----+----+----+----+----+----+----+----+----+
2639 | b9 || b9 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | b0 |
2640 +----++----+----+----+----+----+----+----+----+----+----+
2644 We can actually construct the 27-dimensional Albert algebra,
2645 and we get the right unit element if we recompute it::
2647 sage: J = OctonionHermitianEJA(3, # long time
2648 ....: field=QQ, # long time
2649 ....: orthonormalize=False) # long time
2650 sage: J.one.clear_cache() # long time
2651 sage: J.one() # long time
2653 sage: J.one().to_matrix() # long time
2662 The 2-by-2 algebra is isomorphic to the ten-dimensional Jordan
2663 spin algebra, but just to be sure, we recompute its rank::
2665 sage: J = OctonionHermitianEJA(2, # long time
2666 ....: field=QQ, # long time
2667 ....: orthonormalize=False) # long time
2668 sage: J.rank.clear_cache() # long time
2669 sage: J.rank() # long time
2672 def __init__(self
, n
, field
=AA
, **kwargs
):
2674 # Otherwise we don't get an EJA.
2675 raise ValueError("n cannot exceed 3")
2677 # We know this is a valid EJA, but will double-check
2678 # if the user passes check_axioms=True.
2679 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2681 super().__init
__(self
._denormalized
_basis
(n
,field
),
2682 self
.jordan_product
,
2683 self
.trace_inner_product
,
2687 # TODO: this could be factored out somehow, but is left here
2688 # because the MatrixEJA is not presently a subclass of the
2689 # FDEJA class that defines rank() and one().
2690 self
.rank
.set_cache(n
)
2691 idV
= self
.matrix_space().one()
2692 self
.one
.set_cache(self(idV
))
2696 def _denormalized_basis(cls
, n
, field
):
2698 Returns a basis for the space of octonion Hermitian n-by-n
2703 sage: from mjo.eja.eja_algebra import OctonionHermitianEJA
2707 sage: B = OctonionHermitianEJA._denormalized_basis(3,QQ)
2708 sage: all( M.is_hermitian() for M in B )
2714 from mjo
.octonions
import OctonionMatrixAlgebra
2715 MS
= OctonionMatrixAlgebra(n
, scalars
=field
)
2716 es
= MS
.entry_algebra().gens()
2720 for j
in range(i
+1):
2722 E_ii
= MS
.monomial( (i
,j
,es
[0]) )
2726 E_ij
= MS
.monomial( (i
,j
,e
) )
2728 # If the conjugate has a negative sign in front
2729 # of it, (j,i,ec) won't be a monomial!
2730 if (j
,i
,ec
) in MS
.indices():
2731 E_ij
+= MS
.monomial( (j
,i
,ec
) )
2733 E_ij
-= MS
.monomial( (j
,i
,-ec
) )
2736 return tuple( basis
)
2739 def trace_inner_product(X
,Y
):
2741 The octonions don't know that the reals are embedded in them,
2742 so we have to take the e0 component ourselves.
2746 sage: from mjo.eja.eja_algebra import OctonionHermitianEJA
2750 sage: J = OctonionHermitianEJA(2,field=QQ,orthonormalize=False)
2751 sage: I = J.one().to_matrix()
2752 sage: J.trace_inner_product(I, -I)
2756 return (X
*Y
).trace().real().coefficient(0)
2758 class HadamardEJA(ConcreteEJA
, RationalBasisEJA
):
2760 Return the Euclidean Jordan Algebra corresponding to the set
2761 `R^n` under the Hadamard product.
2763 Note: this is nothing more than the Cartesian product of ``n``
2764 copies of the spin algebra. Once Cartesian product algebras
2765 are implemented, this can go.
2769 sage: from mjo.eja.eja_algebra import HadamardEJA
2773 This multiplication table can be verified by hand::
2775 sage: J = HadamardEJA(3)
2776 sage: b0,b1,b2 = J.gens()
2792 We can change the generator prefix::
2794 sage: HadamardEJA(3, prefix='r').gens()
2798 def __init__(self
, n
, field
=AA
, **kwargs
):
2800 jordan_product
= lambda x
,y
: x
2801 inner_product
= lambda x
,y
: x
2803 def jordan_product(x
,y
):
2805 return P( xi
*yi
for (xi
,yi
) in zip(x
,y
) )
2807 def inner_product(x
,y
):
2810 # New defaults for keyword arguments. Don't orthonormalize
2811 # because our basis is already orthonormal with respect to our
2812 # inner-product. Don't check the axioms, because we know this
2813 # is a valid EJA... but do double-check if the user passes
2814 # check_axioms=True. Note: we DON'T override the "check_field"
2815 # default here, because the user can pass in a field!
2816 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2817 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2819 column_basis
= tuple( b
.column()
2820 for b
in FreeModule(field
, n
).basis() )
2821 super().__init
__(column_basis
,
2827 self
.rank
.set_cache(n
)
2830 self
.one
.set_cache( self
.zero() )
2832 self
.one
.set_cache( sum(self
.gens()) )
2835 def _max_random_instance_size():
2837 The maximum dimension of a random HadamardEJA.
2842 def random_instance(cls
, **kwargs
):
2844 Return a random instance of this type of algebra.
2846 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2847 return cls(n
, **kwargs
)
2850 class BilinearFormEJA(ConcreteEJA
, RationalBasisEJA
):
2852 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2853 with the half-trace inner product and jordan product ``x*y =
2854 (<Bx,y>,y_bar>, x0*y_bar + y0*x_bar)`` where `B = 1 \times B22` is
2855 a symmetric positive-definite "bilinear form" matrix. Its
2856 dimension is the size of `B`, and it has rank two in dimensions
2857 larger than two. It reduces to the ``JordanSpinEJA`` when `B` is
2858 the identity matrix of order ``n``.
2860 We insist that the one-by-one upper-left identity block of `B` be
2861 passed in as well so that we can be passed a matrix of size zero
2862 to construct a trivial algebra.
2866 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
2867 ....: JordanSpinEJA)
2871 When no bilinear form is specified, the identity matrix is used,
2872 and the resulting algebra is the Jordan spin algebra::
2874 sage: B = matrix.identity(AA,3)
2875 sage: J0 = BilinearFormEJA(B)
2876 sage: J1 = JordanSpinEJA(3)
2877 sage: J0.multiplication_table() == J0.multiplication_table()
2880 An error is raised if the matrix `B` does not correspond to a
2881 positive-definite bilinear form::
2883 sage: B = matrix.random(QQ,2,3)
2884 sage: J = BilinearFormEJA(B)
2885 Traceback (most recent call last):
2887 ValueError: bilinear form is not positive-definite
2888 sage: B = matrix.zero(QQ,3)
2889 sage: J = BilinearFormEJA(B)
2890 Traceback (most recent call last):
2892 ValueError: bilinear form is not positive-definite
2896 We can create a zero-dimensional algebra::
2898 sage: B = matrix.identity(AA,0)
2899 sage: J = BilinearFormEJA(B)
2903 We can check the multiplication condition given in the Jordan, von
2904 Neumann, and Wigner paper (and also discussed on my "On the
2905 symmetry..." paper). Note that this relies heavily on the standard
2906 choice of basis, as does anything utilizing the bilinear form
2907 matrix. We opt not to orthonormalize the basis, because if we
2908 did, we would have to normalize the `s_{i}` in a similar manner::
2910 sage: set_random_seed()
2911 sage: n = ZZ.random_element(5)
2912 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2913 sage: B11 = matrix.identity(QQ,1)
2914 sage: B22 = M.transpose()*M
2915 sage: B = block_matrix(2,2,[ [B11,0 ],
2917 sage: J = BilinearFormEJA(B, orthonormalize=False)
2918 sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
2919 sage: V = J.vector_space()
2920 sage: sis = [ J( V([0] + (M.inverse()*ei).list()).column() )
2921 ....: for ei in eis ]
2922 sage: actual = [ sis[i]*sis[j]
2923 ....: for i in range(n-1)
2924 ....: for j in range(n-1) ]
2925 sage: expected = [ J.one() if i == j else J.zero()
2926 ....: for i in range(n-1)
2927 ....: for j in range(n-1) ]
2928 sage: actual == expected
2932 def __init__(self
, B
, field
=AA
, **kwargs
):
2933 # The matrix "B" is supplied by the user in most cases,
2934 # so it makes sense to check whether or not its positive-
2935 # definite unless we are specifically asked not to...
2936 if ("check_axioms" not in kwargs
) or kwargs
["check_axioms"]:
2937 if not B
.is_positive_definite():
2938 raise ValueError("bilinear form is not positive-definite")
2940 # However, all of the other data for this EJA is computed
2941 # by us in manner that guarantees the axioms are
2942 # satisfied. So, again, unless we are specifically asked to
2943 # verify things, we'll skip the rest of the checks.
2944 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2946 def inner_product(x
,y
):
2947 return (y
.T
*B
*x
)[0,0]
2949 def jordan_product(x
,y
):
2955 z0
= inner_product(y
,x
)
2956 zbar
= y0
*xbar
+ x0
*ybar
2957 return P([z0
] + zbar
.list())
2960 column_basis
= tuple( b
.column()
2961 for b
in FreeModule(field
, n
).basis() )
2963 # TODO: I haven't actually checked this, but it seems legit.
2968 super().__init
__(column_basis
,
2972 associative
=associative
,
2975 # The rank of this algebra is two, unless we're in a
2976 # one-dimensional ambient space (because the rank is bounded
2977 # by the ambient dimension).
2978 self
.rank
.set_cache(min(n
,2))
2981 self
.one
.set_cache( self
.zero() )
2983 self
.one
.set_cache( self
.monomial(0) )
2986 def _max_random_instance_size():
2988 The maximum dimension of a random BilinearFormEJA.
2993 def random_instance(cls
, **kwargs
):
2995 Return a random instance of this algebra.
2997 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2999 B
= matrix
.identity(ZZ
, n
)
3000 return cls(B
, **kwargs
)
3002 B11
= matrix
.identity(ZZ
, 1)
3003 M
= matrix
.random(ZZ
, n
-1)
3004 I
= matrix
.identity(ZZ
, n
-1)
3006 while alpha
.is_zero():
3007 alpha
= ZZ
.random_element().abs()
3008 B22
= M
.transpose()*M
+ alpha
*I
3010 from sage
.matrix
.special
import block_matrix
3011 B
= block_matrix(2,2, [ [B11
, ZZ(0) ],
3014 return cls(B
, **kwargs
)
3017 class JordanSpinEJA(BilinearFormEJA
):
3019 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
3020 with the usual inner product and jordan product ``x*y =
3021 (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
3026 sage: from mjo.eja.eja_algebra import JordanSpinEJA
3030 This multiplication table can be verified by hand::
3032 sage: J = JordanSpinEJA(4)
3033 sage: b0,b1,b2,b3 = J.gens()
3049 We can change the generator prefix::
3051 sage: JordanSpinEJA(2, prefix='B').gens()
3056 Ensure that we have the usual inner product on `R^n`::
3058 sage: set_random_seed()
3059 sage: J = JordanSpinEJA.random_instance()
3060 sage: x,y = J.random_elements(2)
3061 sage: actual = x.inner_product(y)
3062 sage: expected = x.to_vector().inner_product(y.to_vector())
3063 sage: actual == expected
3067 def __init__(self
, n
, *args
, **kwargs
):
3068 # This is a special case of the BilinearFormEJA with the
3069 # identity matrix as its bilinear form.
3070 B
= matrix
.identity(ZZ
, n
)
3072 # Don't orthonormalize because our basis is already
3073 # orthonormal with respect to our inner-product.
3074 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
3076 # But also don't pass check_field=False here, because the user
3077 # can pass in a field!
3078 super().__init
__(B
, *args
, **kwargs
)
3081 def _max_random_instance_size():
3083 The maximum dimension of a random JordanSpinEJA.
3088 def random_instance(cls
, **kwargs
):
3090 Return a random instance of this type of algebra.
3092 Needed here to override the implementation for ``BilinearFormEJA``.
3094 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
3095 return cls(n
, **kwargs
)
3098 class TrivialEJA(ConcreteEJA
, RationalBasisEJA
):
3100 The trivial Euclidean Jordan algebra consisting of only a zero element.
3104 sage: from mjo.eja.eja_algebra import TrivialEJA
3108 sage: J = TrivialEJA()
3115 sage: 7*J.one()*12*J.one()
3117 sage: J.one().inner_product(J.one())
3119 sage: J.one().norm()
3121 sage: J.one().subalgebra_generated_by()
3122 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
3127 def __init__(self
, **kwargs
):
3128 jordan_product
= lambda x
,y
: x
3129 inner_product
= lambda x
,y
: 0
3132 # New defaults for keyword arguments
3133 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
3134 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
3136 super().__init
__(basis
,
3142 # The rank is zero using my definition, namely the dimension of the
3143 # largest subalgebra generated by any element.
3144 self
.rank
.set_cache(0)
3145 self
.one
.set_cache( self
.zero() )
3148 def random_instance(cls
, **kwargs
):
3149 # We don't take a "size" argument so the superclass method is
3150 # inappropriate for us.
3151 return cls(**kwargs
)
3154 class CartesianProductEJA(FiniteDimensionalEJA
):
3156 The external (orthogonal) direct sum of two or more Euclidean
3157 Jordan algebras. Every Euclidean Jordan algebra decomposes into an
3158 orthogonal direct sum of simple Euclidean Jordan algebras which is
3159 then isometric to a Cartesian product, so no generality is lost by
3160 providing only this construction.
3164 sage: from mjo.eja.eja_algebra import (random_eja,
3165 ....: CartesianProductEJA,
3167 ....: JordanSpinEJA,
3168 ....: RealSymmetricEJA)
3172 The Jordan product is inherited from our factors and implemented by
3173 our CombinatorialFreeModule Cartesian product superclass::
3175 sage: set_random_seed()
3176 sage: J1 = HadamardEJA(2)
3177 sage: J2 = RealSymmetricEJA(2)
3178 sage: J = cartesian_product([J1,J2])
3179 sage: x,y = J.random_elements(2)
3183 The ability to retrieve the original factors is implemented by our
3184 CombinatorialFreeModule Cartesian product superclass::
3186 sage: J1 = HadamardEJA(2, field=QQ)
3187 sage: J2 = JordanSpinEJA(3, field=QQ)
3188 sage: J = cartesian_product([J1,J2])
3189 sage: J.cartesian_factors()
3190 (Euclidean Jordan algebra of dimension 2 over Rational Field,
3191 Euclidean Jordan algebra of dimension 3 over Rational Field)
3193 You can provide more than two factors::
3195 sage: J1 = HadamardEJA(2)
3196 sage: J2 = JordanSpinEJA(3)
3197 sage: J3 = RealSymmetricEJA(3)
3198 sage: cartesian_product([J1,J2,J3])
3199 Euclidean Jordan algebra of dimension 2 over Algebraic Real
3200 Field (+) Euclidean Jordan algebra of dimension 3 over Algebraic
3201 Real Field (+) Euclidean Jordan algebra of dimension 6 over
3202 Algebraic Real Field
3204 Rank is additive on a Cartesian product::
3206 sage: J1 = HadamardEJA(1)
3207 sage: J2 = RealSymmetricEJA(2)
3208 sage: J = cartesian_product([J1,J2])
3209 sage: J1.rank.clear_cache()
3210 sage: J2.rank.clear_cache()
3211 sage: J.rank.clear_cache()
3214 sage: J.rank() == J1.rank() + J2.rank()
3217 The same rank computation works over the rationals, with whatever
3220 sage: J1 = HadamardEJA(1, field=QQ, orthonormalize=False)
3221 sage: J2 = RealSymmetricEJA(2, field=QQ, orthonormalize=False)
3222 sage: J = cartesian_product([J1,J2])
3223 sage: J1.rank.clear_cache()
3224 sage: J2.rank.clear_cache()
3225 sage: J.rank.clear_cache()
3228 sage: J.rank() == J1.rank() + J2.rank()
3231 The product algebra will be associative if and only if all of its
3232 components are associative::
3234 sage: J1 = HadamardEJA(2)
3235 sage: J1.is_associative()
3237 sage: J2 = HadamardEJA(3)
3238 sage: J2.is_associative()
3240 sage: J3 = RealSymmetricEJA(3)
3241 sage: J3.is_associative()
3243 sage: CP1 = cartesian_product([J1,J2])
3244 sage: CP1.is_associative()
3246 sage: CP2 = cartesian_product([J1,J3])
3247 sage: CP2.is_associative()
3250 Cartesian products of Cartesian products work::
3252 sage: J1 = JordanSpinEJA(1)
3253 sage: J2 = JordanSpinEJA(1)
3254 sage: J3 = JordanSpinEJA(1)
3255 sage: J = cartesian_product([J1,cartesian_product([J2,J3])])
3256 sage: J.multiplication_table()
3257 +----++----+----+----+
3258 | * || b0 | b1 | b2 |
3259 +====++====+====+====+
3260 | b0 || b0 | 0 | 0 |
3261 +----++----+----+----+
3262 | b1 || 0 | b1 | 0 |
3263 +----++----+----+----+
3264 | b2 || 0 | 0 | b2 |
3265 +----++----+----+----+
3266 sage: HadamardEJA(3).multiplication_table()
3267 +----++----+----+----+
3268 | * || b0 | b1 | b2 |
3269 +====++====+====+====+
3270 | b0 || b0 | 0 | 0 |
3271 +----++----+----+----+
3272 | b1 || 0 | b1 | 0 |
3273 +----++----+----+----+
3274 | b2 || 0 | 0 | b2 |
3275 +----++----+----+----+
3279 All factors must share the same base field::
3281 sage: J1 = HadamardEJA(2, field=QQ)
3282 sage: J2 = RealSymmetricEJA(2)
3283 sage: CartesianProductEJA((J1,J2))
3284 Traceback (most recent call last):
3286 ValueError: all factors must share the same base field
3288 The cached unit element is the same one that would be computed::
3290 sage: set_random_seed() # long time
3291 sage: J1 = random_eja() # long time
3292 sage: J2 = random_eja() # long time
3293 sage: J = cartesian_product([J1,J2]) # long time
3294 sage: actual = J.one() # long time
3295 sage: J.one.clear_cache() # long time
3296 sage: expected = J.one() # long time
3297 sage: actual == expected # long time
3301 Element
= FiniteDimensionalEJAElement
3304 def __init__(self
, factors
, **kwargs
):
3309 self
._sets
= factors
3311 field
= factors
[0].base_ring()
3312 if not all( J
.base_ring() == field
for J
in factors
):
3313 raise ValueError("all factors must share the same base field")
3315 associative
= all( f
.is_associative() for f
in factors
)
3317 MS
= self
.matrix_space()
3321 for b
in factors
[i
].matrix_basis():
3326 basis
= tuple( MS(b
) for b
in basis
)
3328 # Define jordan/inner products that operate on that matrix_basis.
3329 def jordan_product(x
,y
):
3331 (factors
[i
](x
[i
])*factors
[i
](y
[i
])).to_matrix()
3335 def inner_product(x
, y
):
3337 factors
[i
](x
[i
]).inner_product(factors
[i
](y
[i
]))
3341 # There's no need to check the field since it already came
3342 # from an EJA. Likewise the axioms are guaranteed to be
3343 # satisfied, unless the guy writing this class sucks.
3345 # If you want the basis to be orthonormalized, orthonormalize
3347 FiniteDimensionalEJA
.__init
__(self
,
3352 orthonormalize
=False,
3353 associative
=associative
,
3354 cartesian_product
=True,
3358 ones
= tuple(J
.one().to_matrix() for J
in factors
)
3359 self
.one
.set_cache(self(ones
))
3360 self
.rank
.set_cache(sum(J
.rank() for J
in factors
))
3362 def cartesian_factors(self
):
3363 # Copy/pasted from CombinatorialFreeModule_CartesianProduct.
3366 def cartesian_factor(self
, i
):
3368 Return the ``i``th factor of this algebra.
3370 return self
._sets
[i
]
3373 # Copy/pasted from CombinatorialFreeModule_CartesianProduct.
3374 from sage
.categories
.cartesian_product
import cartesian_product
3375 return cartesian_product
.symbol
.join("%s" % factor
3376 for factor
in self
._sets
)
3378 def matrix_space(self
):
3380 Return the space that our matrix basis lives in as a Cartesian
3385 sage: from mjo.eja.eja_algebra import (HadamardEJA,
3386 ....: RealSymmetricEJA)
3390 sage: J1 = HadamardEJA(1)
3391 sage: J2 = RealSymmetricEJA(2)
3392 sage: J = cartesian_product([J1,J2])
3393 sage: J.matrix_space()
3394 The Cartesian product of (Full MatrixSpace of 1 by 1 dense
3395 matrices over Algebraic Real Field, Full MatrixSpace of 2
3396 by 2 dense matrices over Algebraic Real Field)
3399 from sage
.categories
.cartesian_product
import cartesian_product
3400 return cartesian_product( [J
.matrix_space()
3401 for J
in self
.cartesian_factors()] )
3404 def cartesian_projection(self
, i
):
3408 sage: from mjo.eja.eja_algebra import (random_eja,
3409 ....: JordanSpinEJA,
3411 ....: RealSymmetricEJA,
3412 ....: ComplexHermitianEJA)
3416 The projection morphisms are Euclidean Jordan algebra
3419 sage: J1 = HadamardEJA(2)
3420 sage: J2 = RealSymmetricEJA(2)
3421 sage: J = cartesian_product([J1,J2])
3422 sage: J.cartesian_projection(0)
3423 Linear operator between finite-dimensional Euclidean Jordan
3424 algebras represented by the matrix:
3427 Domain: Euclidean Jordan algebra of dimension 2 over Algebraic
3428 Real Field (+) Euclidean Jordan algebra of dimension 3 over
3429 Algebraic Real Field
3430 Codomain: Euclidean Jordan algebra of dimension 2 over Algebraic
3432 sage: J.cartesian_projection(1)
3433 Linear operator between finite-dimensional Euclidean Jordan
3434 algebras represented by the matrix:
3438 Domain: Euclidean Jordan algebra of dimension 2 over Algebraic
3439 Real Field (+) Euclidean Jordan algebra of dimension 3 over
3440 Algebraic Real Field
3441 Codomain: Euclidean Jordan algebra of dimension 3 over Algebraic
3444 The projections work the way you'd expect on the vector
3445 representation of an element::
3447 sage: J1 = JordanSpinEJA(2)
3448 sage: J2 = ComplexHermitianEJA(2)
3449 sage: J = cartesian_product([J1,J2])
3450 sage: pi_left = J.cartesian_projection(0)
3451 sage: pi_right = J.cartesian_projection(1)
3452 sage: pi_left(J.one()).to_vector()
3454 sage: pi_right(J.one()).to_vector()
3456 sage: J.one().to_vector()
3461 The answer never changes::
3463 sage: set_random_seed()
3464 sage: J1 = random_eja()
3465 sage: J2 = random_eja()
3466 sage: J = cartesian_product([J1,J2])
3467 sage: P0 = J.cartesian_projection(0)
3468 sage: P1 = J.cartesian_projection(0)
3473 offset
= sum( self
.cartesian_factor(k
).dimension()
3475 Ji
= self
.cartesian_factor(i
)
3476 Pi
= self
._module
_morphism
(lambda j
: Ji
.monomial(j
- offset
),
3479 return FiniteDimensionalEJAOperator(self
,Ji
,Pi
.matrix())
3482 def cartesian_embedding(self
, i
):
3486 sage: from mjo.eja.eja_algebra import (random_eja,
3487 ....: JordanSpinEJA,
3489 ....: RealSymmetricEJA)
3493 The embedding morphisms are Euclidean Jordan algebra
3496 sage: J1 = HadamardEJA(2)
3497 sage: J2 = RealSymmetricEJA(2)
3498 sage: J = cartesian_product([J1,J2])
3499 sage: J.cartesian_embedding(0)
3500 Linear operator between finite-dimensional Euclidean Jordan
3501 algebras represented by the matrix:
3507 Domain: Euclidean Jordan algebra of dimension 2 over
3508 Algebraic Real Field
3509 Codomain: Euclidean Jordan algebra of dimension 2 over
3510 Algebraic Real Field (+) Euclidean Jordan algebra of
3511 dimension 3 over Algebraic Real Field
3512 sage: J.cartesian_embedding(1)
3513 Linear operator between finite-dimensional Euclidean Jordan
3514 algebras represented by the matrix:
3520 Domain: Euclidean Jordan algebra of dimension 3 over
3521 Algebraic Real Field
3522 Codomain: Euclidean Jordan algebra of dimension 2 over
3523 Algebraic Real Field (+) Euclidean Jordan algebra of
3524 dimension 3 over Algebraic Real Field
3526 The embeddings work the way you'd expect on the vector
3527 representation of an element::
3529 sage: J1 = JordanSpinEJA(3)
3530 sage: J2 = RealSymmetricEJA(2)
3531 sage: J = cartesian_product([J1,J2])
3532 sage: iota_left = J.cartesian_embedding(0)
3533 sage: iota_right = J.cartesian_embedding(1)
3534 sage: iota_left(J1.zero()) == J.zero()
3536 sage: iota_right(J2.zero()) == J.zero()
3538 sage: J1.one().to_vector()
3540 sage: iota_left(J1.one()).to_vector()
3542 sage: J2.one().to_vector()
3544 sage: iota_right(J2.one()).to_vector()
3546 sage: J.one().to_vector()
3551 The answer never changes::
3553 sage: set_random_seed()
3554 sage: J1 = random_eja()
3555 sage: J2 = random_eja()
3556 sage: J = cartesian_product([J1,J2])
3557 sage: E0 = J.cartesian_embedding(0)
3558 sage: E1 = J.cartesian_embedding(0)
3562 Composing a projection with the corresponding inclusion should
3563 produce the identity map, and mismatching them should produce
3566 sage: set_random_seed()
3567 sage: J1 = random_eja()
3568 sage: J2 = random_eja()
3569 sage: J = cartesian_product([J1,J2])
3570 sage: iota_left = J.cartesian_embedding(0)
3571 sage: iota_right = J.cartesian_embedding(1)
3572 sage: pi_left = J.cartesian_projection(0)
3573 sage: pi_right = J.cartesian_projection(1)
3574 sage: pi_left*iota_left == J1.one().operator()
3576 sage: pi_right*iota_right == J2.one().operator()
3578 sage: (pi_left*iota_right).is_zero()
3580 sage: (pi_right*iota_left).is_zero()
3584 offset
= sum( self
.cartesian_factor(k
).dimension()
3586 Ji
= self
.cartesian_factor(i
)
3587 Ei
= Ji
._module
_morphism
(lambda j
: self
.monomial(j
+ offset
),
3589 return FiniteDimensionalEJAOperator(Ji
,self
,Ei
.matrix())
3593 FiniteDimensionalEJA
.CartesianProduct
= CartesianProductEJA
3595 class RationalBasisCartesianProductEJA(CartesianProductEJA
,
3598 A separate class for products of algebras for which we know a
3603 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
3604 ....: RealSymmetricEJA)
3608 This gives us fast characteristic polynomial computations in
3609 product algebras, too::
3612 sage: J1 = JordanSpinEJA(2)
3613 sage: J2 = RealSymmetricEJA(3)
3614 sage: J = cartesian_product([J1,J2])
3615 sage: J.characteristic_polynomial_of().degree()
3621 def __init__(self
, algebras
, **kwargs
):
3622 CartesianProductEJA
.__init
__(self
, algebras
, **kwargs
)
3624 self
._rational
_algebra
= None
3625 if self
.vector_space().base_field() is not QQ
:
3626 self
._rational
_algebra
= cartesian_product([
3627 r
._rational
_algebra
for r
in algebras
3631 RationalBasisEJA
.CartesianProduct
= RationalBasisCartesianProductEJA
3633 def random_eja(*args
, **kwargs
):
3634 J1
= ConcreteEJA
.random_instance(*args
, **kwargs
)
3636 # This might make Cartesian products appear roughly as often as
3637 # any other ConcreteEJA.
3638 if ZZ
.random_element(len(ConcreteEJA
.__subclasses
__()) + 1) == 0:
3639 # Use random_eja() again so we can get more than two factors.
3640 J2
= random_eja(*args
, **kwargs
)
3641 J
= cartesian_product([J1
,J2
])