2 Representations and constructions for Euclidean Jordan algebras.
4 A Euclidean Jordan algebra is a Jordan algebra that has some
7 1. It is finite-dimensional.
8 2. Its scalar field is the real numbers.
9 3a. An inner product is defined on it, and...
10 3b. That inner product is compatible with the Jordan product
11 in the sense that `<x*y,z> = <y,x*z>` for all elements
12 `x,y,z` in the algebra.
14 Every Euclidean Jordan algebra is formally-real: for any two elements
15 `x` and `y` in the algebra, `x^{2} + y^{2} = 0` implies that `x = y =
16 0`. Conversely, every finite-dimensional formally-real Jordan algebra
17 can be made into a Euclidean Jordan algebra with an appropriate choice
20 Formally-real Jordan algebras were originally studied as a framework
21 for quantum mechanics. Today, Euclidean Jordan algebras are crucial in
22 symmetric cone optimization, since every symmetric cone arises as the
23 cone of squares in some Euclidean Jordan algebra.
25 It is known that every Euclidean Jordan algebra decomposes into an
26 orthogonal direct sum (essentially, a Cartesian product) of simple
27 algebras, and that moreover, up to Jordan-algebra isomorphism, there
28 are only five families of simple algebras. We provide constructions
29 for these simple algebras:
31 * :class:`BilinearFormEJA`
32 * :class:`RealSymmetricEJA`
33 * :class:`ComplexHermitianEJA`
34 * :class:`QuaternionHermitianEJA`
35 * :class:`OctonionHermitianEJA`
37 In addition to these, we provide a few other example constructions,
39 * :class:`JordanSpinEJA`
40 * :class:`HadamardEJA`
43 * :class:`ComplexSkewSymmetricEJA`
45 The Jordan spin algebra is a bilinear form algebra where the bilinear
46 form is the identity. The Hadamard EJA is simply a Cartesian product
47 of one-dimensional spin algebras. The Albert EJA is simply a special
48 case of the :class:`OctonionHermitianEJA` where the matrices are
49 three-by-three and the resulting space has dimension 27. And
50 last/least, the trivial EJA is exactly what you think it is; it could
51 also be obtained by constructing a dimension-zero instance of any of
52 the other algebras. Cartesian products of these are also supported
53 using the usual ``cartesian_product()`` function; as a result, we
54 support (up to isomorphism) all Euclidean Jordan algebras.
56 At a minimum, the following are required to construct a Euclidean
59 * A basis of matrices, column vectors, or MatrixAlgebra elements
60 * A Jordan product defined on the basis
61 * Its inner product defined on the basis
63 The real numbers form a Euclidean Jordan algebra when both the Jordan
64 and inner products are the usual multiplication. We use this as our
65 example, and demonstrate a few ways to construct an EJA.
67 First, we can use one-by-one SageMath matrices with algebraic real
68 entries to represent real numbers. We define the Jordan and inner
69 products to be essentially real-number multiplication, with the only
70 difference being that the Jordan product again returns a one-by-one
71 matrix, whereas the inner product must return a scalar. Our basis for
72 the one-by-one matrices is of course the set consisting of a single
73 matrix with its sole entry non-zero::
75 sage: from mjo.eja.eja_algebra import FiniteDimensionalEJA
76 sage: jp = lambda X,Y: X*Y
77 sage: ip = lambda X,Y: X[0,0]*Y[0,0]
78 sage: b1 = matrix(AA, [[1]])
79 sage: J1 = FiniteDimensionalEJA((b1,), jp, ip)
81 Euclidean Jordan algebra of dimension 1 over Algebraic Real Field
83 In fact, any positive scalar multiple of that inner-product would work::
85 sage: ip2 = lambda X,Y: 16*ip(X,Y)
86 sage: J2 = FiniteDimensionalEJA((b1,), jp, ip2)
88 Euclidean Jordan algebra of dimension 1 over Algebraic Real Field
90 But beware that your basis will be orthonormalized _with respect to the
91 given inner-product_ unless you pass ``orthonormalize=False`` to the
92 constructor. For example::
94 sage: J3 = FiniteDimensionalEJA((b1,), jp, ip2, orthonormalize=False)
96 Euclidean Jordan algebra of dimension 1 over Algebraic Real Field
98 To see the difference, you can take the first and only basis element
99 of the resulting algebra, and ask for it to be converted back into
102 sage: J1.basis()[0].to_matrix()
104 sage: J2.basis()[0].to_matrix()
106 sage: J3.basis()[0].to_matrix()
109 Since square roots are used in that process, the default scalar field
110 that we use is the field of algebraic real numbers, ``AA``. You can
111 also Use rational numbers, but only if you either pass
112 ``orthonormalize=False`` or know that orthonormalizing your basis
113 won't stray beyond the rational numbers. The example above would
114 have worked only because ``sqrt(16) == 4`` is rational.
116 Another option for your basis is to use elemebts of a
117 :class:`MatrixAlgebra`::
119 sage: from mjo.matrix_algebra import MatrixAlgebra
120 sage: A = MatrixAlgebra(1,AA,AA)
121 sage: J4 = FiniteDimensionalEJA(A.gens(), jp, ip)
123 Euclidean Jordan algebra of dimension 1 over Algebraic Real Field
124 sage: J4.basis()[0].to_matrix()
129 An easier way to view the entire EJA basis in its original (but
130 perhaps orthonormalized) matrix form is to use the ``matrix_basis``
133 sage: J4.matrix_basis()
138 In particular, a :class:`MatrixAlgebra` is needed to work around the
139 fact that matrices in SageMath must have entries in the same
140 (commutative and associative) ring as its scalars. There are many
141 Euclidean Jordan algebras whose elements are matrices that violate
142 those assumptions. The complex, quaternion, and octonion Hermitian
143 matrices all have entries in a ring (the complex numbers, quaternions,
144 or octonions...) that differs from the algebra's scalar ring (the real
145 numbers). Quaternions are also non-commutative; the octonions are
146 neither commutative nor associative.
150 sage: from mjo.eja.eja_algebra import random_eja
155 Euclidean Jordan algebra of dimension...
158 from sage
.algebras
.quatalg
.quaternion_algebra
import QuaternionAlgebra
159 from sage
.categories
.magmatic_algebras
import MagmaticAlgebras
160 from sage
.categories
.sets_cat
import cartesian_product
161 from sage
.combinat
.free_module
import CombinatorialFreeModule
162 from sage
.matrix
.constructor
import matrix
163 from sage
.matrix
.matrix_space
import MatrixSpace
164 from sage
.misc
.cachefunc
import cached_method
165 from sage
.misc
.table
import table
166 from sage
.modules
.free_module
import FreeModule
, VectorSpace
167 from sage
.rings
.all
import (ZZ
, QQ
, AA
, QQbar
, RR
, RLF
, CLF
,
170 from mjo
.eja
.eja_element
import (CartesianProductEJAElement
,
171 FiniteDimensionalEJAElement
)
172 from mjo
.eja
.eja_operator
import FiniteDimensionalEJAOperator
173 from mjo
.eja
.eja_utils
import _all2list
175 def EuclideanJordanAlgebras(field
):
177 The category of Euclidean Jordan algebras over ``field``, which
178 must be a subfield of the real numbers. For now this is just a
179 convenient wrapper around all of the other category axioms that
182 category
= MagmaticAlgebras(field
).FiniteDimensional()
183 category
= category
.WithBasis().Unital().Commutative()
186 class FiniteDimensionalEJA(CombinatorialFreeModule
):
188 A finite-dimensional Euclidean Jordan algebra.
192 - ``basis`` -- a tuple; a tuple of basis elements in "matrix
193 form," which must be the same form as the arguments to
194 ``jordan_product`` and ``inner_product``. In reality, "matrix
195 form" can be either vectors, matrices, or a Cartesian product
196 (ordered tuple) of vectors or matrices. All of these would
197 ideally be vector spaces in sage with no special-casing
198 needed; but in reality we turn vectors into column-matrices
199 and Cartesian products `(a,b)` into column matrices
200 `(a,b)^{T}` after converting `a` and `b` themselves.
202 - ``jordan_product`` -- a function; afunction of two ``basis``
203 elements (in matrix form) that returns their jordan product,
204 also in matrix form; this will be applied to ``basis`` to
205 compute a multiplication table for the algebra.
207 - ``inner_product`` -- a function; a function of two ``basis``
208 elements (in matrix form) that returns their inner
209 product. This will be applied to ``basis`` to compute an
210 inner-product table (basically a matrix) for this algebra.
212 - ``matrix_space`` -- the space that your matrix basis lives in,
213 or ``None`` (the default). So long as your basis does not have
214 length zero you can omit this. But in trivial algebras, it is
217 - ``field`` -- a subfield of the reals (default: ``AA``); the scalar
218 field for the algebra.
220 - ``orthonormalize`` -- boolean (default: ``True``); whether or
221 not to orthonormalize the basis. Doing so is expensive and
222 generally rules out using the rationals as your ``field``, but
223 is required for spectral decompositions.
227 sage: from mjo.eja.eja_algebra import random_eja
231 We should compute that an element subalgebra is associative even
232 if we circumvent the element method::
234 sage: J = random_eja(field=QQ,orthonormalize=False)
235 sage: x = J.random_element()
236 sage: A = x.subalgebra_generated_by(orthonormalize=False)
237 sage: basis = tuple(b.superalgebra_element() for b in A.basis())
238 sage: J.subalgebra(basis, orthonormalize=False).is_associative()
241 Element
= FiniteDimensionalEJAElement
244 def _check_input_field(field
):
245 if not field
.is_subring(RR
):
246 # Note: this does return true for the real algebraic
247 # field, the rationals, and any quadratic field where
248 # we've specified a real embedding.
249 raise ValueError("scalar field is not real")
252 def _check_input_axioms(basis
, jordan_product
, inner_product
):
253 if not all( jordan_product(bi
,bj
) == jordan_product(bj
,bi
)
256 raise ValueError("Jordan product is not commutative")
258 if not all( inner_product(bi
,bj
) == inner_product(bj
,bi
)
261 raise ValueError("inner-product is not commutative")
278 self
._check
_input
_field
(field
)
281 # Check commutativity of the Jordan and inner-products.
282 # This has to be done before we build the multiplication
283 # and inner-product tables/matrices, because we take
284 # advantage of symmetry in the process.
285 self
._check
_input
_axioms
(basis
, jordan_product
, inner_product
)
288 # All zero- and one-dimensional algebras are just the real
289 # numbers with (some positive multiples of) the usual
290 # multiplication as its Jordan and inner-product.
292 if associative
is None:
293 # We should figure it out. As with check_axioms, we have to do
294 # this without the help of the _jordan_product_is_associative()
295 # method because we need to know the category before we
296 # initialize the algebra.
297 associative
= all( jordan_product(jordan_product(bi
,bj
),bk
)
299 jordan_product(bi
,jordan_product(bj
,bk
))
304 category
= EuclideanJordanAlgebras(field
)
307 # Element subalgebras can take advantage of this.
308 category
= category
.Associative()
310 # Call the superclass constructor so that we can use its from_vector()
311 # method to build our multiplication table.
312 CombinatorialFreeModule
.__init
__(self
,
319 # Now comes all of the hard work. We'll be constructing an
320 # ambient vector space V that our (vectorized) basis lives in,
321 # as well as a subspace W of V spanned by those (vectorized)
322 # basis elements. The W-coordinates are the coefficients that
323 # we see in things like x = 1*b1 + 2*b2.
327 degree
= len(_all2list(basis
[0]))
329 # Build an ambient space that fits our matrix basis when
330 # written out as "long vectors."
331 V
= VectorSpace(field
, degree
)
333 # The matrix that will hold the orthonormal -> unorthonormal
334 # coordinate transformation. Default to an identity matrix of
335 # the appropriate size to avoid special cases for None
337 self
._deortho
_matrix
= matrix
.identity(field
,n
)
340 # Save a copy of the un-orthonormalized basis for later.
341 # Convert it to ambient V (vector) coordinates while we're
342 # at it, because we'd have to do it later anyway.
343 deortho_vector_basis
= tuple( V(_all2list(b
)) for b
in basis
)
345 from mjo
.eja
.eja_utils
import gram_schmidt
346 basis
= tuple(gram_schmidt(basis
, inner_product
))
348 # Save the (possibly orthonormalized) matrix basis for
349 # later, as well as the space that its elements live in.
350 # In most cases we can deduce the matrix space, but when
351 # n == 0 (that is, there are no basis elements) we cannot.
352 self
._matrix
_basis
= basis
353 if matrix_space
is None:
354 self
._matrix
_space
= self
._matrix
_basis
[0].parent()
356 self
._matrix
_space
= matrix_space
358 # Now create the vector space for the algebra, which will have
359 # its own set of non-ambient coordinates (in terms of the
361 vector_basis
= tuple( V(_all2list(b
)) for b
in basis
)
363 # Save the span of our matrix basis (when written out as long
364 # vectors) because otherwise we'll have to reconstruct it
365 # every time we want to coerce a matrix into the algebra.
366 self
._matrix
_span
= V
.span_of_basis( vector_basis
, check
=check_axioms
)
369 # Now "self._matrix_span" is the vector space of our
370 # algebra coordinates. The variables "X0", "X1",... refer
371 # to the entries of vectors in self._matrix_span. Thus to
372 # convert back and forth between the orthonormal
373 # coordinates and the given ones, we need to stick the
374 # original basis in self._matrix_span.
375 U
= V
.span_of_basis( deortho_vector_basis
, check
=check_axioms
)
376 self
._deortho
_matrix
= matrix
.column( U
.coordinate_vector(q
)
377 for q
in vector_basis
)
380 # Now we actually compute the multiplication and inner-product
381 # tables/matrices using the possibly-orthonormalized basis.
382 self
._inner
_product
_matrix
= matrix
.identity(field
, n
)
384 self
._multiplication
_table
= [ [zed
for j
in range(i
+1)]
387 # Note: the Jordan and inner-products are defined in terms
388 # of the ambient basis. It's important that their arguments
389 # are in ambient coordinates as well.
392 # ortho basis w.r.t. ambient coords
396 # The jordan product returns a matrixy answer, so we
397 # have to convert it to the algebra coordinates.
398 elt
= jordan_product(q_i
, q_j
)
399 elt
= self
._matrix
_span
.coordinate_vector(V(_all2list(elt
)))
400 self
._multiplication
_table
[i
][j
] = self
.from_vector(elt
)
402 if not orthonormalize
:
403 # If we're orthonormalizing the basis with respect
404 # to an inner-product, then the inner-product
405 # matrix with respect to the resulting basis is
406 # just going to be the identity.
407 ip
= inner_product(q_i
, q_j
)
408 self
._inner
_product
_matrix
[i
,j
] = ip
409 self
._inner
_product
_matrix
[j
,i
] = ip
411 self
._inner
_product
_matrix
._cache
= {'hermitian': True}
412 self
._inner
_product
_matrix
.set_immutable()
415 if not self
._is
_jordanian
():
416 raise ValueError("Jordan identity does not hold")
417 if not self
._inner
_product
_is
_associative
():
418 raise ValueError("inner product is not associative")
421 def _coerce_map_from_base_ring(self
):
423 Disable the map from the base ring into the algebra.
425 Performing a nonsense conversion like this automatically
426 is counterpedagogical. The fallback is to try the usual
427 element constructor, which should also fail.
431 sage: from mjo.eja.eja_algebra import random_eja
435 sage: J = random_eja()
437 Traceback (most recent call last):
439 ValueError: not an element of this algebra
445 def product_on_basis(self
, i
, j
):
447 Returns the Jordan product of the `i` and `j`th basis elements.
449 This completely defines the Jordan product on the algebra, and
450 is used direclty by our superclass machinery to implement
455 sage: from mjo.eja.eja_algebra import random_eja
459 sage: J = random_eja()
460 sage: n = J.dimension()
463 sage: bi_bj = J.zero()*J.zero()
465 ....: i = ZZ.random_element(n)
466 ....: j = ZZ.random_element(n)
467 ....: bi = J.monomial(i)
468 ....: bj = J.monomial(j)
469 ....: bi_bj = J.product_on_basis(i,j)
474 # We only stored the lower-triangular portion of the
475 # multiplication table.
477 return self
._multiplication
_table
[i
][j
]
479 return self
._multiplication
_table
[j
][i
]
481 def inner_product(self
, x
, y
):
483 The inner product associated with this Euclidean Jordan algebra.
485 Defaults to the trace inner product, but can be overridden by
486 subclasses if they are sure that the necessary properties are
491 sage: from mjo.eja.eja_algebra import (random_eja,
493 ....: BilinearFormEJA)
497 Our inner product is "associative," which means the following for
498 a symmetric bilinear form::
500 sage: J = random_eja()
501 sage: x,y,z = J.random_elements(3)
502 sage: (x*y).inner_product(z) == y.inner_product(x*z)
507 Ensure that this is the usual inner product for the algebras
510 sage: J = HadamardEJA.random_instance()
511 sage: x,y = J.random_elements(2)
512 sage: actual = x.inner_product(y)
513 sage: expected = x.to_vector().inner_product(y.to_vector())
514 sage: actual == expected
517 Ensure that this is one-half of the trace inner-product in a
518 BilinearFormEJA that isn't just the reals (when ``n`` isn't
519 one). This is in Faraut and Koranyi, and also my "On the
522 sage: J = BilinearFormEJA.random_instance()
523 sage: n = J.dimension()
524 sage: x = J.random_element()
525 sage: y = J.random_element()
526 sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
530 B
= self
._inner
_product
_matrix
531 return (B
*x
.to_vector()).inner_product(y
.to_vector())
534 def is_associative(self
):
536 Return whether or not this algebra's Jordan product is associative.
540 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
544 sage: J = ComplexHermitianEJA(3, field=QQ, orthonormalize=False)
545 sage: J.is_associative()
547 sage: x = sum(J.gens())
548 sage: A = x.subalgebra_generated_by(orthonormalize=False)
549 sage: A.is_associative()
553 return "Associative" in self
.category().axioms()
555 def _is_commutative(self
):
557 Whether or not this algebra's multiplication table is commutative.
559 This method should of course always return ``True``, unless
560 this algebra was constructed with ``check_axioms=False`` and
561 passed an invalid multiplication table.
563 return all( x
*y
== y
*x
for x
in self
.gens() for y
in self
.gens() )
565 def _is_jordanian(self
):
567 Whether or not this algebra's multiplication table respects the
568 Jordan identity `(x^{2})(xy) = x(x^{2}y)`.
570 We only check one arrangement of `x` and `y`, so for a
571 ``True`` result to be truly true, you should also check
572 :meth:`_is_commutative`. This method should of course always
573 return ``True``, unless this algebra was constructed with
574 ``check_axioms=False`` and passed an invalid multiplication table.
576 return all( (self
.monomial(i
)**2)*(self
.monomial(i
)*self
.monomial(j
))
578 (self
.monomial(i
))*((self
.monomial(i
)**2)*self
.monomial(j
))
579 for i
in range(self
.dimension())
580 for j
in range(self
.dimension()) )
582 def _jordan_product_is_associative(self
):
584 Return whether or not this algebra's Jordan product is
585 associative; that is, whether or not `x*(y*z) = (x*y)*z`
588 This method should agree with :meth:`is_associative` unless
589 you lied about the value of the ``associative`` parameter
590 when you constructed the algebra.
594 sage: from mjo.eja.eja_algebra import (random_eja,
595 ....: RealSymmetricEJA,
596 ....: ComplexHermitianEJA,
597 ....: QuaternionHermitianEJA)
601 sage: J = RealSymmetricEJA(4, orthonormalize=False)
602 sage: J._jordan_product_is_associative()
604 sage: x = sum(J.gens())
605 sage: A = x.subalgebra_generated_by()
606 sage: A._jordan_product_is_associative()
611 sage: J = ComplexHermitianEJA(2,field=QQ,orthonormalize=False)
612 sage: J._jordan_product_is_associative()
614 sage: x = sum(J.gens())
615 sage: A = x.subalgebra_generated_by(orthonormalize=False)
616 sage: A._jordan_product_is_associative()
621 sage: J = QuaternionHermitianEJA(2)
622 sage: J._jordan_product_is_associative()
624 sage: x = sum(J.gens())
625 sage: A = x.subalgebra_generated_by()
626 sage: A._jordan_product_is_associative()
631 The values we've presupplied to the constructors agree with
634 sage: J = random_eja()
635 sage: J.is_associative() == J._jordan_product_is_associative()
641 # Used to check whether or not something is zero.
644 # I don't know of any examples that make this magnitude
645 # necessary because I don't know how to make an
646 # associative algebra when the element subalgebra
647 # construction is unreliable (as it is over RDF; we can't
648 # find the degree of an element because we can't compute
649 # the rank of a matrix). But even multiplication of floats
650 # is non-associative, so *some* epsilon is needed... let's
651 # just take the one from _inner_product_is_associative?
654 for i
in range(self
.dimension()):
655 for j
in range(self
.dimension()):
656 for k
in range(self
.dimension()):
660 diff
= (x
*y
)*z
- x
*(y
*z
)
662 if diff
.norm() > epsilon
:
667 def _inner_product_is_associative(self
):
669 Return whether or not this algebra's inner product `B` is
670 associative; that is, whether or not `B(xy,z) = B(x,yz)`.
672 This method should of course always return ``True``, unless
673 this algebra was constructed with ``check_axioms=False`` and
674 passed an invalid Jordan or inner-product.
678 # Used to check whether or not something is zero.
681 # This choice is sufficient to allow the construction of
682 # QuaternionHermitianEJA(2, field=RDF) with check_axioms=True.
685 for i
in range(self
.dimension()):
686 for j
in range(self
.dimension()):
687 for k
in range(self
.dimension()):
691 diff
= (x
*y
).inner_product(z
) - x
.inner_product(y
*z
)
693 if diff
.abs() > epsilon
:
698 def _element_constructor_(self
, elt
):
700 Construct an element of this algebra or a subalgebra from its
701 EJA element, vector, or matrix representation.
703 This gets called only after the parent element _call_ method
704 fails to find a coercion for the argument.
708 sage: from mjo.eja.eja_algebra import (random_eja,
711 ....: RealSymmetricEJA)
715 The identity in `S^n` is converted to the identity in the EJA::
717 sage: J = RealSymmetricEJA(3)
718 sage: I = matrix.identity(QQ,3)
719 sage: J(I) == J.one()
722 This skew-symmetric matrix can't be represented in the EJA::
724 sage: J = RealSymmetricEJA(3)
725 sage: A = matrix(QQ,3, lambda i,j: i-j)
727 Traceback (most recent call last):
729 ValueError: not an element of this algebra
731 Tuples work as well, provided that the matrix basis for the
732 algebra consists of them::
734 sage: J1 = HadamardEJA(3)
735 sage: J2 = RealSymmetricEJA(2)
736 sage: J = cartesian_product([J1,J2])
737 sage: J( (J1.matrix_basis()[1], J2.matrix_basis()[2]) )
740 Subalgebra elements are embedded into the superalgebra::
742 sage: J = JordanSpinEJA(3)
745 sage: x = sum(J.gens())
746 sage: A = x.subalgebra_generated_by()
752 Ensure that we can convert any element back and forth
753 faithfully between its matrix and algebra representations::
755 sage: J = random_eja()
756 sage: x = J.random_element()
757 sage: J(x.to_matrix()) == x
760 We cannot coerce elements between algebras just because their
761 matrix representations are compatible::
763 sage: J1 = HadamardEJA(3)
764 sage: J2 = JordanSpinEJA(3)
766 Traceback (most recent call last):
768 ValueError: not an element of this algebra
770 Traceback (most recent call last):
772 ValueError: not an element of this algebra
775 msg
= "not an element of this algebra"
776 if elt
in self
.base_ring():
777 # Ensure that no base ring -> algebra coercion is performed
778 # by this method. There's some stupidity in sage that would
779 # otherwise propagate to this method; for example, sage thinks
780 # that the integer 3 belongs to the space of 2-by-2 matrices.
781 raise ValueError(msg
)
783 if hasattr(elt
, 'superalgebra_element'):
784 # Handle subalgebra elements
785 if elt
.parent().superalgebra() == self
:
786 return elt
.superalgebra_element()
788 if hasattr(elt
, 'sparse_vector'):
789 # Convert a vector into a column-matrix. We check for
790 # "sparse_vector" and not "column" because matrices also
791 # have a "column" method.
794 if elt
not in self
.matrix_space():
795 raise ValueError(msg
)
797 # Thanks for nothing! Matrix spaces aren't vector spaces in
798 # Sage, so we have to figure out its matrix-basis coordinates
799 # ourselves. We use the basis space's ring instead of the
800 # element's ring because the basis space might be an algebraic
801 # closure whereas the base ring of the 3-by-3 identity matrix
802 # could be QQ instead of QQbar.
804 # And, we also have to handle Cartesian product bases (when
805 # the matrix basis consists of tuples) here. The "good news"
806 # is that we're already converting everything to long vectors,
807 # and that strategy works for tuples as well.
809 elt
= self
._matrix
_span
.ambient_vector_space()(_all2list(elt
))
812 coords
= self
._matrix
_span
.coordinate_vector(elt
)
813 except ArithmeticError: # vector is not in free module
814 raise ValueError(msg
)
816 return self
.from_vector(coords
)
820 Return a string representation of ``self``.
824 sage: from mjo.eja.eja_algebra import JordanSpinEJA
828 Ensure that it says what we think it says::
830 sage: JordanSpinEJA(2, field=AA)
831 Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
832 sage: JordanSpinEJA(3, field=RDF)
833 Euclidean Jordan algebra of dimension 3 over Real Double Field
836 fmt
= "Euclidean Jordan algebra of dimension {} over {}"
837 return fmt
.format(self
.dimension(), self
.base_ring())
841 def characteristic_polynomial_of(self
):
843 Return the algebra's "characteristic polynomial of" function,
844 which is itself a multivariate polynomial that, when evaluated
845 at the coordinates of some algebra element, returns that
846 element's characteristic polynomial.
848 The resulting polynomial has `n+1` variables, where `n` is the
849 dimension of this algebra. The first `n` variables correspond to
850 the coordinates of an algebra element: when evaluated at the
851 coordinates of an algebra element with respect to a certain
852 basis, the result is a univariate polynomial (in the one
853 remaining variable ``t``), namely the characteristic polynomial
858 sage: from mjo.eja.eja_algebra import JordanSpinEJA, TrivialEJA
862 The characteristic polynomial in the spin algebra is given in
863 Alizadeh, Example 11.11::
865 sage: J = JordanSpinEJA(3)
866 sage: p = J.characteristic_polynomial_of(); p
867 X0^2 - X1^2 - X2^2 + (-2*t)*X0 + t^2
868 sage: xvec = J.one().to_vector()
872 By definition, the characteristic polynomial is a monic
873 degree-zero polynomial in a rank-zero algebra. Note that
874 Cayley-Hamilton is indeed satisfied since the polynomial
875 ``1`` evaluates to the identity element of the algebra on
878 sage: J = TrivialEJA()
879 sage: J.characteristic_polynomial_of()
886 # The list of coefficient polynomials a_0, a_1, a_2, ..., a_(r-1).
887 a
= self
._charpoly
_coefficients
()
889 # We go to a bit of trouble here to reorder the
890 # indeterminates, so that it's easier to evaluate the
891 # characteristic polynomial at x's coordinates and get back
892 # something in terms of t, which is what we want.
893 S
= PolynomialRing(self
.base_ring(),'t')
897 S
= PolynomialRing(S
, R
.variable_names())
900 return (t
**r
+ sum( a
[k
]*(t
**k
) for k
in range(r
) ))
902 def coordinate_polynomial_ring(self
):
904 The multivariate polynomial ring in which this algebra's
905 :meth:`characteristic_polynomial_of` lives.
909 sage: from mjo.eja.eja_algebra import (HadamardEJA,
910 ....: RealSymmetricEJA)
914 sage: J = HadamardEJA(2)
915 sage: J.coordinate_polynomial_ring()
916 Multivariate Polynomial Ring in X0, X1...
917 sage: J = RealSymmetricEJA(3,field=QQ,orthonormalize=False)
918 sage: J.coordinate_polynomial_ring()
919 Multivariate Polynomial Ring in X0, X1, X2, X3, X4, X5...
922 var_names
= tuple( "X%d" % z
for z
in range(self
.dimension()) )
923 return PolynomialRing(self
.base_ring(), var_names
)
925 def inner_product(self
, x
, y
):
927 The inner product associated with this Euclidean Jordan algebra.
929 Defaults to the trace inner product, but can be overridden by
930 subclasses if they are sure that the necessary properties are
935 sage: from mjo.eja.eja_algebra import (random_eja,
937 ....: BilinearFormEJA)
941 Our inner product is "associative," which means the following for
942 a symmetric bilinear form::
944 sage: J = random_eja()
945 sage: x,y,z = J.random_elements(3)
946 sage: (x*y).inner_product(z) == y.inner_product(x*z)
951 Ensure that this is the usual inner product for the algebras
954 sage: J = HadamardEJA.random_instance()
955 sage: x,y = J.random_elements(2)
956 sage: actual = x.inner_product(y)
957 sage: expected = x.to_vector().inner_product(y.to_vector())
958 sage: actual == expected
961 Ensure that this is one-half of the trace inner-product in a
962 BilinearFormEJA that isn't just the reals (when ``n`` isn't
963 one). This is in Faraut and Koranyi, and also my "On the
966 sage: J = BilinearFormEJA.random_instance()
967 sage: n = J.dimension()
968 sage: x = J.random_element()
969 sage: y = J.random_element()
970 sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
973 B
= self
._inner
_product
_matrix
974 return (B
*x
.to_vector()).inner_product(y
.to_vector())
977 def is_trivial(self
):
979 Return whether or not this algebra is trivial.
981 A trivial algebra contains only the zero element.
985 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
990 sage: J = ComplexHermitianEJA(3)
996 sage: J = TrivialEJA()
1001 return self
.dimension() == 0
1004 def multiplication_table(self
):
1006 Return a visual representation of this algebra's multiplication
1007 table (on basis elements).
1011 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1015 sage: J = JordanSpinEJA(4)
1016 sage: J.multiplication_table()
1017 +----++----+----+----+----+
1018 | * || b0 | b1 | b2 | b3 |
1019 +====++====+====+====+====+
1020 | b0 || b0 | b1 | b2 | b3 |
1021 +----++----+----+----+----+
1022 | b1 || b1 | b0 | 0 | 0 |
1023 +----++----+----+----+----+
1024 | b2 || b2 | 0 | b0 | 0 |
1025 +----++----+----+----+----+
1026 | b3 || b3 | 0 | 0 | b0 |
1027 +----++----+----+----+----+
1030 n
= self
.dimension()
1031 # Prepend the header row.
1032 M
= [["*"] + list(self
.gens())]
1034 # And to each subsequent row, prepend an entry that belongs to
1035 # the left-side "header column."
1036 M
+= [ [self
.monomial(i
)] + [ self
.monomial(i
)*self
.monomial(j
)
1040 return table(M
, header_row
=True, header_column
=True, frame
=True)
1043 def matrix_basis(self
):
1045 Return an (often more natural) representation of this algebras
1046 basis as an ordered tuple of matrices.
1048 Every finite-dimensional Euclidean Jordan Algebra is a, up to
1049 Jordan isomorphism, a direct sum of five simple
1050 algebras---four of which comprise Hermitian matrices. And the
1051 last type of algebra can of course be thought of as `n`-by-`1`
1052 column matrices (ambiguusly called column vectors) to avoid
1053 special cases. As a result, matrices (and column vectors) are
1054 a natural representation format for Euclidean Jordan algebra
1057 But, when we construct an algebra from a basis of matrices,
1058 those matrix representations are lost in favor of coordinate
1059 vectors *with respect to* that basis. We could eventually
1060 convert back if we tried hard enough, but having the original
1061 representations handy is valuable enough that we simply store
1062 them and return them from this method.
1064 Why implement this for non-matrix algebras? Avoiding special
1065 cases for the :class:`BilinearFormEJA` pays with simplicity in
1066 its own right. But mainly, we would like to be able to assume
1067 that elements of a :class:`CartesianProductEJA` can be displayed
1068 nicely, without having to have special classes for direct sums
1069 one of whose components was a matrix algebra.
1073 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1074 ....: RealSymmetricEJA)
1078 sage: J = RealSymmetricEJA(2)
1080 Finite family {0: b0, 1: b1, 2: b2}
1081 sage: J.matrix_basis()
1083 [1 0] [ 0 0.7071067811865475?] [0 0]
1084 [0 0], [0.7071067811865475? 0], [0 1]
1089 sage: J = JordanSpinEJA(2)
1091 Finite family {0: b0, 1: b1}
1092 sage: J.matrix_basis()
1098 return self
._matrix
_basis
1101 def matrix_space(self
):
1103 Return the matrix space in which this algebra's elements live, if
1104 we think of them as matrices (including column vectors of the
1107 "By default" this will be an `n`-by-`1` column-matrix space,
1108 except when the algebra is trivial. There it's `n`-by-`n`
1109 (where `n` is zero), to ensure that two elements of the matrix
1110 space (empty matrices) can be multiplied. For algebras of
1111 matrices, this returns the space in which their
1112 real embeddings live.
1116 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
1117 ....: JordanSpinEJA,
1118 ....: QuaternionHermitianEJA,
1123 By default, the matrix representation is just a column-matrix
1124 equivalent to the vector representation::
1126 sage: J = JordanSpinEJA(3)
1127 sage: J.matrix_space()
1128 Full MatrixSpace of 3 by 1 dense matrices over Algebraic
1131 The matrix representation in the trivial algebra is
1132 zero-by-zero instead of the usual `n`-by-one::
1134 sage: J = TrivialEJA()
1135 sage: J.matrix_space()
1136 Full MatrixSpace of 0 by 0 dense matrices over Algebraic
1139 The matrix space for complex/quaternion Hermitian matrix EJA
1140 is the space in which their real-embeddings live, not the
1141 original complex/quaternion matrix space::
1143 sage: J = ComplexHermitianEJA(2,field=QQ,orthonormalize=False)
1144 sage: J.matrix_space()
1145 Module of 2 by 2 matrices with entries in Algebraic Field over
1146 the scalar ring Rational Field
1147 sage: J = QuaternionHermitianEJA(1,field=QQ,orthonormalize=False)
1148 sage: J.matrix_space()
1149 Module of 1 by 1 matrices with entries in Quaternion
1150 Algebra (-1, -1) with base ring Rational Field over
1151 the scalar ring Rational Field
1154 return self
._matrix
_space
1160 Return the unit element of this algebra.
1164 sage: from mjo.eja.eja_algebra import (HadamardEJA,
1169 We can compute unit element in the Hadamard EJA::
1171 sage: J = HadamardEJA(5)
1173 b0 + b1 + b2 + b3 + b4
1175 The unit element in the Hadamard EJA is inherited in the
1176 subalgebras generated by its elements::
1178 sage: J = HadamardEJA(5)
1180 b0 + b1 + b2 + b3 + b4
1181 sage: x = sum(J.gens())
1182 sage: A = x.subalgebra_generated_by(orthonormalize=False)
1185 sage: A.one().superalgebra_element()
1186 b0 + b1 + b2 + b3 + b4
1190 The identity element acts like the identity, regardless of
1191 whether or not we orthonormalize::
1193 sage: J = random_eja()
1194 sage: x = J.random_element()
1195 sage: J.one()*x == x and x*J.one() == x
1197 sage: A = x.subalgebra_generated_by()
1198 sage: y = A.random_element()
1199 sage: A.one()*y == y and y*A.one() == y
1204 sage: J = random_eja(field=QQ, orthonormalize=False)
1205 sage: x = J.random_element()
1206 sage: J.one()*x == x and x*J.one() == x
1208 sage: A = x.subalgebra_generated_by(orthonormalize=False)
1209 sage: y = A.random_element()
1210 sage: A.one()*y == y and y*A.one() == y
1213 The matrix of the unit element's operator is the identity,
1214 regardless of the base field and whether or not we
1217 sage: J = random_eja()
1218 sage: actual = J.one().operator().matrix()
1219 sage: expected = matrix.identity(J.base_ring(), J.dimension())
1220 sage: actual == expected
1222 sage: x = J.random_element()
1223 sage: A = x.subalgebra_generated_by()
1224 sage: actual = A.one().operator().matrix()
1225 sage: expected = matrix.identity(A.base_ring(), A.dimension())
1226 sage: actual == expected
1231 sage: J = random_eja(field=QQ, orthonormalize=False)
1232 sage: actual = J.one().operator().matrix()
1233 sage: expected = matrix.identity(J.base_ring(), J.dimension())
1234 sage: actual == expected
1236 sage: x = J.random_element()
1237 sage: A = x.subalgebra_generated_by(orthonormalize=False)
1238 sage: actual = A.one().operator().matrix()
1239 sage: expected = matrix.identity(A.base_ring(), A.dimension())
1240 sage: actual == expected
1243 Ensure that the cached unit element (often precomputed by
1244 hand) agrees with the computed one::
1246 sage: J = random_eja()
1247 sage: cached = J.one()
1248 sage: J.one.clear_cache()
1249 sage: J.one() == cached
1254 sage: J = random_eja(field=QQ, orthonormalize=False)
1255 sage: cached = J.one()
1256 sage: J.one.clear_cache()
1257 sage: J.one() == cached
1261 # We can brute-force compute the matrices of the operators
1262 # that correspond to the basis elements of this algebra.
1263 # If some linear combination of those basis elements is the
1264 # algebra identity, then the same linear combination of
1265 # their matrices has to be the identity matrix.
1267 # Of course, matrices aren't vectors in sage, so we have to
1268 # appeal to the "long vectors" isometry.
1270 V
= VectorSpace(self
.base_ring(), self
.dimension()**2)
1271 oper_vecs
= [ V(g
.operator().matrix().list()) for g
in self
.gens() ]
1273 # Now we use basic linear algebra to find the coefficients,
1274 # of the matrices-as-vectors-linear-combination, which should
1275 # work for the original algebra basis too.
1276 A
= matrix(self
.base_ring(), oper_vecs
)
1278 # We used the isometry on the left-hand side already, but we
1279 # still need to do it for the right-hand side. Recall that we
1280 # wanted something that summed to the identity matrix.
1281 b
= V( matrix
.identity(self
.base_ring(), self
.dimension()).list() )
1283 # Now if there's an identity element in the algebra, this
1284 # should work. We solve on the left to avoid having to
1285 # transpose the matrix "A".
1286 return self
.from_vector(A
.solve_left(b
))
1289 def peirce_decomposition(self
, c
):
1291 The Peirce decomposition of this algebra relative to the
1294 In the future, this can be extended to a complete system of
1295 orthogonal idempotents.
1299 - ``c`` -- an idempotent of this algebra.
1303 A triple (J0, J5, J1) containing two subalgebras and one subspace
1306 - ``J0`` -- the algebra on the eigenspace of ``c.operator()``
1307 corresponding to the eigenvalue zero.
1309 - ``J5`` -- the eigenspace (NOT a subalgebra) of ``c.operator()``
1310 corresponding to the eigenvalue one-half.
1312 - ``J1`` -- the algebra on the eigenspace of ``c.operator()``
1313 corresponding to the eigenvalue one.
1315 These are the only possible eigenspaces for that operator, and this
1316 algebra is a direct sum of them. The spaces ``J0`` and ``J1`` are
1317 orthogonal, and are subalgebras of this algebra with the appropriate
1322 sage: from mjo.eja.eja_algebra import random_eja, RealSymmetricEJA
1326 The canonical example comes from the symmetric matrices, which
1327 decompose into diagonal and off-diagonal parts::
1329 sage: J = RealSymmetricEJA(3)
1330 sage: C = matrix(QQ, [ [1,0,0],
1334 sage: J0,J5,J1 = J.peirce_decomposition(c)
1336 Euclidean Jordan algebra of dimension 1...
1338 Vector space of degree 6 and dimension 2...
1340 Euclidean Jordan algebra of dimension 3...
1341 sage: J0.one().to_matrix()
1345 sage: orig_df = AA.options.display_format
1346 sage: AA.options.display_format = 'radical'
1347 sage: J.from_vector(J5.basis()[0]).to_matrix()
1351 sage: J.from_vector(J5.basis()[1]).to_matrix()
1355 sage: AA.options.display_format = orig_df
1356 sage: J1.one().to_matrix()
1363 Every algebra decomposes trivially with respect to its identity
1366 sage: J = random_eja()
1367 sage: J0,J5,J1 = J.peirce_decomposition(J.one())
1368 sage: J0.dimension() == 0 and J5.dimension() == 0
1370 sage: J1.superalgebra() == J and J1.dimension() == J.dimension()
1373 The decomposition is into eigenspaces, and its components are
1374 therefore necessarily orthogonal. Moreover, the identity
1375 elements in the two subalgebras are the projections onto their
1376 respective subspaces of the superalgebra's identity element::
1378 sage: J = random_eja()
1379 sage: x = J.random_element()
1380 sage: if not J.is_trivial():
1381 ....: while x.is_nilpotent():
1382 ....: x = J.random_element()
1383 sage: c = x.subalgebra_idempotent()
1384 sage: J0,J5,J1 = J.peirce_decomposition(c)
1386 sage: for (w,y,z) in zip(J0.basis(), J5.basis(), J1.basis()):
1387 ....: w = w.superalgebra_element()
1388 ....: y = J.from_vector(y)
1389 ....: z = z.superalgebra_element()
1390 ....: ipsum += w.inner_product(y).abs()
1391 ....: ipsum += w.inner_product(z).abs()
1392 ....: ipsum += y.inner_product(z).abs()
1395 sage: J1(c) == J1.one()
1397 sage: J0(J.one() - c) == J0.one()
1401 if not c
.is_idempotent():
1402 raise ValueError("element is not idempotent: %s" % c
)
1404 # Default these to what they should be if they turn out to be
1405 # trivial, because eigenspaces_left() won't return eigenvalues
1406 # corresponding to trivial spaces (e.g. it returns only the
1407 # eigenspace corresponding to lambda=1 if you take the
1408 # decomposition relative to the identity element).
1409 trivial
= self
.subalgebra((), check_axioms
=False)
1410 J0
= trivial
# eigenvalue zero
1411 J5
= VectorSpace(self
.base_ring(), 0) # eigenvalue one-half
1412 J1
= trivial
# eigenvalue one
1414 for (eigval
, eigspace
) in c
.operator().matrix().right_eigenspaces():
1415 if eigval
== ~
(self
.base_ring()(2)):
1418 gens
= tuple( self
.from_vector(b
) for b
in eigspace
.basis() )
1419 subalg
= self
.subalgebra(gens
, check_axioms
=False)
1425 raise ValueError("unexpected eigenvalue: %s" % eigval
)
1430 def random_element(self
, thorough
=False):
1432 Return a random element of this algebra.
1434 Our algebra superclass method only returns a linear
1435 combination of at most two basis elements. We instead
1436 want the vector space "random element" method that
1437 returns a more diverse selection.
1441 - ``thorough`` -- (boolean; default False) whether or not we
1442 should generate irrational coefficients for the random
1443 element when our base ring is irrational; this slows the
1444 algebra operations to a crawl, but any truly random method
1448 # For a general base ring... maybe we can trust this to do the
1449 # right thing? Unlikely, but.
1450 V
= self
.vector_space()
1451 if self
.base_ring() is AA
and not thorough
:
1452 # Now that AA generates actually random random elements
1453 # (post Trac 30875), we only need to de-thorough the
1454 # randomness when asked to.
1455 V
= V
.change_ring(QQ
)
1457 v
= V
.random_element()
1458 return self
.from_vector(V
.coordinate_vector(v
))
1460 def random_elements(self
, count
, thorough
=False):
1462 Return ``count`` random elements as a tuple.
1466 - ``thorough`` -- (boolean; default False) whether or not we
1467 should generate irrational coefficients for the random
1468 elements when our base ring is irrational; this slows the
1469 algebra operations to a crawl, but any truly random method
1474 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1478 sage: J = JordanSpinEJA(3)
1479 sage: x,y,z = J.random_elements(3)
1480 sage: all( [ x in J, y in J, z in J ])
1482 sage: len( J.random_elements(10) ) == 10
1486 return tuple( self
.random_element(thorough
)
1487 for idx
in range(count
) )
1490 def operator_polynomial_matrix(self
):
1492 Return the matrix of polynomials (over this algebra's
1493 :meth:`coordinate_polynomial_ring`) that, when evaluated at
1494 the basis coordinates of an element `x`, produces the basis
1495 representation of `L_{x}`.
1499 sage: from mjo.eja.eja_algebra import (HadamardEJA,
1500 ....: JordanSpinEJA)
1504 sage: J = HadamardEJA(4)
1505 sage: L_x = J.operator_polynomial_matrix()
1512 sage: d = zip(J.coordinate_polynomial_ring().gens(), x.to_vector())
1513 sage: L_x.subs(dict(d))
1521 sage: J = JordanSpinEJA(4)
1522 sage: L_x = J.operator_polynomial_matrix()
1529 sage: d = zip(J.coordinate_polynomial_ring().gens(), x.to_vector())
1530 sage: L_x.subs(dict(d))
1537 R
= self
.coordinate_polynomial_ring()
1540 # From a result in my book, these are the entries of the
1541 # basis representation of L_x.
1542 return sum( v
*self
.monomial(k
).operator().matrix()[i
,j
]
1543 for (k
,v
) in enumerate(R
.gens()) )
1545 n
= self
.dimension()
1546 return matrix(R
, n
, n
, L_x_i_j
)
1549 def _charpoly_coefficients(self
):
1551 The `r` polynomial coefficients of the "characteristic polynomial
1556 sage: from mjo.eja.eja_algebra import random_eja
1560 The theory shows that these are all homogeneous polynomials of
1563 sage: J = random_eja()
1564 sage: all(p.is_homogeneous() for p in J._charpoly_coefficients())
1568 n
= self
.dimension()
1569 R
= self
.coordinate_polynomial_ring()
1570 F
= R
.fraction_field()
1572 L_x
= self
.operator_polynomial_matrix()
1575 if self
.rank
.is_in_cache():
1577 # There's no need to pad the system with redundant
1578 # columns if we *know* they'll be redundant.
1581 # Compute an extra power in case the rank is equal to
1582 # the dimension (otherwise, we would stop at x^(r-1)).
1583 x_powers
= [ (L_x
**k
)*self
.one().to_vector()
1584 for k
in range(n
+1) ]
1585 A
= matrix
.column(F
, x_powers
[:n
])
1586 AE
= A
.extended_echelon_form()
1593 # The theory says that only the first "r" coefficients are
1594 # nonzero, and they actually live in the original polynomial
1595 # ring and not the fraction field. We negate them because in
1596 # the actual characteristic polynomial, they get moved to the
1597 # other side where x^r lives. We don't bother to trim A_rref
1598 # down to a square matrix and solve the resulting system,
1599 # because the upper-left r-by-r portion of A_rref is
1600 # guaranteed to be the identity matrix, so e.g.
1602 # A_rref.solve_right(Y)
1604 # would just be returning Y.
1605 return (-E
*b
)[:r
].change_ring(R
)
1610 Return the rank of this EJA.
1612 This is a cached method because we know the rank a priori for
1613 all of the algebras we can construct. Thus we can avoid the
1614 expensive ``_charpoly_coefficients()`` call unless we truly
1615 need to compute the whole characteristic polynomial.
1619 sage: from mjo.eja.eja_algebra import (HadamardEJA,
1620 ....: JordanSpinEJA,
1621 ....: RealSymmetricEJA,
1622 ....: ComplexHermitianEJA,
1623 ....: QuaternionHermitianEJA,
1628 The rank of the Jordan spin algebra is always two::
1630 sage: JordanSpinEJA(2).rank()
1632 sage: JordanSpinEJA(3).rank()
1634 sage: JordanSpinEJA(4).rank()
1637 The rank of the `n`-by-`n` Hermitian real, complex, or
1638 quaternion matrices is `n`::
1640 sage: RealSymmetricEJA(4).rank()
1642 sage: ComplexHermitianEJA(3).rank()
1644 sage: QuaternionHermitianEJA(2).rank()
1649 Ensure that every EJA that we know how to construct has a
1650 positive integer rank, unless the algebra is trivial in
1651 which case its rank will be zero::
1653 sage: J = random_eja()
1657 sage: r > 0 or (r == 0 and J.is_trivial())
1660 Ensure that computing the rank actually works, since the ranks
1661 of all simple algebras are known and will be cached by default::
1663 sage: J = random_eja() # long time
1664 sage: cached = J.rank() # long time
1665 sage: J.rank.clear_cache() # long time
1666 sage: J.rank() == cached # long time
1670 return len(self
._charpoly
_coefficients
())
1673 def subalgebra(self
, basis
, **kwargs
):
1675 Create a subalgebra of this algebra from the given basis.
1677 from mjo
.eja
.eja_subalgebra
import FiniteDimensionalEJASubalgebra
1678 return FiniteDimensionalEJASubalgebra(self
, basis
, **kwargs
)
1681 def vector_space(self
):
1683 Return the vector space that underlies this algebra.
1687 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1691 sage: J = RealSymmetricEJA(2)
1692 sage: J.vector_space()
1693 Vector space of dimension 3 over...
1696 return self
.zero().to_vector().parent().ambient_vector_space()
1700 class RationalBasisEJA(FiniteDimensionalEJA
):
1702 Algebras whose supplied basis elements have all rational entries.
1706 sage: from mjo.eja.eja_algebra import BilinearFormEJA
1710 The supplied basis is orthonormalized by default::
1712 sage: B = matrix(QQ, [[1, 0, 0], [0, 25, -32], [0, -32, 41]])
1713 sage: J = BilinearFormEJA(B)
1714 sage: J.matrix_basis()
1731 # Abuse the check_field parameter to check that the entries of
1732 # out basis (in ambient coordinates) are in the field QQ.
1733 # Use _all2list to get the vector coordinates of octonion
1734 # entries and not the octonions themselves (which are not
1736 if not all( all(b_i
in QQ
for b_i
in _all2list(b
))
1738 raise TypeError("basis not rational")
1740 super().__init
__(basis
,
1744 check_field
=check_field
,
1747 self
._rational
_algebra
= None
1749 # There's no point in constructing the extra algebra if this
1750 # one is already rational.
1752 # Note: the same Jordan and inner-products work here,
1753 # because they are necessarily defined with respect to
1754 # ambient coordinates and not any particular basis.
1755 self
._rational
_algebra
= FiniteDimensionalEJA(
1760 matrix_space
=self
.matrix_space(),
1761 associative
=self
.is_associative(),
1762 orthonormalize
=False,
1766 def rational_algebra(self
):
1767 # Using None as a flag here (rather than just assigning "self"
1768 # to self._rational_algebra by default) feels a little bit
1769 # more sane to me in a garbage-collected environment.
1770 if self
._rational
_algebra
is None:
1773 return self
._rational
_algebra
1776 def _charpoly_coefficients(self
):
1780 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
1781 ....: JordanSpinEJA)
1785 The base ring of the resulting polynomial coefficients is what
1786 it should be, and not the rationals (unless the algebra was
1787 already over the rationals)::
1789 sage: J = JordanSpinEJA(3)
1790 sage: J._charpoly_coefficients()
1791 (X0^2 - X1^2 - X2^2, -2*X0)
1792 sage: a0 = J._charpoly_coefficients()[0]
1794 Algebraic Real Field
1795 sage: a0.base_ring()
1796 Algebraic Real Field
1799 if self
.rational_algebra() is self
:
1800 # Bypass the hijinks if they won't benefit us.
1801 return super()._charpoly
_coefficients
()
1803 # Do the computation over the rationals. The answer will be
1804 # the same, because all we've done is a change of basis.
1805 # Then, change back from QQ to our real base ring
1806 a
= ( a_i
.change_ring(self
.base_ring())
1807 for a_i
in self
.rational_algebra()._charpoly
_coefficients
() )
1809 # Otherwise, convert the coordinate variables back to the
1810 # deorthonormalized ones.
1811 R
= self
.coordinate_polynomial_ring()
1812 from sage
.modules
.free_module_element
import vector
1813 X
= vector(R
, R
.gens())
1814 BX
= self
._deortho
_matrix
*X
1816 subs_dict
= { X[i]: BX[i] for i in range(len(X)) }
1817 return tuple( a_i
.subs(subs_dict
) for a_i
in a
)
1819 class ConcreteEJA(FiniteDimensionalEJA
):
1821 A class for the Euclidean Jordan algebras that we know by name.
1823 These are the Jordan algebras whose basis, multiplication table,
1824 rank, and so on are known a priori. More to the point, they are
1825 the Euclidean Jordan algebras for which we are able to conjure up
1826 a "random instance."
1830 sage: from mjo.eja.eja_algebra import ConcreteEJA
1834 Our basis is normalized with respect to the algebra's inner
1835 product, unless we specify otherwise::
1837 sage: J = ConcreteEJA.random_instance()
1838 sage: all( b.norm() == 1 for b in J.gens() )
1841 Since our basis is orthonormal with respect to the algebra's inner
1842 product, and since we know that this algebra is an EJA, any
1843 left-multiplication operator's matrix will be symmetric because
1844 natural->EJA basis representation is an isometry and within the
1845 EJA the operator is self-adjoint by the Jordan axiom::
1847 sage: J = ConcreteEJA.random_instance()
1848 sage: x = J.random_element()
1849 sage: x.operator().is_self_adjoint()
1854 def _max_random_instance_dimension():
1856 The maximum dimension of any random instance. Ten dimensions seems
1857 to be about the point where everything takes a turn for the
1858 worse. And dimension ten (but not nine) allows the 4-by-4 real
1859 Hermitian matrices, the 2-by-2 quaternion Hermitian matrices,
1860 and the 2-by-2 octonion Hermitian matrices.
1865 def _max_random_instance_size(max_dimension
):
1867 Return an integer "size" that is an upper bound on the size of
1868 this algebra when it is used in a random test case. This size
1869 (which can be passed to the algebra's constructor) is itself
1870 based on the ``max_dimension`` parameter.
1872 This method must be implemented in each subclass.
1874 raise NotImplementedError
1877 def random_instance(cls
, max_dimension
=None, *args
, **kwargs
):
1879 Return a random instance of this type of algebra whose dimension
1880 is less than or equal to the lesser of ``max_dimension`` and
1881 the value returned by ``_max_random_instance_dimension()``. If
1882 the dimension bound is omitted, then only the
1883 ``_max_random_instance_dimension()`` is used as a bound.
1885 This method should be implemented in each subclass.
1889 sage: from mjo.eja.eja_algebra import ConcreteEJA
1893 Both the class bound and the ``max_dimension`` argument are upper
1894 bounds on the dimension of the algebra returned::
1896 sage: from sage.misc.prandom import choice
1897 sage: eja_class = choice(ConcreteEJA.__subclasses__())
1898 sage: class_max_d = eja_class._max_random_instance_dimension()
1899 sage: J = eja_class.random_instance(max_dimension=20,
1901 ....: orthonormalize=False)
1902 sage: J.dimension() <= class_max_d
1904 sage: J = eja_class.random_instance(max_dimension=2,
1906 ....: orthonormalize=False)
1907 sage: J.dimension() <= 2
1911 from sage
.misc
.prandom
import choice
1912 eja_class
= choice(cls
.__subclasses
__())
1914 # These all bubble up to the RationalBasisEJA superclass
1915 # constructor, so any (kw)args valid there are also valid
1917 return eja_class
.random_instance(max_dimension
, *args
, **kwargs
)
1920 class HermitianMatrixEJA(FiniteDimensionalEJA
):
1922 def _denormalized_basis(A
):
1924 Returns a basis for the given Hermitian matrix space.
1926 Why do we embed these? Basically, because all of numerical linear
1927 algebra assumes that you're working with vectors consisting of `n`
1928 entries from a field and scalars from the same field. There's no way
1929 to tell SageMath that (for example) the vectors contain complex
1930 numbers, while the scalar field is real.
1934 sage: from mjo.hurwitz import (ComplexMatrixAlgebra,
1935 ....: QuaternionMatrixAlgebra,
1936 ....: OctonionMatrixAlgebra)
1937 sage: from mjo.eja.eja_algebra import HermitianMatrixEJA
1941 sage: n = ZZ.random_element(1,5)
1942 sage: A = MatrixSpace(QQ, n)
1943 sage: B = HermitianMatrixEJA._denormalized_basis(A)
1944 sage: all( M.is_hermitian() for M in B)
1949 sage: n = ZZ.random_element(1,5)
1950 sage: A = ComplexMatrixAlgebra(n, scalars=QQ)
1951 sage: B = HermitianMatrixEJA._denormalized_basis(A)
1952 sage: all( M.is_hermitian() for M in B)
1957 sage: n = ZZ.random_element(1,5)
1958 sage: A = QuaternionMatrixAlgebra(n, scalars=QQ)
1959 sage: B = HermitianMatrixEJA._denormalized_basis(A)
1960 sage: all( M.is_hermitian() for M in B )
1965 sage: n = ZZ.random_element(1,5)
1966 sage: A = OctonionMatrixAlgebra(n, scalars=QQ)
1967 sage: B = HermitianMatrixEJA._denormalized_basis(A)
1968 sage: all( M.is_hermitian() for M in B )
1972 # These work for real MatrixSpace, whose monomials only have
1973 # two coordinates (because the last one would always be "1").
1974 es
= A
.base_ring().gens()
1975 gen
= lambda A
,m
: A
.monomial(m
[:2])
1977 if hasattr(A
, 'entry_algebra_gens'):
1978 # We've got a MatrixAlgebra, and its monomials will have
1979 # three coordinates.
1980 es
= A
.entry_algebra_gens()
1981 gen
= lambda A
,m
: A
.monomial(m
)
1984 for i
in range(A
.nrows()):
1985 for j
in range(i
+1):
1987 E_ii
= gen(A
, (i
,j
,es
[0]))
1991 E_ij
= gen(A
, (i
,j
,e
))
1992 E_ij
+= E_ij
.conjugate_transpose()
1995 return tuple( basis
)
1998 def jordan_product(X
,Y
):
1999 return (X
*Y
+ Y
*X
)/2
2002 def trace_inner_product(X
,Y
):
2004 A trace inner-product for matrices that aren't embedded in the
2005 reals. It takes MATRICES as arguments, not EJA elements.
2009 sage: from mjo.eja.eja_algebra import (RealSymmetricEJA,
2010 ....: ComplexHermitianEJA,
2011 ....: QuaternionHermitianEJA,
2012 ....: OctonionHermitianEJA)
2016 sage: J = RealSymmetricEJA(2,field=QQ,orthonormalize=False)
2017 sage: I = J.one().to_matrix()
2018 sage: J.trace_inner_product(I, -I)
2023 sage: J = ComplexHermitianEJA(2,field=QQ,orthonormalize=False)
2024 sage: I = J.one().to_matrix()
2025 sage: J.trace_inner_product(I, -I)
2030 sage: J = QuaternionHermitianEJA(2,field=QQ,orthonormalize=False)
2031 sage: I = J.one().to_matrix()
2032 sage: J.trace_inner_product(I, -I)
2037 sage: J = OctonionHermitianEJA(2,field=QQ,orthonormalize=False)
2038 sage: I = J.one().to_matrix()
2039 sage: J.trace_inner_product(I, -I)
2044 if hasattr(tr
, 'coefficient'):
2045 # Works for octonions, and has to come first because they
2046 # also have a "real()" method that doesn't return an
2047 # element of the scalar ring.
2048 return tr
.coefficient(0)
2049 elif hasattr(tr
, 'coefficient_tuple'):
2050 # Works for quaternions.
2051 return tr
.coefficient_tuple()[0]
2053 # Works for real and complex numbers.
2057 def __init__(self
, matrix_space
, **kwargs
):
2058 # We know this is a valid EJA, but will double-check
2059 # if the user passes check_axioms=True.
2060 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2062 super().__init
__(self
._denormalized
_basis
(matrix_space
),
2063 self
.jordan_product
,
2064 self
.trace_inner_product
,
2065 field
=matrix_space
.base_ring(),
2066 matrix_space
=matrix_space
,
2069 self
.rank
.set_cache(matrix_space
.nrows())
2070 self
.one
.set_cache( self(matrix_space
.one()) )
2072 class RealSymmetricEJA(HermitianMatrixEJA
, RationalBasisEJA
, ConcreteEJA
):
2074 The rank-n simple EJA consisting of real symmetric n-by-n
2075 matrices, the usual symmetric Jordan product, and the trace inner
2076 product. It has dimension `(n^2 + n)/2` over the reals.
2080 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
2084 sage: J = RealSymmetricEJA(2)
2085 sage: b0, b1, b2 = J.gens()
2093 In theory, our "field" can be any subfield of the reals::
2095 sage: RealSymmetricEJA(2, field=RDF, check_axioms=True)
2096 Euclidean Jordan algebra of dimension 3 over Real Double Field
2097 sage: RealSymmetricEJA(2, field=RR, check_axioms=True)
2098 Euclidean Jordan algebra of dimension 3 over Real Field with
2099 53 bits of precision
2103 The dimension of this algebra is `(n^2 + n) / 2`::
2105 sage: d = RealSymmetricEJA._max_random_instance_dimension()
2106 sage: n = RealSymmetricEJA._max_random_instance_size(d)
2107 sage: J = RealSymmetricEJA(n)
2108 sage: J.dimension() == (n^2 + n)/2
2111 The Jordan multiplication is what we think it is::
2113 sage: J = RealSymmetricEJA.random_instance()
2114 sage: x,y = J.random_elements(2)
2115 sage: actual = (x*y).to_matrix()
2116 sage: X = x.to_matrix()
2117 sage: Y = y.to_matrix()
2118 sage: expected = (X*Y + Y*X)/2
2119 sage: actual == expected
2121 sage: J(expected) == x*y
2124 We can change the generator prefix::
2126 sage: RealSymmetricEJA(3, prefix='q').gens()
2127 (q0, q1, q2, q3, q4, q5)
2129 We can construct the (trivial) algebra of rank zero::
2131 sage: RealSymmetricEJA(0)
2132 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2136 def _max_random_instance_size(max_dimension
):
2137 # Obtained by solving d = (n^2 + n)/2.
2138 # The ZZ-int-ZZ thing is just "floor."
2139 return ZZ(int(ZZ(8*max_dimension
+ 1).sqrt()/2 - 1/2))
2142 def random_instance(cls
, max_dimension
=None, *args
, **kwargs
):
2144 Return a random instance of this type of algebra.
2146 class_max_d
= cls
._max
_random
_instance
_dimension
()
2147 if (max_dimension
is None or max_dimension
> class_max_d
):
2148 max_dimension
= class_max_d
2149 max_size
= cls
._max
_random
_instance
_size
(max_dimension
)
2150 n
= ZZ
.random_element(max_size
+ 1)
2151 return cls(n
, **kwargs
)
2153 def __init__(self
, n
, field
=AA
, **kwargs
):
2154 A
= MatrixSpace(field
, n
)
2155 super().__init
__(A
, **kwargs
)
2157 from mjo
.eja
.eja_cache
import real_symmetric_eja_coeffs
2158 a
= real_symmetric_eja_coeffs(self
)
2160 self
.rational_algebra()._charpoly
_coefficients
.set_cache(a
)
2164 class ComplexHermitianEJA(HermitianMatrixEJA
, RationalBasisEJA
, ConcreteEJA
):
2166 The rank-n simple EJA consisting of complex Hermitian n-by-n
2167 matrices over the real numbers, the usual symmetric Jordan product,
2168 and the real-part-of-trace inner product. It has dimension `n^2` over
2173 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
2177 In theory, our "field" can be any subfield of the reals, but we
2178 can't use inexact real fields at the moment because SageMath
2179 doesn't know how to convert their elements into complex numbers,
2180 or even into algebraic reals::
2183 Traceback (most recent call last):
2185 TypeError: Illegal initializer for algebraic number
2187 Traceback (most recent call last):
2189 TypeError: Illegal initializer for algebraic number
2191 This causes the following error when we try to scale a matrix of
2192 complex numbers by an inexact real number::
2194 sage: ComplexHermitianEJA(2,field=RR)
2195 Traceback (most recent call last):
2197 TypeError: Unable to coerce entries (=(1.00000000000000,
2198 -0.000000000000000)) to coefficients in Algebraic Real Field
2202 The dimension of this algebra is `n^2`::
2204 sage: d = ComplexHermitianEJA._max_random_instance_dimension()
2205 sage: n = ComplexHermitianEJA._max_random_instance_size(d)
2206 sage: J = ComplexHermitianEJA(n)
2207 sage: J.dimension() == n^2
2210 The Jordan multiplication is what we think it is::
2212 sage: J = ComplexHermitianEJA.random_instance()
2213 sage: x,y = J.random_elements(2)
2214 sage: actual = (x*y).to_matrix()
2215 sage: X = x.to_matrix()
2216 sage: Y = y.to_matrix()
2217 sage: expected = (X*Y + Y*X)/2
2218 sage: actual == expected
2220 sage: J(expected) == x*y
2223 We can change the generator prefix::
2225 sage: ComplexHermitianEJA(2, prefix='z').gens()
2228 We can construct the (trivial) algebra of rank zero::
2230 sage: ComplexHermitianEJA(0)
2231 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2234 def __init__(self
, n
, field
=AA
, **kwargs
):
2235 from mjo
.hurwitz
import ComplexMatrixAlgebra
2236 A
= ComplexMatrixAlgebra(n
, scalars
=field
)
2237 super().__init
__(A
, **kwargs
)
2239 from mjo
.eja
.eja_cache
import complex_hermitian_eja_coeffs
2240 a
= complex_hermitian_eja_coeffs(self
)
2242 self
.rational_algebra()._charpoly
_coefficients
.set_cache(a
)
2245 def _max_random_instance_size(max_dimension
):
2246 # Obtained by solving d = n^2.
2247 # The ZZ-int-ZZ thing is just "floor."
2248 return ZZ(int(ZZ(max_dimension
).sqrt()))
2251 def random_instance(cls
, max_dimension
=None, *args
, **kwargs
):
2253 Return a random instance of this type of algebra.
2255 class_max_d
= cls
._max
_random
_instance
_dimension
()
2256 if (max_dimension
is None or max_dimension
> class_max_d
):
2257 max_dimension
= class_max_d
2258 max_size
= cls
._max
_random
_instance
_size
(max_dimension
)
2259 n
= ZZ
.random_element(max_size
+ 1)
2260 return cls(n
, **kwargs
)
2263 class QuaternionHermitianEJA(HermitianMatrixEJA
, RationalBasisEJA
, ConcreteEJA
):
2265 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
2266 matrices, the usual symmetric Jordan product, and the
2267 real-part-of-trace inner product. It has dimension `2n^2 - n` over
2272 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2276 In theory, our "field" can be any subfield of the reals::
2278 sage: QuaternionHermitianEJA(2, field=RDF, check_axioms=True)
2279 Euclidean Jordan algebra of dimension 6 over Real Double Field
2280 sage: QuaternionHermitianEJA(2, field=RR, check_axioms=True)
2281 Euclidean Jordan algebra of dimension 6 over Real Field with
2282 53 bits of precision
2286 The dimension of this algebra is `2*n^2 - n`::
2288 sage: d = QuaternionHermitianEJA._max_random_instance_dimension()
2289 sage: n = QuaternionHermitianEJA._max_random_instance_size(d)
2290 sage: J = QuaternionHermitianEJA(n)
2291 sage: J.dimension() == 2*(n^2) - n
2294 The Jordan multiplication is what we think it is::
2296 sage: J = QuaternionHermitianEJA.random_instance()
2297 sage: x,y = J.random_elements(2)
2298 sage: actual = (x*y).to_matrix()
2299 sage: X = x.to_matrix()
2300 sage: Y = y.to_matrix()
2301 sage: expected = (X*Y + Y*X)/2
2302 sage: actual == expected
2304 sage: J(expected) == x*y
2307 We can change the generator prefix::
2309 sage: QuaternionHermitianEJA(2, prefix='a').gens()
2310 (a0, a1, a2, a3, a4, a5)
2312 We can construct the (trivial) algebra of rank zero::
2314 sage: QuaternionHermitianEJA(0)
2315 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2318 def __init__(self
, n
, field
=AA
, **kwargs
):
2319 from mjo
.hurwitz
import QuaternionMatrixAlgebra
2320 A
= QuaternionMatrixAlgebra(n
, scalars
=field
)
2321 super().__init
__(A
, **kwargs
)
2323 from mjo
.eja
.eja_cache
import quaternion_hermitian_eja_coeffs
2324 a
= quaternion_hermitian_eja_coeffs(self
)
2326 self
.rational_algebra()._charpoly
_coefficients
.set_cache(a
)
2331 def _max_random_instance_size(max_dimension
):
2333 The maximum rank of a random QuaternionHermitianEJA.
2335 # Obtained by solving d = 2n^2 - n.
2336 # The ZZ-int-ZZ thing is just "floor."
2337 return ZZ(int(ZZ(8*max_dimension
+ 1).sqrt()/4 + 1/4))
2340 def random_instance(cls
, max_dimension
=None, *args
, **kwargs
):
2342 Return a random instance of this type of algebra.
2344 class_max_d
= cls
._max
_random
_instance
_dimension
()
2345 if (max_dimension
is None or max_dimension
> class_max_d
):
2346 max_dimension
= class_max_d
2347 max_size
= cls
._max
_random
_instance
_size
(max_dimension
)
2348 n
= ZZ
.random_element(max_size
+ 1)
2349 return cls(n
, **kwargs
)
2351 class OctonionHermitianEJA(HermitianMatrixEJA
, RationalBasisEJA
, ConcreteEJA
):
2355 sage: from mjo.eja.eja_algebra import (FiniteDimensionalEJA,
2356 ....: OctonionHermitianEJA)
2357 sage: from mjo.hurwitz import Octonions, OctonionMatrixAlgebra
2361 The 3-by-3 algebra satisfies the axioms of an EJA::
2363 sage: OctonionHermitianEJA(3, # long time
2364 ....: field=QQ, # long time
2365 ....: orthonormalize=False, # long time
2366 ....: check_axioms=True) # long time
2367 Euclidean Jordan algebra of dimension 27 over Rational Field
2369 After a change-of-basis, the 2-by-2 algebra has the same
2370 multiplication table as the ten-dimensional Jordan spin algebra::
2372 sage: A = OctonionMatrixAlgebra(2,Octonions(QQ),QQ)
2373 sage: b = OctonionHermitianEJA._denormalized_basis(A)
2374 sage: basis = (b[0] + b[9],) + b[1:9] + (b[0] - b[9],)
2375 sage: jp = OctonionHermitianEJA.jordan_product
2376 sage: ip = OctonionHermitianEJA.trace_inner_product
2377 sage: J = FiniteDimensionalEJA(basis,
2381 ....: orthonormalize=False)
2382 sage: J.multiplication_table()
2383 +----++----+----+----+----+----+----+----+----+----+----+
2384 | * || b0 | b1 | b2 | b3 | b4 | b5 | b6 | b7 | b8 | b9 |
2385 +====++====+====+====+====+====+====+====+====+====+====+
2386 | b0 || b0 | b1 | b2 | b3 | b4 | b5 | b6 | b7 | b8 | b9 |
2387 +----++----+----+----+----+----+----+----+----+----+----+
2388 | b1 || b1 | b0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
2389 +----++----+----+----+----+----+----+----+----+----+----+
2390 | b2 || b2 | 0 | b0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
2391 +----++----+----+----+----+----+----+----+----+----+----+
2392 | b3 || b3 | 0 | 0 | b0 | 0 | 0 | 0 | 0 | 0 | 0 |
2393 +----++----+----+----+----+----+----+----+----+----+----+
2394 | b4 || b4 | 0 | 0 | 0 | b0 | 0 | 0 | 0 | 0 | 0 |
2395 +----++----+----+----+----+----+----+----+----+----+----+
2396 | b5 || b5 | 0 | 0 | 0 | 0 | b0 | 0 | 0 | 0 | 0 |
2397 +----++----+----+----+----+----+----+----+----+----+----+
2398 | b6 || b6 | 0 | 0 | 0 | 0 | 0 | b0 | 0 | 0 | 0 |
2399 +----++----+----+----+----+----+----+----+----+----+----+
2400 | b7 || b7 | 0 | 0 | 0 | 0 | 0 | 0 | b0 | 0 | 0 |
2401 +----++----+----+----+----+----+----+----+----+----+----+
2402 | b8 || b8 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | b0 | 0 |
2403 +----++----+----+----+----+----+----+----+----+----+----+
2404 | b9 || b9 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | b0 |
2405 +----++----+----+----+----+----+----+----+----+----+----+
2409 We can actually construct the 27-dimensional Albert algebra,
2410 and we get the right unit element if we recompute it::
2412 sage: J = OctonionHermitianEJA(3, # long time
2413 ....: field=QQ, # long time
2414 ....: orthonormalize=False) # long time
2415 sage: J.one.clear_cache() # long time
2416 sage: J.one() # long time
2418 sage: J.one().to_matrix() # long time
2427 The 2-by-2 algebra is isomorphic to the ten-dimensional Jordan
2428 spin algebra, but just to be sure, we recompute its rank::
2430 sage: J = OctonionHermitianEJA(2, # long time
2431 ....: field=QQ, # long time
2432 ....: orthonormalize=False) # long time
2433 sage: J.rank.clear_cache() # long time
2434 sage: J.rank() # long time
2439 def _max_random_instance_size(max_dimension
):
2441 The maximum rank of a random OctonionHermitianEJA.
2443 # There's certainly a formula for this, but with only four
2444 # cases to worry about, I'm not that motivated to derive it.
2445 if max_dimension
>= 27:
2447 elif max_dimension
>= 10:
2449 elif max_dimension
>= 1:
2455 def random_instance(cls
, max_dimension
=None, *args
, **kwargs
):
2457 Return a random instance of this type of algebra.
2459 class_max_d
= cls
._max
_random
_instance
_dimension
()
2460 if (max_dimension
is None or max_dimension
> class_max_d
):
2461 max_dimension
= class_max_d
2462 max_size
= cls
._max
_random
_instance
_size
(max_dimension
)
2463 n
= ZZ
.random_element(max_size
+ 1)
2464 return cls(n
, **kwargs
)
2466 def __init__(self
, n
, field
=AA
, **kwargs
):
2468 # Otherwise we don't get an EJA.
2469 raise ValueError("n cannot exceed 3")
2471 # We know this is a valid EJA, but will double-check
2472 # if the user passes check_axioms=True.
2473 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2475 from mjo
.hurwitz
import OctonionMatrixAlgebra
2476 A
= OctonionMatrixAlgebra(n
, scalars
=field
)
2477 super().__init
__(A
, **kwargs
)
2479 from mjo
.eja
.eja_cache
import octonion_hermitian_eja_coeffs
2480 a
= octonion_hermitian_eja_coeffs(self
)
2482 self
.rational_algebra()._charpoly
_coefficients
.set_cache(a
)
2485 class AlbertEJA(OctonionHermitianEJA
):
2487 The Albert algebra is the algebra of three-by-three Hermitian
2488 matrices whose entries are octonions.
2492 sage: from mjo.eja.eja_algebra import AlbertEJA
2496 sage: AlbertEJA(field=QQ, orthonormalize=False)
2497 Euclidean Jordan algebra of dimension 27 over Rational Field
2498 sage: AlbertEJA() # long time
2499 Euclidean Jordan algebra of dimension 27 over Algebraic Real Field
2502 def __init__(self
, *args
, **kwargs
):
2503 super().__init
__(3, *args
, **kwargs
)
2506 class HadamardEJA(RationalBasisEJA
, ConcreteEJA
):
2508 Return the Euclidean Jordan algebra on `R^n` with the Hadamard
2509 (pointwise real-number multiplication) Jordan product and the
2510 usual inner-product.
2512 This is nothing more than the Cartesian product of ``n`` copies of
2513 the one-dimensional Jordan spin algebra, and is the most common
2514 example of a non-simple Euclidean Jordan algebra.
2518 sage: from mjo.eja.eja_algebra import HadamardEJA
2522 This multiplication table can be verified by hand::
2524 sage: J = HadamardEJA(3)
2525 sage: b0,b1,b2 = J.gens()
2541 We can change the generator prefix::
2543 sage: HadamardEJA(3, prefix='r').gens()
2546 def __init__(self
, n
, field
=AA
, **kwargs
):
2547 MS
= MatrixSpace(field
, n
, 1)
2550 jordan_product
= lambda x
,y
: x
2551 inner_product
= lambda x
,y
: x
2553 def jordan_product(x
,y
):
2554 return MS( xi
*yi
for (xi
,yi
) in zip(x
,y
) )
2556 def inner_product(x
,y
):
2559 # New defaults for keyword arguments. Don't orthonormalize
2560 # because our basis is already orthonormal with respect to our
2561 # inner-product. Don't check the axioms, because we know this
2562 # is a valid EJA... but do double-check if the user passes
2563 # check_axioms=True. Note: we DON'T override the "check_field"
2564 # default here, because the user can pass in a field!
2565 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2566 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2568 column_basis
= tuple( MS(b
) for b
in FreeModule(field
, n
).basis() )
2569 super().__init
__(column_basis
,
2576 self
.rank
.set_cache(n
)
2578 self
.one
.set_cache( self
.sum(self
.gens()) )
2581 def _max_random_instance_dimension():
2583 There's no reason to go higher than five here. That's
2584 enough to get the point across.
2589 def _max_random_instance_size(max_dimension
):
2591 The maximum size (=dimension) of a random HadamardEJA.
2593 return max_dimension
2596 def random_instance(cls
, max_dimension
=None, *args
, **kwargs
):
2598 Return a random instance of this type of algebra.
2600 class_max_d
= cls
._max
_random
_instance
_dimension
()
2601 if (max_dimension
is None or max_dimension
> class_max_d
):
2602 max_dimension
= class_max_d
2603 max_size
= cls
._max
_random
_instance
_size
(max_dimension
)
2604 n
= ZZ
.random_element(max_size
+ 1)
2605 return cls(n
, **kwargs
)
2608 class BilinearFormEJA(RationalBasisEJA
, ConcreteEJA
):
2610 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2611 with the half-trace inner product and jordan product ``x*y =
2612 (<Bx,y>,y_bar>, x0*y_bar + y0*x_bar)`` where `B = 1 \times B22` is
2613 a symmetric positive-definite "bilinear form" matrix. Its
2614 dimension is the size of `B`, and it has rank two in dimensions
2615 larger than two. It reduces to the ``JordanSpinEJA`` when `B` is
2616 the identity matrix of order ``n``.
2618 We insist that the one-by-one upper-left identity block of `B` be
2619 passed in as well so that we can be passed a matrix of size zero
2620 to construct a trivial algebra.
2624 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
2625 ....: JordanSpinEJA)
2629 When no bilinear form is specified, the identity matrix is used,
2630 and the resulting algebra is the Jordan spin algebra::
2632 sage: B = matrix.identity(AA,3)
2633 sage: J0 = BilinearFormEJA(B)
2634 sage: J1 = JordanSpinEJA(3)
2635 sage: J0.multiplication_table() == J0.multiplication_table()
2638 An error is raised if the matrix `B` does not correspond to a
2639 positive-definite bilinear form::
2641 sage: B = matrix.random(QQ,2,3)
2642 sage: J = BilinearFormEJA(B)
2643 Traceback (most recent call last):
2645 ValueError: bilinear form is not positive-definite
2646 sage: B = matrix.zero(QQ,3)
2647 sage: J = BilinearFormEJA(B)
2648 Traceback (most recent call last):
2650 ValueError: bilinear form is not positive-definite
2654 We can create a zero-dimensional algebra::
2656 sage: B = matrix.identity(AA,0)
2657 sage: J = BilinearFormEJA(B)
2661 We can check the multiplication condition given in the Jordan, von
2662 Neumann, and Wigner paper (and also discussed on my "On the
2663 symmetry..." paper). Note that this relies heavily on the standard
2664 choice of basis, as does anything utilizing the bilinear form
2665 matrix. We opt not to orthonormalize the basis, because if we
2666 did, we would have to normalize the `s_{i}` in a similar manner::
2668 sage: n = ZZ.random_element(5)
2669 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2670 sage: B11 = matrix.identity(QQ,1)
2671 sage: B22 = M.transpose()*M
2672 sage: B = block_matrix(2,2,[ [B11,0 ],
2674 sage: J = BilinearFormEJA(B, orthonormalize=False)
2675 sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
2676 sage: V = J.vector_space()
2677 sage: sis = [ J( V([0] + (M.inverse()*ei).list()).column() )
2678 ....: for ei in eis ]
2679 sage: actual = [ sis[i]*sis[j]
2680 ....: for i in range(n-1)
2681 ....: for j in range(n-1) ]
2682 sage: expected = [ J.one() if i == j else J.zero()
2683 ....: for i in range(n-1)
2684 ....: for j in range(n-1) ]
2685 sage: actual == expected
2689 def __init__(self
, B
, field
=AA
, **kwargs
):
2690 # The matrix "B" is supplied by the user in most cases,
2691 # so it makes sense to check whether or not its positive-
2692 # definite unless we are specifically asked not to...
2693 if ("check_axioms" not in kwargs
) or kwargs
["check_axioms"]:
2694 if not B
.is_positive_definite():
2695 raise ValueError("bilinear form is not positive-definite")
2697 # However, all of the other data for this EJA is computed
2698 # by us in manner that guarantees the axioms are
2699 # satisfied. So, again, unless we are specifically asked to
2700 # verify things, we'll skip the rest of the checks.
2701 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2704 MS
= MatrixSpace(field
, n
, 1)
2706 def inner_product(x
,y
):
2707 return (y
.T
*B
*x
)[0,0]
2709 def jordan_product(x
,y
):
2714 z0
= inner_product(y
,x
)
2715 zbar
= y0
*xbar
+ x0
*ybar
2716 return MS([z0
] + zbar
.list())
2718 column_basis
= tuple( MS(b
) for b
in FreeModule(field
, n
).basis() )
2720 # TODO: I haven't actually checked this, but it seems legit.
2725 super().__init
__(column_basis
,
2730 associative
=associative
,
2733 # The rank of this algebra is two, unless we're in a
2734 # one-dimensional ambient space (because the rank is bounded
2735 # by the ambient dimension).
2736 self
.rank
.set_cache(min(n
,2))
2738 self
.one
.set_cache( self
.zero() )
2740 self
.one
.set_cache( self
.monomial(0) )
2743 def _max_random_instance_dimension():
2745 There's no reason to go higher than five here. That's
2746 enough to get the point across.
2751 def _max_random_instance_size(max_dimension
):
2753 The maximum size (=dimension) of a random BilinearFormEJA.
2755 return max_dimension
2758 def random_instance(cls
, max_dimension
=None, *args
, **kwargs
):
2760 Return a random instance of this algebra.
2762 class_max_d
= cls
._max
_random
_instance
_dimension
()
2763 if (max_dimension
is None or max_dimension
> class_max_d
):
2764 max_dimension
= class_max_d
2765 max_size
= cls
._max
_random
_instance
_size
(max_dimension
)
2766 n
= ZZ
.random_element(max_size
+ 1)
2769 B
= matrix
.identity(ZZ
, n
)
2770 return cls(B
, **kwargs
)
2772 B11
= matrix
.identity(ZZ
, 1)
2773 M
= matrix
.random(ZZ
, n
-1)
2774 I
= matrix
.identity(ZZ
, n
-1)
2776 while alpha
.is_zero():
2777 alpha
= ZZ
.random_element().abs()
2779 B22
= M
.transpose()*M
+ alpha
*I
2781 from sage
.matrix
.special
import block_matrix
2782 B
= block_matrix(2,2, [ [B11
, ZZ(0) ],
2785 return cls(B
, **kwargs
)
2788 class JordanSpinEJA(BilinearFormEJA
):
2790 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2791 with the usual inner product and jordan product ``x*y =
2792 (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2797 sage: from mjo.eja.eja_algebra import JordanSpinEJA
2801 This multiplication table can be verified by hand::
2803 sage: J = JordanSpinEJA(4)
2804 sage: b0,b1,b2,b3 = J.gens()
2820 We can change the generator prefix::
2822 sage: JordanSpinEJA(2, prefix='B').gens()
2827 Ensure that we have the usual inner product on `R^n`::
2829 sage: J = JordanSpinEJA.random_instance()
2830 sage: x,y = J.random_elements(2)
2831 sage: actual = x.inner_product(y)
2832 sage: expected = x.to_vector().inner_product(y.to_vector())
2833 sage: actual == expected
2837 def __init__(self
, n
, *args
, **kwargs
):
2838 # This is a special case of the BilinearFormEJA with the
2839 # identity matrix as its bilinear form.
2840 B
= matrix
.identity(ZZ
, n
)
2842 # Don't orthonormalize because our basis is already
2843 # orthonormal with respect to our inner-product.
2844 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2846 # But also don't pass check_field=False here, because the user
2847 # can pass in a field!
2848 super().__init
__(B
, *args
, **kwargs
)
2851 def random_instance(cls
, max_dimension
=None, *args
, **kwargs
):
2853 Return a random instance of this type of algebra.
2855 Needed here to override the implementation for ``BilinearFormEJA``.
2857 class_max_d
= cls
._max
_random
_instance
_dimension
()
2858 if (max_dimension
is None or max_dimension
> class_max_d
):
2859 max_dimension
= class_max_d
2860 max_size
= cls
._max
_random
_instance
_size
(max_dimension
)
2861 n
= ZZ
.random_element(max_size
+ 1)
2862 return cls(n
, **kwargs
)
2865 class TrivialEJA(RationalBasisEJA
, ConcreteEJA
):
2867 The trivial Euclidean Jordan algebra consisting of only a zero element.
2871 sage: from mjo.eja.eja_algebra import TrivialEJA
2875 sage: J = TrivialEJA()
2882 sage: 7*J.one()*12*J.one()
2884 sage: J.one().inner_product(J.one())
2886 sage: J.one().norm()
2888 sage: J.one().subalgebra_generated_by()
2889 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2894 def __init__(self
, field
=AA
, **kwargs
):
2895 jordan_product
= lambda x
,y
: x
2896 inner_product
= lambda x
,y
: field
.zero()
2898 MS
= MatrixSpace(field
,0)
2900 # New defaults for keyword arguments
2901 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2902 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2904 super().__init
__(basis
,
2912 # The rank is zero using my definition, namely the dimension of the
2913 # largest subalgebra generated by any element.
2914 self
.rank
.set_cache(0)
2915 self
.one
.set_cache( self
.zero() )
2918 def random_instance(cls
, max_dimension
=None, *args
, **kwargs
):
2919 # We don't take a "size" argument so the superclass method is
2920 # inappropriate for us. The ``max_dimension`` argument is
2921 # included so that if this method is called generically with a
2922 # ``max_dimension=<whatever>`` argument, we don't try to pass
2923 # it on to the algebra constructor.
2924 return cls(**kwargs
)
2927 class CartesianProductEJA(FiniteDimensionalEJA
):
2929 The external (orthogonal) direct sum of two or more Euclidean
2930 Jordan algebras. Every Euclidean Jordan algebra decomposes into an
2931 orthogonal direct sum of simple Euclidean Jordan algebras which is
2932 then isometric to a Cartesian product, so no generality is lost by
2933 providing only this construction.
2937 sage: from mjo.eja.eja_algebra import (random_eja,
2938 ....: CartesianProductEJA,
2939 ....: ComplexHermitianEJA,
2941 ....: JordanSpinEJA,
2942 ....: RealSymmetricEJA)
2946 The Jordan product is inherited from our factors and implemented by
2947 our CombinatorialFreeModule Cartesian product superclass::
2949 sage: J1 = HadamardEJA(2)
2950 sage: J2 = RealSymmetricEJA(2)
2951 sage: J = cartesian_product([J1,J2])
2952 sage: x,y = J.random_elements(2)
2956 The ability to retrieve the original factors is implemented by our
2957 CombinatorialFreeModule Cartesian product superclass::
2959 sage: J1 = HadamardEJA(2, field=QQ)
2960 sage: J2 = JordanSpinEJA(3, field=QQ)
2961 sage: J = cartesian_product([J1,J2])
2962 sage: J.cartesian_factors()
2963 (Euclidean Jordan algebra of dimension 2 over Rational Field,
2964 Euclidean Jordan algebra of dimension 3 over Rational Field)
2966 You can provide more than two factors::
2968 sage: J1 = HadamardEJA(2)
2969 sage: J2 = JordanSpinEJA(3)
2970 sage: J3 = RealSymmetricEJA(3)
2971 sage: cartesian_product([J1,J2,J3])
2972 Euclidean Jordan algebra of dimension 2 over Algebraic Real
2973 Field (+) Euclidean Jordan algebra of dimension 3 over Algebraic
2974 Real Field (+) Euclidean Jordan algebra of dimension 6 over
2975 Algebraic Real Field
2977 Rank is additive on a Cartesian product::
2979 sage: J1 = HadamardEJA(1)
2980 sage: J2 = RealSymmetricEJA(2)
2981 sage: J = cartesian_product([J1,J2])
2982 sage: J1.rank.clear_cache()
2983 sage: J2.rank.clear_cache()
2984 sage: J.rank.clear_cache()
2987 sage: J.rank() == J1.rank() + J2.rank()
2990 The same rank computation works over the rationals, with whatever
2993 sage: J1 = HadamardEJA(1, field=QQ, orthonormalize=False)
2994 sage: J2 = RealSymmetricEJA(2, field=QQ, orthonormalize=False)
2995 sage: J = cartesian_product([J1,J2])
2996 sage: J1.rank.clear_cache()
2997 sage: J2.rank.clear_cache()
2998 sage: J.rank.clear_cache()
3001 sage: J.rank() == J1.rank() + J2.rank()
3004 The product algebra will be associative if and only if all of its
3005 components are associative::
3007 sage: J1 = HadamardEJA(2)
3008 sage: J1.is_associative()
3010 sage: J2 = HadamardEJA(3)
3011 sage: J2.is_associative()
3013 sage: J3 = RealSymmetricEJA(3)
3014 sage: J3.is_associative()
3016 sage: CP1 = cartesian_product([J1,J2])
3017 sage: CP1.is_associative()
3019 sage: CP2 = cartesian_product([J1,J3])
3020 sage: CP2.is_associative()
3023 Cartesian products of Cartesian products work::
3025 sage: J1 = JordanSpinEJA(1)
3026 sage: J2 = JordanSpinEJA(1)
3027 sage: J3 = JordanSpinEJA(1)
3028 sage: J = cartesian_product([J1,cartesian_product([J2,J3])])
3029 sage: J.multiplication_table()
3030 +----++----+----+----+
3031 | * || b0 | b1 | b2 |
3032 +====++====+====+====+
3033 | b0 || b0 | 0 | 0 |
3034 +----++----+----+----+
3035 | b1 || 0 | b1 | 0 |
3036 +----++----+----+----+
3037 | b2 || 0 | 0 | b2 |
3038 +----++----+----+----+
3039 sage: HadamardEJA(3).multiplication_table()
3040 +----++----+----+----+
3041 | * || b0 | b1 | b2 |
3042 +====++====+====+====+
3043 | b0 || b0 | 0 | 0 |
3044 +----++----+----+----+
3045 | b1 || 0 | b1 | 0 |
3046 +----++----+----+----+
3047 | b2 || 0 | 0 | b2 |
3048 +----++----+----+----+
3050 The "matrix space" of a Cartesian product always consists of
3051 ordered pairs (or triples, or...) whose components are the
3052 matrix spaces of its factors::
3054 sage: J1 = HadamardEJA(2)
3055 sage: J2 = ComplexHermitianEJA(2)
3056 sage: J = cartesian_product([J1,J2])
3057 sage: J.matrix_space()
3058 The Cartesian product of (Full MatrixSpace of 2 by 1 dense
3059 matrices over Algebraic Real Field, Module of 2 by 2 matrices
3060 with entries in Algebraic Field over the scalar ring Algebraic
3062 sage: J.one().to_matrix()[0]
3065 sage: J.one().to_matrix()[1]
3074 All factors must share the same base field::
3076 sage: J1 = HadamardEJA(2, field=QQ)
3077 sage: J2 = RealSymmetricEJA(2)
3078 sage: CartesianProductEJA((J1,J2))
3079 Traceback (most recent call last):
3081 ValueError: all factors must share the same base field
3083 The cached unit element is the same one that would be computed::
3085 sage: J1 = random_eja() # long time
3086 sage: J2 = random_eja() # long time
3087 sage: J = cartesian_product([J1,J2]) # long time
3088 sage: actual = J.one() # long time
3089 sage: J.one.clear_cache() # long time
3090 sage: expected = J.one() # long time
3091 sage: actual == expected # long time
3094 Element
= CartesianProductEJAElement
3095 def __init__(self
, factors
, **kwargs
):
3100 self
._sets
= factors
3102 field
= factors
[0].base_ring()
3103 if not all( J
.base_ring() == field
for J
in factors
):
3104 raise ValueError("all factors must share the same base field")
3106 # Figure out the category to use.
3107 associative
= all( f
.is_associative() for f
in factors
)
3108 category
= EuclideanJordanAlgebras(field
)
3109 if associative
: category
= category
.Associative()
3110 category
= category
.join([category
, category
.CartesianProducts()])
3112 # Compute my matrix space. We don't simply use the
3113 # ``cartesian_product()`` functor here because it acts
3114 # differently on SageMath MatrixSpaces and our custom
3115 # MatrixAlgebras, which are CombinatorialFreeModules. We
3116 # always want the result to be represented (and indexed) as an
3117 # ordered tuple. This category isn't perfect, but is good
3118 # enough for what we need to do.
3119 MS_cat
= MagmaticAlgebras(field
).FiniteDimensional().WithBasis()
3120 MS_cat
= MS_cat
.Unital().CartesianProducts()
3121 MS_factors
= tuple( J
.matrix_space() for J
in factors
)
3122 from sage
.sets
.cartesian_product
import CartesianProduct
3123 self
._matrix
_space
= CartesianProduct(MS_factors
, MS_cat
)
3125 self
._matrix
_basis
= []
3126 zero
= self
._matrix
_space
.zero()
3128 for b
in factors
[i
].matrix_basis():
3131 self
._matrix
_basis
.append(z
)
3133 self
._matrix
_basis
= tuple( self
._matrix
_space
(b
)
3134 for b
in self
._matrix
_basis
)
3135 n
= len(self
._matrix
_basis
)
3137 # We already have what we need for the super-superclass constructor.
3138 CombinatorialFreeModule
.__init
__(self
,
3145 # Now create the vector space for the algebra, which will have
3146 # its own set of non-ambient coordinates (in terms of the
3148 degree
= sum( f
._matrix
_span
.ambient_vector_space().degree()
3150 V
= VectorSpace(field
, degree
)
3151 vector_basis
= tuple( V(_all2list(b
)) for b
in self
._matrix
_basis
)
3153 # Save the span of our matrix basis (when written out as long
3154 # vectors) because otherwise we'll have to reconstruct it
3155 # every time we want to coerce a matrix into the algebra.
3156 self
._matrix
_span
= V
.span_of_basis( vector_basis
, check
=False)
3158 # Since we don't (re)orthonormalize the basis, the FDEJA
3159 # constructor is going to set self._deortho_matrix to the
3160 # identity matrix. Here we set it to the correct value using
3161 # the deortho matrices from our factors.
3162 self
._deortho
_matrix
= matrix
.block_diagonal(
3163 [J
._deortho
_matrix
for J
in factors
]
3166 self
._inner
_product
_matrix
= matrix
.block_diagonal(
3167 [J
._inner
_product
_matrix
for J
in factors
]
3169 self
._inner
_product
_matrix
._cache
= {'hermitian': True}
3170 self
._inner
_product
_matrix
.set_immutable()
3172 # Building the multiplication table is a bit more tricky
3173 # because we have to embed the entries of the factors'
3174 # multiplication tables into the product EJA.
3176 self
._multiplication
_table
= [ [zed
for j
in range(i
+1)]
3179 # Keep track of an offset that tallies the dimensions of all
3180 # previous factors. If the second factor is dim=2 and if the
3181 # first one is dim=3, then we want to skip the first 3x3 block
3182 # when copying the multiplication table for the second factor.
3185 phi_f
= self
.cartesian_embedding(f
)
3186 factor_dim
= factors
[f
].dimension()
3187 for i
in range(factor_dim
):
3188 for j
in range(i
+1):
3189 f_ij
= factors
[f
]._multiplication
_table
[i
][j
]
3191 self
._multiplication
_table
[offset
+i
][offset
+j
] = e
3192 offset
+= factor_dim
3194 self
.rank
.set_cache(sum(J
.rank() for J
in factors
))
3195 ones
= tuple(J
.one().to_matrix() for J
in factors
)
3196 self
.one
.set_cache(self(ones
))
3198 def _sets_keys(self
):
3203 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
3204 ....: RealSymmetricEJA)
3208 The superclass uses ``_sets_keys()`` to implement its
3209 ``cartesian_factors()`` method::
3211 sage: J1 = RealSymmetricEJA(2,
3213 ....: orthonormalize=False,
3215 sage: J2 = ComplexHermitianEJA(2,field=QQ,orthonormalize=False)
3216 sage: J = cartesian_product([J1,J2])
3217 sage: x = sum(i*J.gens()[i] for i in range(len(J.gens())))
3218 sage: x.cartesian_factors()
3219 (a1 + 2*a2, 3*b0 + 4*b1 + 5*b2 + 6*b3)
3222 # Copy/pasted from CombinatorialFreeModule_CartesianProduct,
3223 # but returning a tuple instead of a list.
3224 return tuple(range(len(self
.cartesian_factors())))
3226 def cartesian_factors(self
):
3227 # Copy/pasted from CombinatorialFreeModule_CartesianProduct.
3230 def cartesian_factor(self
, i
):
3232 Return the ``i``th factor of this algebra.
3234 return self
._sets
[i
]
3237 # Copy/pasted from CombinatorialFreeModule_CartesianProduct.
3238 from sage
.categories
.cartesian_product
import cartesian_product
3239 return cartesian_product
.symbol
.join("%s" % factor
3240 for factor
in self
._sets
)
3244 def cartesian_projection(self
, i
):
3248 sage: from mjo.eja.eja_algebra import (random_eja,
3249 ....: JordanSpinEJA,
3251 ....: RealSymmetricEJA,
3252 ....: ComplexHermitianEJA)
3256 The projection morphisms are Euclidean Jordan algebra
3259 sage: J1 = HadamardEJA(2)
3260 sage: J2 = RealSymmetricEJA(2)
3261 sage: J = cartesian_product([J1,J2])
3262 sage: J.cartesian_projection(0)
3263 Linear operator between finite-dimensional Euclidean Jordan
3264 algebras represented by the matrix:
3267 Domain: Euclidean Jordan algebra of dimension 2 over Algebraic
3268 Real Field (+) Euclidean Jordan algebra of dimension 3 over
3269 Algebraic Real Field
3270 Codomain: Euclidean Jordan algebra of dimension 2 over Algebraic
3272 sage: J.cartesian_projection(1)
3273 Linear operator between finite-dimensional Euclidean Jordan
3274 algebras represented by the matrix:
3278 Domain: Euclidean Jordan algebra of dimension 2 over Algebraic
3279 Real Field (+) Euclidean Jordan algebra of dimension 3 over
3280 Algebraic Real Field
3281 Codomain: Euclidean Jordan algebra of dimension 3 over Algebraic
3284 The projections work the way you'd expect on the vector
3285 representation of an element::
3287 sage: J1 = JordanSpinEJA(2)
3288 sage: J2 = ComplexHermitianEJA(2)
3289 sage: J = cartesian_product([J1,J2])
3290 sage: pi_left = J.cartesian_projection(0)
3291 sage: pi_right = J.cartesian_projection(1)
3292 sage: pi_left(J.one()).to_vector()
3294 sage: pi_right(J.one()).to_vector()
3296 sage: J.one().to_vector()
3301 The answer never changes::
3303 sage: J1 = random_eja()
3304 sage: J2 = random_eja()
3305 sage: J = cartesian_product([J1,J2])
3306 sage: P0 = J.cartesian_projection(0)
3307 sage: P1 = J.cartesian_projection(0)
3312 offset
= sum( self
.cartesian_factor(k
).dimension()
3314 Ji
= self
.cartesian_factor(i
)
3315 Pi
= self
._module
_morphism
(lambda j
: Ji
.monomial(j
- offset
),
3318 return FiniteDimensionalEJAOperator(self
,Ji
,Pi
.matrix())
3321 def cartesian_embedding(self
, i
):
3325 sage: from mjo.eja.eja_algebra import (random_eja,
3326 ....: JordanSpinEJA,
3328 ....: RealSymmetricEJA)
3332 The embedding morphisms are Euclidean Jordan algebra
3335 sage: J1 = HadamardEJA(2)
3336 sage: J2 = RealSymmetricEJA(2)
3337 sage: J = cartesian_product([J1,J2])
3338 sage: J.cartesian_embedding(0)
3339 Linear operator between finite-dimensional Euclidean Jordan
3340 algebras represented by the matrix:
3346 Domain: Euclidean Jordan algebra of dimension 2 over
3347 Algebraic Real Field
3348 Codomain: Euclidean Jordan algebra of dimension 2 over
3349 Algebraic Real Field (+) Euclidean Jordan algebra of
3350 dimension 3 over Algebraic Real Field
3351 sage: J.cartesian_embedding(1)
3352 Linear operator between finite-dimensional Euclidean Jordan
3353 algebras represented by the matrix:
3359 Domain: Euclidean Jordan algebra of dimension 3 over
3360 Algebraic Real Field
3361 Codomain: Euclidean Jordan algebra of dimension 2 over
3362 Algebraic Real Field (+) Euclidean Jordan algebra of
3363 dimension 3 over Algebraic Real Field
3365 The embeddings work the way you'd expect on the vector
3366 representation of an element::
3368 sage: J1 = JordanSpinEJA(3)
3369 sage: J2 = RealSymmetricEJA(2)
3370 sage: J = cartesian_product([J1,J2])
3371 sage: iota_left = J.cartesian_embedding(0)
3372 sage: iota_right = J.cartesian_embedding(1)
3373 sage: iota_left(J1.zero()) == J.zero()
3375 sage: iota_right(J2.zero()) == J.zero()
3377 sage: J1.one().to_vector()
3379 sage: iota_left(J1.one()).to_vector()
3381 sage: J2.one().to_vector()
3383 sage: iota_right(J2.one()).to_vector()
3385 sage: J.one().to_vector()
3390 The answer never changes::
3392 sage: J1 = random_eja()
3393 sage: J2 = random_eja()
3394 sage: J = cartesian_product([J1,J2])
3395 sage: E0 = J.cartesian_embedding(0)
3396 sage: E1 = J.cartesian_embedding(0)
3400 Composing a projection with the corresponding inclusion should
3401 produce the identity map, and mismatching them should produce
3404 sage: J1 = random_eja()
3405 sage: J2 = random_eja()
3406 sage: J = cartesian_product([J1,J2])
3407 sage: iota_left = J.cartesian_embedding(0)
3408 sage: iota_right = J.cartesian_embedding(1)
3409 sage: pi_left = J.cartesian_projection(0)
3410 sage: pi_right = J.cartesian_projection(1)
3411 sage: pi_left*iota_left == J1.one().operator()
3413 sage: pi_right*iota_right == J2.one().operator()
3415 sage: (pi_left*iota_right).is_zero()
3417 sage: (pi_right*iota_left).is_zero()
3421 offset
= sum( self
.cartesian_factor(k
).dimension()
3423 Ji
= self
.cartesian_factor(i
)
3424 Ei
= Ji
._module
_morphism
(lambda j
: self
.monomial(j
+ offset
),
3426 return FiniteDimensionalEJAOperator(Ji
,self
,Ei
.matrix())
3430 FiniteDimensionalEJA
.CartesianProduct
= CartesianProductEJA
3432 class RationalBasisCartesianProductEJA(CartesianProductEJA
,
3435 A separate class for products of algebras for which we know a
3440 sage: from mjo.eja.eja_algebra import (FiniteDimensionalEJA,
3442 ....: JordanSpinEJA,
3443 ....: RealSymmetricEJA)
3447 This gives us fast characteristic polynomial computations in
3448 product algebras, too::
3451 sage: J1 = JordanSpinEJA(2)
3452 sage: J2 = RealSymmetricEJA(3)
3453 sage: J = cartesian_product([J1,J2])
3454 sage: J.characteristic_polynomial_of().degree()
3461 The ``cartesian_product()`` function only uses the first factor to
3462 decide where the result will live; thus we have to be careful to
3463 check that all factors do indeed have a ``rational_algebra()`` method
3464 before we construct an algebra that claims to have a rational basis::
3466 sage: J1 = HadamardEJA(2)
3467 sage: jp = lambda X,Y: X*Y
3468 sage: ip = lambda X,Y: X[0,0]*Y[0,0]
3469 sage: b1 = matrix(QQ, [[1]])
3470 sage: J2 = FiniteDimensionalEJA((b1,), jp, ip)
3471 sage: cartesian_product([J2,J1]) # factor one not RationalBasisEJA
3472 Euclidean Jordan algebra of dimension 1 over Algebraic Real
3473 Field (+) Euclidean Jordan algebra of dimension 2 over Algebraic
3475 sage: cartesian_product([J1,J2]) # factor one is RationalBasisEJA
3476 Traceback (most recent call last):
3478 ValueError: factor not a RationalBasisEJA
3481 def __init__(self
, algebras
, **kwargs
):
3482 if not all( hasattr(r
, "rational_algebra") for r
in algebras
):
3483 raise ValueError("factor not a RationalBasisEJA")
3485 CartesianProductEJA
.__init
__(self
, algebras
, **kwargs
)
3488 def rational_algebra(self
):
3489 if self
.base_ring() is QQ
:
3492 return cartesian_product([
3493 r
.rational_algebra() for r
in self
.cartesian_factors()
3497 RationalBasisEJA
.CartesianProduct
= RationalBasisCartesianProductEJA
3499 def random_eja(max_dimension
=None, *args
, **kwargs
):
3504 sage: from mjo.eja.eja_algebra import random_eja
3508 sage: n = ZZ.random_element(1,5)
3509 sage: J = random_eja(max_dimension=n, field=QQ, orthonormalize=False)
3510 sage: J.dimension() <= n
3514 # Use the ConcreteEJA default as the total upper bound (regardless
3515 # of any whether or not any individual factors set a lower limit).
3516 if max_dimension
is None:
3517 max_dimension
= ConcreteEJA
._max
_random
_instance
_dimension
()
3518 J1
= ConcreteEJA
.random_instance(max_dimension
, *args
, **kwargs
)
3521 # Roll the dice to see if we attempt a Cartesian product.
3522 dice_roll
= ZZ
.random_element(len(ConcreteEJA
.__subclasses
__()) + 1)
3523 new_max_dimension
= max_dimension
- J1
.dimension()
3524 if new_max_dimension
== 0 or dice_roll
!= 0:
3525 # If it's already as big as we're willing to tolerate, just
3526 # return it and don't worry about Cartesian products.
3529 # Use random_eja() again so we can get more than two factors
3530 # if the sub-call also Decides on a cartesian product.
3531 J2
= random_eja(new_max_dimension
, *args
, **kwargs
)
3532 return cartesian_product([J1
,J2
])
3535 class ComplexSkewSymmetricEJA(RationalBasisEJA
, ConcreteEJA
):
3537 The skew-symmetric EJA of order `2n` described in Faraut and
3538 Koranyi's Exercise III.1.b. It has dimension `2n^2 - n`.
3540 It is (not obviously) isomorphic to the QuaternionHermitianEJA of
3541 order `n`, as can be inferred by comparing rank/dimension or
3542 explicitly from their "characteristic polynomial of" functions,
3543 which just so happen to align nicely.
3547 sage: from mjo.eja.eja_algebra import (ComplexSkewSymmetricEJA,
3548 ....: QuaternionHermitianEJA)
3549 sage: from mjo.eja.eja_operator import FiniteDimensionalEJAOperator
3553 This EJA is isomorphic to the quaternions::
3555 sage: J = ComplexSkewSymmetricEJA(2, field=QQ, orthonormalize=False)
3556 sage: K = QuaternionHermitianEJA(2, field=QQ, orthonormalize=False)
3557 sage: jordan_isom_matrix = matrix.diagonal(QQ,[-1,1,1,1,1,-1])
3558 sage: phi = FiniteDimensionalEJAOperator(J,K,jordan_isom_matrix)
3559 sage: all( phi(x*y) == phi(x)*phi(y)
3560 ....: for x in J.gens()
3561 ....: for y in J.gens() )
3563 sage: x,y = J.random_elements(2)
3564 sage: phi(x*y) == phi(x)*phi(y)
3569 Random elements should satisfy the same conditions that the basis
3572 sage: K = ComplexSkewSymmetricEJA.random_instance(field=QQ,
3573 ....: orthonormalize=False)
3574 sage: x,y = K.random_elements(2)
3576 sage: x = x.to_matrix()
3577 sage: y = y.to_matrix()
3578 sage: z = z.to_matrix()
3579 sage: all( e.is_skew_symmetric() for e in (x,y,z) )
3581 sage: J = -K.one().to_matrix()
3582 sage: all( e*J == J*e.conjugate() for e in (x,y,z) )
3585 The power law in Faraut & Koranyi's II.7.a is satisfied.
3586 We're in a subalgebra of theirs, but powers are still
3589 sage: K = ComplexSkewSymmetricEJA.random_instance(field=QQ,
3590 ....: orthonormalize=False)
3591 sage: x = K.random_element()
3592 sage: k = ZZ.random_element(5)
3594 sage: J = -K.one().to_matrix()
3595 sage: expected = K(-J*(J*x.to_matrix())^k)
3596 sage: actual == expected
3601 def _max_random_instance_size(max_dimension
):
3602 # Obtained by solving d = 2n^2 - n, which comes from noticing
3603 # that, in 2x2 block form, any element of this algebra has a
3604 # free skew-symmetric top-left block, a Hermitian top-right
3605 # block, and two bottom blocks that are determined by the top.
3606 # The ZZ-int-ZZ thing is just "floor."
3607 return ZZ(int(ZZ(8*max_dimension
+ 1).sqrt()/4 + 1/4))
3610 def random_instance(cls
, max_dimension
=None, *args
, **kwargs
):
3612 Return a random instance of this type of algebra.
3614 class_max_d
= cls
._max
_random
_instance
_dimension
()
3615 if (max_dimension
is None or max_dimension
> class_max_d
):
3616 max_dimension
= class_max_d
3617 max_size
= cls
._max
_random
_instance
_size
(max_dimension
)
3618 n
= ZZ
.random_element(max_size
+ 1)
3619 return cls(n
, **kwargs
)
3622 def _denormalized_basis(A
):
3626 sage: from mjo.hurwitz import ComplexMatrixAlgebra
3627 sage: from mjo.eja.eja_algebra import ComplexSkewSymmetricEJA
3631 The basis elements are all skew-Hermitian::
3633 sage: d_max = ComplexSkewSymmetricEJA._max_random_instance_dimension()
3634 sage: n_max = ComplexSkewSymmetricEJA._max_random_instance_size(d_max)
3635 sage: n = ZZ.random_element(n_max + 1)
3636 sage: A = ComplexMatrixAlgebra(2*n, scalars=QQ)
3637 sage: B = ComplexSkewSymmetricEJA._denormalized_basis(A)
3638 sage: all( M.is_skew_symmetric() for M in B)
3641 The basis elements ``b`` all satisfy ``b*J == J*b.conjugate()``,
3642 as in the definition of the algebra::
3644 sage: d_max = ComplexSkewSymmetricEJA._max_random_instance_dimension()
3645 sage: n_max = ComplexSkewSymmetricEJA._max_random_instance_size(d_max)
3646 sage: n = ZZ.random_element(n_max + 1)
3647 sage: A = ComplexMatrixAlgebra(2*n, scalars=QQ)
3648 sage: I_n = matrix.identity(ZZ, n)
3649 sage: J = matrix.block(ZZ, 2, 2, (0, I_n, -I_n, 0), subdivide=False)
3650 sage: J = A.from_list(J.rows())
3651 sage: B = ComplexSkewSymmetricEJA._denormalized_basis(A)
3652 sage: all( b*J == J*b.conjugate() for b in B )
3656 es
= A
.entry_algebra_gens()
3657 gen
= lambda A
,m
: A
.monomial(m
)
3661 # The size of the blocks. We're going to treat these thing as
3662 # 2x2 block matrices,
3665 # [ -x2-conj x1-conj ]
3667 # where x1 is skew-symmetric and x2 is Hermitian.
3671 # We only loop through the top half of the matrix, because the
3672 # bottom can be constructed from the top.
3674 # First do the top-left block, which is skew-symmetric.
3675 # We can compute the bottom-right block in the process.
3676 for j
in range(i
+1):
3678 # Skew-symmetry implies zeros for (i == j).
3680 # Top-left block's entry.
3681 E_ij
= gen(A
, (i
,j
,e
))
3682 E_ij
-= gen(A
, (j
,i
,e
))
3684 # Bottom-right block's entry.
3685 F_ij
= gen(A
, (i
+m
,j
+m
,e
)).conjugate()
3686 F_ij
-= gen(A
, (j
+m
,i
+m
,e
)).conjugate()
3688 basis
.append(E_ij
+ F_ij
)
3690 # Now do the top-right block, which is Hermitian, and compute
3691 # the bottom-left block along the way.
3692 for j
in range(m
,i
+m
+1):
3694 # Hermitian matrices have real diagonal entries.
3695 # Top-right block's entry.
3696 E_ii
= gen(A
, (i
,j
,es
[0]))
3698 # Bottom-left block's entry. Don't conjugate
3700 E_ii
-= gen(A
, (i
+m
,j
-m
,es
[0]))
3704 # Top-right block's entry. BEWARE! We're not
3705 # reflecting across the main diagonal as in
3706 # (i,j)~(j,i). We're only reflecting across
3707 # the diagonal for the top-right block.
3708 E_ij
= gen(A
, (i
,j
,e
))
3710 # Shift it back to non-offset coords, transpose,
3711 # conjugate, and put it back:
3713 # (i,j) -> (i,j-m) -> (j-m, i) -> (j-m, i+m)
3714 E_ij
+= gen(A
, (j
-m
,i
+m
,e
)).conjugate()
3716 # Bottom-left's block's below-diagonal entry.
3717 # Just shift the top-right coords down m and
3719 F_ij
= -gen(A
, (i
+m
,j
-m
,e
)).conjugate()
3720 F_ij
+= -gen(A
, (j
,i
,e
)) # double-conjugate cancels
3722 basis
.append(E_ij
+ F_ij
)
3724 return tuple( basis
)
3728 def _J_matrix(matrix_space
):
3729 n
= matrix_space
.nrows() // 2
3730 F
= matrix_space
.base_ring()
3731 I_n
= matrix
.identity(F
, n
)
3732 J
= matrix
.block(F
, 2, 2, (0, I_n
, -I_n
, 0), subdivide
=False)
3733 return matrix_space
.from_list(J
.rows())
3736 return ComplexSkewSymmetricEJA
._J
_matrix
(self
.matrix_space())
3738 def __init__(self
, n
, field
=AA
, **kwargs
):
3739 # New code; always check the axioms.
3740 #if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
3742 from mjo
.hurwitz
import ComplexMatrixAlgebra
3743 A
= ComplexMatrixAlgebra(2*n
, scalars
=field
)
3744 J
= ComplexSkewSymmetricEJA
._J
_matrix
(A
)
3746 def jordan_product(X
,Y
):
3747 return (X
*J
*Y
+ Y
*J
*X
)/2
3749 def inner_product(X
,Y
):
3750 return (X
.conjugate_transpose()*Y
).trace().real()
3752 super().__init
__(self
._denormalized
_basis
(A
),
3759 # This algebra is conjectured (by me) to be isomorphic to
3760 # the quaternion Hermitian EJA of size n, and the rank
3761 # would follow from that.
3762 #self.rank.set_cache(n)
3763 self
.one
.set_cache( self(-J
) )