2 Representations and constructions for Euclidean Jordan algebras.
4 A Euclidean Jordan algebra is a Jordan algebra that has some
7 1. It is finite-dimensional.
8 2. Its scalar field is the real numbers.
9 3a. An inner product is defined on it, and...
10 3b. That inner product is compatible with the Jordan product
11 in the sense that `<x*y,z> = <y,x*z>` for all elements
12 `x,y,z` in the algebra.
14 Every Euclidean Jordan algebra is formally-real: for any two elements
15 `x` and `y` in the algebra, `x^{2} + y^{2} = 0` implies that `x = y =
16 0`. Conversely, every finite-dimensional formally-real Jordan algebra
17 can be made into a Euclidean Jordan algebra with an appropriate choice
20 Formally-real Jordan algebras were originally studied as a framework
21 for quantum mechanics. Today, Euclidean Jordan algebras are crucial in
22 symmetric cone optimization, since every symmetric cone arises as the
23 cone of squares in some Euclidean Jordan algebra.
25 It is known that every Euclidean Jordan algebra decomposes into an
26 orthogonal direct sum (essentially, a Cartesian product) of simple
27 algebras, and that moreover, up to Jordan-algebra isomorphism, there
28 are only five families of simple algebras. We provide constructions
29 for these simple algebras:
31 * :class:`BilinearFormEJA`
32 * :class:`RealSymmetricEJA`
33 * :class:`ComplexHermitianEJA`
34 * :class:`QuaternionHermitianEJA`
35 * :class:`OctonionHermitianEJA`
37 In addition to these, we provide two other example constructions,
39 * :class:`JordanSpinEJA`
40 * :class:`HadamardEJA`
44 The Jordan spin algebra is a bilinear form algebra where the bilinear
45 form is the identity. The Hadamard EJA is simply a Cartesian product
46 of one-dimensional spin algebras. The Albert EJA is simply a special
47 case of the :class:`OctonionHermitianEJA` where the matrices are
48 three-by-three and the resulting space has dimension 27. And
49 last/least, the trivial EJA is exactly what you think it is; it could
50 also be obtained by constructing a dimension-zero instance of any of
51 the other algebras. Cartesian products of these are also supported
52 using the usual ``cartesian_product()`` function; as a result, we
53 support (up to isomorphism) all Euclidean Jordan algebras.
57 sage: from mjo.eja.eja_algebra import random_eja
62 Euclidean Jordan algebra of dimension...
65 from sage
.algebras
.quatalg
.quaternion_algebra
import QuaternionAlgebra
66 from sage
.categories
.magmatic_algebras
import MagmaticAlgebras
67 from sage
.categories
.sets_cat
import cartesian_product
68 from sage
.combinat
.free_module
import CombinatorialFreeModule
69 from sage
.matrix
.constructor
import matrix
70 from sage
.matrix
.matrix_space
import MatrixSpace
71 from sage
.misc
.cachefunc
import cached_method
72 from sage
.misc
.table
import table
73 from sage
.modules
.free_module
import FreeModule
, VectorSpace
74 from sage
.rings
.all
import (ZZ
, QQ
, AA
, QQbar
, RR
, RLF
, CLF
,
77 from mjo
.eja
.eja_element
import FiniteDimensionalEJAElement
78 from mjo
.eja
.eja_operator
import FiniteDimensionalEJAOperator
79 from mjo
.eja
.eja_utils
import _all2list
, _mat2vec
81 class FiniteDimensionalEJA(CombinatorialFreeModule
):
83 A finite-dimensional Euclidean Jordan algebra.
87 - ``basis`` -- a tuple; a tuple of basis elements in "matrix
88 form," which must be the same form as the arguments to
89 ``jordan_product`` and ``inner_product``. In reality, "matrix
90 form" can be either vectors, matrices, or a Cartesian product
91 (ordered tuple) of vectors or matrices. All of these would
92 ideally be vector spaces in sage with no special-casing
93 needed; but in reality we turn vectors into column-matrices
94 and Cartesian products `(a,b)` into column matrices
95 `(a,b)^{T}` after converting `a` and `b` themselves.
97 - ``jordan_product`` -- a function; afunction of two ``basis``
98 elements (in matrix form) that returns their jordan product,
99 also in matrix form; this will be applied to ``basis`` to
100 compute a multiplication table for the algebra.
102 - ``inner_product`` -- a function; a function of two ``basis``
103 elements (in matrix form) that returns their inner
104 product. This will be applied to ``basis`` to compute an
105 inner-product table (basically a matrix) for this algebra.
107 - ``field`` -- a subfield of the reals (default: ``AA``); the scalar
108 field for the algebra.
110 - ``orthonormalize`` -- boolean (default: ``True``); whether or
111 not to orthonormalize the basis. Doing so is expensive and
112 generally rules out using the rationals as your ``field``, but
113 is required for spectral decompositions.
117 sage: from mjo.eja.eja_algebra import random_eja
121 We should compute that an element subalgebra is associative even
122 if we circumvent the element method::
124 sage: set_random_seed()
125 sage: J = random_eja(field=QQ,orthonormalize=False)
126 sage: x = J.random_element()
127 sage: A = x.subalgebra_generated_by(orthonormalize=False)
128 sage: basis = tuple(b.superalgebra_element() for b in A.basis())
129 sage: J.subalgebra(basis, orthonormalize=False).is_associative()
133 Element
= FiniteDimensionalEJAElement
142 cartesian_product
=False,
150 if not field
.is_subring(RR
):
151 # Note: this does return true for the real algebraic
152 # field, the rationals, and any quadratic field where
153 # we've specified a real embedding.
154 raise ValueError("scalar field is not real")
157 # Check commutativity of the Jordan and inner-products.
158 # This has to be done before we build the multiplication
159 # and inner-product tables/matrices, because we take
160 # advantage of symmetry in the process.
161 if not all( jordan_product(bi
,bj
) == jordan_product(bj
,bi
)
164 raise ValueError("Jordan product is not commutative")
166 if not all( inner_product(bi
,bj
) == inner_product(bj
,bi
)
169 raise ValueError("inner-product is not commutative")
172 category
= MagmaticAlgebras(field
).FiniteDimensional()
173 category
= category
.WithBasis().Unital().Commutative()
176 # All zero- and one-dimensional algebras are just the real
177 # numbers with (some positive multiples of) the usual
178 # multiplication as its Jordan and inner-product.
180 if associative
is None:
181 # We should figure it out. As with check_axioms, we have to do
182 # this without the help of the _jordan_product_is_associative()
183 # method because we need to know the category before we
184 # initialize the algebra.
185 associative
= all( jordan_product(jordan_product(bi
,bj
),bk
)
187 jordan_product(bi
,jordan_product(bj
,bk
))
193 # Element subalgebras can take advantage of this.
194 category
= category
.Associative()
195 if cartesian_product
:
196 # Use join() here because otherwise we only get the
197 # "Cartesian product of..." and not the things themselves.
198 category
= category
.join([category
,
199 category
.CartesianProducts()])
201 # Call the superclass constructor so that we can use its from_vector()
202 # method to build our multiplication table.
203 CombinatorialFreeModule
.__init
__(self
,
210 # Now comes all of the hard work. We'll be constructing an
211 # ambient vector space V that our (vectorized) basis lives in,
212 # as well as a subspace W of V spanned by those (vectorized)
213 # basis elements. The W-coordinates are the coefficients that
214 # we see in things like x = 1*b1 + 2*b2.
219 degree
= len(_all2list(basis
[0]))
221 # Build an ambient space that fits our matrix basis when
222 # written out as "long vectors."
223 V
= VectorSpace(field
, degree
)
225 # The matrix that will hole the orthonormal -> unorthonormal
226 # coordinate transformation.
227 self
._deortho
_matrix
= None
230 # Save a copy of the un-orthonormalized basis for later.
231 # Convert it to ambient V (vector) coordinates while we're
232 # at it, because we'd have to do it later anyway.
233 deortho_vector_basis
= tuple( V(_all2list(b
)) for b
in basis
)
235 from mjo
.eja
.eja_utils
import gram_schmidt
236 basis
= tuple(gram_schmidt(basis
, inner_product
))
238 # Save the (possibly orthonormalized) matrix basis for
240 self
._matrix
_basis
= basis
242 # Now create the vector space for the algebra, which will have
243 # its own set of non-ambient coordinates (in terms of the
245 vector_basis
= tuple( V(_all2list(b
)) for b
in basis
)
246 W
= V
.span_of_basis( vector_basis
, check
=check_axioms
)
249 # Now "W" is the vector space of our algebra coordinates. The
250 # variables "X1", "X2",... refer to the entries of vectors in
251 # W. Thus to convert back and forth between the orthonormal
252 # coordinates and the given ones, we need to stick the original
254 U
= V
.span_of_basis( deortho_vector_basis
, check
=check_axioms
)
255 self
._deortho
_matrix
= matrix( U
.coordinate_vector(q
)
256 for q
in vector_basis
)
259 # Now we actually compute the multiplication and inner-product
260 # tables/matrices using the possibly-orthonormalized basis.
261 self
._inner
_product
_matrix
= matrix
.identity(field
, n
)
262 self
._multiplication
_table
= [ [0 for j
in range(i
+1)]
265 # Note: the Jordan and inner-products are defined in terms
266 # of the ambient basis. It's important that their arguments
267 # are in ambient coordinates as well.
270 # ortho basis w.r.t. ambient coords
274 # The jordan product returns a matrixy answer, so we
275 # have to convert it to the algebra coordinates.
276 elt
= jordan_product(q_i
, q_j
)
277 elt
= W
.coordinate_vector(V(_all2list(elt
)))
278 self
._multiplication
_table
[i
][j
] = self
.from_vector(elt
)
280 if not orthonormalize
:
281 # If we're orthonormalizing the basis with respect
282 # to an inner-product, then the inner-product
283 # matrix with respect to the resulting basis is
284 # just going to be the identity.
285 ip
= inner_product(q_i
, q_j
)
286 self
._inner
_product
_matrix
[i
,j
] = ip
287 self
._inner
_product
_matrix
[j
,i
] = ip
289 self
._inner
_product
_matrix
._cache
= {'hermitian': True}
290 self
._inner
_product
_matrix
.set_immutable()
293 if not self
._is
_jordanian
():
294 raise ValueError("Jordan identity does not hold")
295 if not self
._inner
_product
_is
_associative
():
296 raise ValueError("inner product is not associative")
299 def _coerce_map_from_base_ring(self
):
301 Disable the map from the base ring into the algebra.
303 Performing a nonsense conversion like this automatically
304 is counterpedagogical. The fallback is to try the usual
305 element constructor, which should also fail.
309 sage: from mjo.eja.eja_algebra import random_eja
313 sage: set_random_seed()
314 sage: J = random_eja()
316 Traceback (most recent call last):
318 ValueError: not an element of this algebra
324 def product_on_basis(self
, i
, j
):
326 Returns the Jordan product of the `i` and `j`th basis elements.
328 This completely defines the Jordan product on the algebra, and
329 is used direclty by our superclass machinery to implement
334 sage: from mjo.eja.eja_algebra import random_eja
338 sage: set_random_seed()
339 sage: J = random_eja()
340 sage: n = J.dimension()
343 sage: bi_bj = J.zero()*J.zero()
345 ....: i = ZZ.random_element(n)
346 ....: j = ZZ.random_element(n)
347 ....: bi = J.monomial(i)
348 ....: bj = J.monomial(j)
349 ....: bi_bj = J.product_on_basis(i,j)
354 # We only stored the lower-triangular portion of the
355 # multiplication table.
357 return self
._multiplication
_table
[i
][j
]
359 return self
._multiplication
_table
[j
][i
]
361 def inner_product(self
, x
, y
):
363 The inner product associated with this Euclidean Jordan algebra.
365 Defaults to the trace inner product, but can be overridden by
366 subclasses if they are sure that the necessary properties are
371 sage: from mjo.eja.eja_algebra import (random_eja,
373 ....: BilinearFormEJA)
377 Our inner product is "associative," which means the following for
378 a symmetric bilinear form::
380 sage: set_random_seed()
381 sage: J = random_eja()
382 sage: x,y,z = J.random_elements(3)
383 sage: (x*y).inner_product(z) == y.inner_product(x*z)
388 Ensure that this is the usual inner product for the algebras
391 sage: set_random_seed()
392 sage: J = HadamardEJA.random_instance()
393 sage: x,y = J.random_elements(2)
394 sage: actual = x.inner_product(y)
395 sage: expected = x.to_vector().inner_product(y.to_vector())
396 sage: actual == expected
399 Ensure that this is one-half of the trace inner-product in a
400 BilinearFormEJA that isn't just the reals (when ``n`` isn't
401 one). This is in Faraut and Koranyi, and also my "On the
404 sage: set_random_seed()
405 sage: J = BilinearFormEJA.random_instance()
406 sage: n = J.dimension()
407 sage: x = J.random_element()
408 sage: y = J.random_element()
409 sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
413 B
= self
._inner
_product
_matrix
414 return (B
*x
.to_vector()).inner_product(y
.to_vector())
417 def is_associative(self
):
419 Return whether or not this algebra's Jordan product is associative.
423 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
427 sage: J = ComplexHermitianEJA(3, field=QQ, orthonormalize=False)
428 sage: J.is_associative()
430 sage: x = sum(J.gens())
431 sage: A = x.subalgebra_generated_by(orthonormalize=False)
432 sage: A.is_associative()
436 return "Associative" in self
.category().axioms()
438 def _is_commutative(self
):
440 Whether or not this algebra's multiplication table is commutative.
442 This method should of course always return ``True``, unless
443 this algebra was constructed with ``check_axioms=False`` and
444 passed an invalid multiplication table.
446 return all( x
*y
== y
*x
for x
in self
.gens() for y
in self
.gens() )
448 def _is_jordanian(self
):
450 Whether or not this algebra's multiplication table respects the
451 Jordan identity `(x^{2})(xy) = x(x^{2}y)`.
453 We only check one arrangement of `x` and `y`, so for a
454 ``True`` result to be truly true, you should also check
455 :meth:`_is_commutative`. This method should of course always
456 return ``True``, unless this algebra was constructed with
457 ``check_axioms=False`` and passed an invalid multiplication table.
459 return all( (self
.monomial(i
)**2)*(self
.monomial(i
)*self
.monomial(j
))
461 (self
.monomial(i
))*((self
.monomial(i
)**2)*self
.monomial(j
))
462 for i
in range(self
.dimension())
463 for j
in range(self
.dimension()) )
465 def _jordan_product_is_associative(self
):
467 Return whether or not this algebra's Jordan product is
468 associative; that is, whether or not `x*(y*z) = (x*y)*z`
471 This method should agree with :meth:`is_associative` unless
472 you lied about the value of the ``associative`` parameter
473 when you constructed the algebra.
477 sage: from mjo.eja.eja_algebra import (random_eja,
478 ....: RealSymmetricEJA,
479 ....: ComplexHermitianEJA,
480 ....: QuaternionHermitianEJA)
484 sage: J = RealSymmetricEJA(4, orthonormalize=False)
485 sage: J._jordan_product_is_associative()
487 sage: x = sum(J.gens())
488 sage: A = x.subalgebra_generated_by()
489 sage: A._jordan_product_is_associative()
494 sage: J = ComplexHermitianEJA(2,field=QQ,orthonormalize=False)
495 sage: J._jordan_product_is_associative()
497 sage: x = sum(J.gens())
498 sage: A = x.subalgebra_generated_by(orthonormalize=False)
499 sage: A._jordan_product_is_associative()
504 sage: J = QuaternionHermitianEJA(2)
505 sage: J._jordan_product_is_associative()
507 sage: x = sum(J.gens())
508 sage: A = x.subalgebra_generated_by()
509 sage: A._jordan_product_is_associative()
514 The values we've presupplied to the constructors agree with
517 sage: set_random_seed()
518 sage: J = random_eja()
519 sage: J.is_associative() == J._jordan_product_is_associative()
525 # Used to check whether or not something is zero.
528 # I don't know of any examples that make this magnitude
529 # necessary because I don't know how to make an
530 # associative algebra when the element subalgebra
531 # construction is unreliable (as it is over RDF; we can't
532 # find the degree of an element because we can't compute
533 # the rank of a matrix). But even multiplication of floats
534 # is non-associative, so *some* epsilon is needed... let's
535 # just take the one from _inner_product_is_associative?
538 for i
in range(self
.dimension()):
539 for j
in range(self
.dimension()):
540 for k
in range(self
.dimension()):
544 diff
= (x
*y
)*z
- x
*(y
*z
)
546 if diff
.norm() > epsilon
:
551 def _inner_product_is_associative(self
):
553 Return whether or not this algebra's inner product `B` is
554 associative; that is, whether or not `B(xy,z) = B(x,yz)`.
556 This method should of course always return ``True``, unless
557 this algebra was constructed with ``check_axioms=False`` and
558 passed an invalid Jordan or inner-product.
562 # Used to check whether or not something is zero.
565 # This choice is sufficient to allow the construction of
566 # QuaternionHermitianEJA(2, field=RDF) with check_axioms=True.
569 for i
in range(self
.dimension()):
570 for j
in range(self
.dimension()):
571 for k
in range(self
.dimension()):
575 diff
= (x
*y
).inner_product(z
) - x
.inner_product(y
*z
)
577 if diff
.abs() > epsilon
:
582 def _element_constructor_(self
, elt
):
584 Construct an element of this algebra from its vector or matrix
587 This gets called only after the parent element _call_ method
588 fails to find a coercion for the argument.
592 sage: from mjo.eja.eja_algebra import (random_eja,
595 ....: RealSymmetricEJA)
599 The identity in `S^n` is converted to the identity in the EJA::
601 sage: J = RealSymmetricEJA(3)
602 sage: I = matrix.identity(QQ,3)
603 sage: J(I) == J.one()
606 This skew-symmetric matrix can't be represented in the EJA::
608 sage: J = RealSymmetricEJA(3)
609 sage: A = matrix(QQ,3, lambda i,j: i-j)
611 Traceback (most recent call last):
613 ValueError: not an element of this algebra
615 Tuples work as well, provided that the matrix basis for the
616 algebra consists of them::
618 sage: J1 = HadamardEJA(3)
619 sage: J2 = RealSymmetricEJA(2)
620 sage: J = cartesian_product([J1,J2])
621 sage: J( (J1.matrix_basis()[1], J2.matrix_basis()[2]) )
626 Ensure that we can convert any element back and forth
627 faithfully between its matrix and algebra representations::
629 sage: set_random_seed()
630 sage: J = random_eja()
631 sage: x = J.random_element()
632 sage: J(x.to_matrix()) == x
635 We cannot coerce elements between algebras just because their
636 matrix representations are compatible::
638 sage: J1 = HadamardEJA(3)
639 sage: J2 = JordanSpinEJA(3)
641 Traceback (most recent call last):
643 ValueError: not an element of this algebra
645 Traceback (most recent call last):
647 ValueError: not an element of this algebra
649 msg
= "not an element of this algebra"
650 if elt
in self
.base_ring():
651 # Ensure that no base ring -> algebra coercion is performed
652 # by this method. There's some stupidity in sage that would
653 # otherwise propagate to this method; for example, sage thinks
654 # that the integer 3 belongs to the space of 2-by-2 matrices.
655 raise ValueError(msg
)
658 # Try to convert a vector into a column-matrix...
660 except (AttributeError, TypeError):
661 # and ignore failure, because we weren't really expecting
662 # a vector as an argument anyway.
665 if elt
not in self
.matrix_space():
666 raise ValueError(msg
)
668 # Thanks for nothing! Matrix spaces aren't vector spaces in
669 # Sage, so we have to figure out its matrix-basis coordinates
670 # ourselves. We use the basis space's ring instead of the
671 # element's ring because the basis space might be an algebraic
672 # closure whereas the base ring of the 3-by-3 identity matrix
673 # could be QQ instead of QQbar.
675 # And, we also have to handle Cartesian product bases (when
676 # the matrix basis consists of tuples) here. The "good news"
677 # is that we're already converting everything to long vectors,
678 # and that strategy works for tuples as well.
680 # We pass check=False because the matrix basis is "guaranteed"
681 # to be linearly independent... right? Ha ha.
683 V
= VectorSpace(self
.base_ring(), len(elt
))
684 W
= V
.span_of_basis( (V(_all2list(s
)) for s
in self
.matrix_basis()),
688 coords
= W
.coordinate_vector(V(elt
))
689 except ArithmeticError: # vector is not in free module
690 raise ValueError(msg
)
692 return self
.from_vector(coords
)
696 Return a string representation of ``self``.
700 sage: from mjo.eja.eja_algebra import JordanSpinEJA
704 Ensure that it says what we think it says::
706 sage: JordanSpinEJA(2, field=AA)
707 Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
708 sage: JordanSpinEJA(3, field=RDF)
709 Euclidean Jordan algebra of dimension 3 over Real Double Field
712 fmt
= "Euclidean Jordan algebra of dimension {} over {}"
713 return fmt
.format(self
.dimension(), self
.base_ring())
717 def characteristic_polynomial_of(self
):
719 Return the algebra's "characteristic polynomial of" function,
720 which is itself a multivariate polynomial that, when evaluated
721 at the coordinates of some algebra element, returns that
722 element's characteristic polynomial.
724 The resulting polynomial has `n+1` variables, where `n` is the
725 dimension of this algebra. The first `n` variables correspond to
726 the coordinates of an algebra element: when evaluated at the
727 coordinates of an algebra element with respect to a certain
728 basis, the result is a univariate polynomial (in the one
729 remaining variable ``t``), namely the characteristic polynomial
734 sage: from mjo.eja.eja_algebra import JordanSpinEJA, TrivialEJA
738 The characteristic polynomial in the spin algebra is given in
739 Alizadeh, Example 11.11::
741 sage: J = JordanSpinEJA(3)
742 sage: p = J.characteristic_polynomial_of(); p
743 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
744 sage: xvec = J.one().to_vector()
748 By definition, the characteristic polynomial is a monic
749 degree-zero polynomial in a rank-zero algebra. Note that
750 Cayley-Hamilton is indeed satisfied since the polynomial
751 ``1`` evaluates to the identity element of the algebra on
754 sage: J = TrivialEJA()
755 sage: J.characteristic_polynomial_of()
762 # The list of coefficient polynomials a_0, a_1, a_2, ..., a_(r-1).
763 a
= self
._charpoly
_coefficients
()
765 # We go to a bit of trouble here to reorder the
766 # indeterminates, so that it's easier to evaluate the
767 # characteristic polynomial at x's coordinates and get back
768 # something in terms of t, which is what we want.
769 S
= PolynomialRing(self
.base_ring(),'t')
773 S
= PolynomialRing(S
, R
.variable_names())
776 return (t
**r
+ sum( a
[k
]*(t
**k
) for k
in range(r
) ))
778 def coordinate_polynomial_ring(self
):
780 The multivariate polynomial ring in which this algebra's
781 :meth:`characteristic_polynomial_of` lives.
785 sage: from mjo.eja.eja_algebra import (HadamardEJA,
786 ....: RealSymmetricEJA)
790 sage: J = HadamardEJA(2)
791 sage: J.coordinate_polynomial_ring()
792 Multivariate Polynomial Ring in X1, X2...
793 sage: J = RealSymmetricEJA(3,field=QQ,orthonormalize=False)
794 sage: J.coordinate_polynomial_ring()
795 Multivariate Polynomial Ring in X1, X2, X3, X4, X5, X6...
798 var_names
= tuple( "X%d" % z
for z
in range(1, self
.dimension()+1) )
799 return PolynomialRing(self
.base_ring(), var_names
)
801 def inner_product(self
, x
, y
):
803 The inner product associated with this Euclidean Jordan algebra.
805 Defaults to the trace inner product, but can be overridden by
806 subclasses if they are sure that the necessary properties are
811 sage: from mjo.eja.eja_algebra import (random_eja,
813 ....: BilinearFormEJA)
817 Our inner product is "associative," which means the following for
818 a symmetric bilinear form::
820 sage: set_random_seed()
821 sage: J = random_eja()
822 sage: x,y,z = J.random_elements(3)
823 sage: (x*y).inner_product(z) == y.inner_product(x*z)
828 Ensure that this is the usual inner product for the algebras
831 sage: set_random_seed()
832 sage: J = HadamardEJA.random_instance()
833 sage: x,y = J.random_elements(2)
834 sage: actual = x.inner_product(y)
835 sage: expected = x.to_vector().inner_product(y.to_vector())
836 sage: actual == expected
839 Ensure that this is one-half of the trace inner-product in a
840 BilinearFormEJA that isn't just the reals (when ``n`` isn't
841 one). This is in Faraut and Koranyi, and also my "On the
844 sage: set_random_seed()
845 sage: J = BilinearFormEJA.random_instance()
846 sage: n = J.dimension()
847 sage: x = J.random_element()
848 sage: y = J.random_element()
849 sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
852 B
= self
._inner
_product
_matrix
853 return (B
*x
.to_vector()).inner_product(y
.to_vector())
856 def is_trivial(self
):
858 Return whether or not this algebra is trivial.
860 A trivial algebra contains only the zero element.
864 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
869 sage: J = ComplexHermitianEJA(3)
875 sage: J = TrivialEJA()
880 return self
.dimension() == 0
883 def multiplication_table(self
):
885 Return a visual representation of this algebra's multiplication
886 table (on basis elements).
890 sage: from mjo.eja.eja_algebra import JordanSpinEJA
894 sage: J = JordanSpinEJA(4)
895 sage: J.multiplication_table()
896 +----++----+----+----+----+
897 | * || b0 | b1 | b2 | b3 |
898 +====++====+====+====+====+
899 | b0 || b0 | b1 | b2 | b3 |
900 +----++----+----+----+----+
901 | b1 || b1 | b0 | 0 | 0 |
902 +----++----+----+----+----+
903 | b2 || b2 | 0 | b0 | 0 |
904 +----++----+----+----+----+
905 | b3 || b3 | 0 | 0 | b0 |
906 +----++----+----+----+----+
910 # Prepend the header row.
911 M
= [["*"] + list(self
.gens())]
913 # And to each subsequent row, prepend an entry that belongs to
914 # the left-side "header column."
915 M
+= [ [self
.monomial(i
)] + [ self
.monomial(i
)*self
.monomial(j
)
919 return table(M
, header_row
=True, header_column
=True, frame
=True)
922 def matrix_basis(self
):
924 Return an (often more natural) representation of this algebras
925 basis as an ordered tuple of matrices.
927 Every finite-dimensional Euclidean Jordan Algebra is a, up to
928 Jordan isomorphism, a direct sum of five simple
929 algebras---four of which comprise Hermitian matrices. And the
930 last type of algebra can of course be thought of as `n`-by-`1`
931 column matrices (ambiguusly called column vectors) to avoid
932 special cases. As a result, matrices (and column vectors) are
933 a natural representation format for Euclidean Jordan algebra
936 But, when we construct an algebra from a basis of matrices,
937 those matrix representations are lost in favor of coordinate
938 vectors *with respect to* that basis. We could eventually
939 convert back if we tried hard enough, but having the original
940 representations handy is valuable enough that we simply store
941 them and return them from this method.
943 Why implement this for non-matrix algebras? Avoiding special
944 cases for the :class:`BilinearFormEJA` pays with simplicity in
945 its own right. But mainly, we would like to be able to assume
946 that elements of a :class:`CartesianProductEJA` can be displayed
947 nicely, without having to have special classes for direct sums
948 one of whose components was a matrix algebra.
952 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
953 ....: RealSymmetricEJA)
957 sage: J = RealSymmetricEJA(2)
959 Finite family {0: b0, 1: b1, 2: b2}
960 sage: J.matrix_basis()
962 [1 0] [ 0 0.7071067811865475?] [0 0]
963 [0 0], [0.7071067811865475? 0], [0 1]
968 sage: J = JordanSpinEJA(2)
970 Finite family {0: b0, 1: b1}
971 sage: J.matrix_basis()
977 return self
._matrix
_basis
980 def matrix_space(self
):
982 Return the matrix space in which this algebra's elements live, if
983 we think of them as matrices (including column vectors of the
986 "By default" this will be an `n`-by-`1` column-matrix space,
987 except when the algebra is trivial. There it's `n`-by-`n`
988 (where `n` is zero), to ensure that two elements of the matrix
989 space (empty matrices) can be multiplied. For algebras of
990 matrices, this returns the space in which their
991 real embeddings live.
995 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
997 ....: QuaternionHermitianEJA,
1002 By default, the matrix representation is just a column-matrix
1003 equivalent to the vector representation::
1005 sage: J = JordanSpinEJA(3)
1006 sage: J.matrix_space()
1007 Full MatrixSpace of 3 by 1 dense matrices over Algebraic
1010 The matrix representation in the trivial algebra is
1011 zero-by-zero instead of the usual `n`-by-one::
1013 sage: J = TrivialEJA()
1014 sage: J.matrix_space()
1015 Full MatrixSpace of 0 by 0 dense matrices over Algebraic
1018 The matrix space for complex/quaternion Hermitian matrix EJA
1019 is the space in which their real-embeddings live, not the
1020 original complex/quaternion matrix space::
1022 sage: J = ComplexHermitianEJA(2,field=QQ,orthonormalize=False)
1023 sage: J.matrix_space()
1024 Module of 2 by 2 matrices with entries in Algebraic Field over
1025 the scalar ring Rational Field
1026 sage: J = QuaternionHermitianEJA(1,field=QQ,orthonormalize=False)
1027 sage: J.matrix_space()
1028 Module of 1 by 1 matrices with entries in Quaternion
1029 Algebra (-1, -1) with base ring Rational Field over
1030 the scalar ring Rational Field
1033 if self
.is_trivial():
1034 return MatrixSpace(self
.base_ring(), 0)
1036 return self
.matrix_basis()[0].parent()
1042 Return the unit element of this algebra.
1046 sage: from mjo.eja.eja_algebra import (HadamardEJA,
1051 We can compute unit element in the Hadamard EJA::
1053 sage: J = HadamardEJA(5)
1055 b0 + b1 + b2 + b3 + b4
1057 The unit element in the Hadamard EJA is inherited in the
1058 subalgebras generated by its elements::
1060 sage: J = HadamardEJA(5)
1062 b0 + b1 + b2 + b3 + b4
1063 sage: x = sum(J.gens())
1064 sage: A = x.subalgebra_generated_by(orthonormalize=False)
1067 sage: A.one().superalgebra_element()
1068 b0 + b1 + b2 + b3 + b4
1072 The identity element acts like the identity, regardless of
1073 whether or not we orthonormalize::
1075 sage: set_random_seed()
1076 sage: J = random_eja()
1077 sage: x = J.random_element()
1078 sage: J.one()*x == x and x*J.one() == x
1080 sage: A = x.subalgebra_generated_by()
1081 sage: y = A.random_element()
1082 sage: A.one()*y == y and y*A.one() == y
1087 sage: set_random_seed()
1088 sage: J = random_eja(field=QQ, orthonormalize=False)
1089 sage: x = J.random_element()
1090 sage: J.one()*x == x and x*J.one() == x
1092 sage: A = x.subalgebra_generated_by(orthonormalize=False)
1093 sage: y = A.random_element()
1094 sage: A.one()*y == y and y*A.one() == y
1097 The matrix of the unit element's operator is the identity,
1098 regardless of the base field and whether or not we
1101 sage: set_random_seed()
1102 sage: J = random_eja()
1103 sage: actual = J.one().operator().matrix()
1104 sage: expected = matrix.identity(J.base_ring(), J.dimension())
1105 sage: actual == expected
1107 sage: x = J.random_element()
1108 sage: A = x.subalgebra_generated_by()
1109 sage: actual = A.one().operator().matrix()
1110 sage: expected = matrix.identity(A.base_ring(), A.dimension())
1111 sage: actual == expected
1116 sage: set_random_seed()
1117 sage: J = random_eja(field=QQ, orthonormalize=False)
1118 sage: actual = J.one().operator().matrix()
1119 sage: expected = matrix.identity(J.base_ring(), J.dimension())
1120 sage: actual == expected
1122 sage: x = J.random_element()
1123 sage: A = x.subalgebra_generated_by(orthonormalize=False)
1124 sage: actual = A.one().operator().matrix()
1125 sage: expected = matrix.identity(A.base_ring(), A.dimension())
1126 sage: actual == expected
1129 Ensure that the cached unit element (often precomputed by
1130 hand) agrees with the computed one::
1132 sage: set_random_seed()
1133 sage: J = random_eja()
1134 sage: cached = J.one()
1135 sage: J.one.clear_cache()
1136 sage: J.one() == cached
1141 sage: set_random_seed()
1142 sage: J = random_eja(field=QQ, orthonormalize=False)
1143 sage: cached = J.one()
1144 sage: J.one.clear_cache()
1145 sage: J.one() == cached
1149 # We can brute-force compute the matrices of the operators
1150 # that correspond to the basis elements of this algebra.
1151 # If some linear combination of those basis elements is the
1152 # algebra identity, then the same linear combination of
1153 # their matrices has to be the identity matrix.
1155 # Of course, matrices aren't vectors in sage, so we have to
1156 # appeal to the "long vectors" isometry.
1157 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
1159 # Now we use basic linear algebra to find the coefficients,
1160 # of the matrices-as-vectors-linear-combination, which should
1161 # work for the original algebra basis too.
1162 A
= matrix(self
.base_ring(), oper_vecs
)
1164 # We used the isometry on the left-hand side already, but we
1165 # still need to do it for the right-hand side. Recall that we
1166 # wanted something that summed to the identity matrix.
1167 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
1169 # Now if there's an identity element in the algebra, this
1170 # should work. We solve on the left to avoid having to
1171 # transpose the matrix "A".
1172 return self
.from_vector(A
.solve_left(b
))
1175 def peirce_decomposition(self
, c
):
1177 The Peirce decomposition of this algebra relative to the
1180 In the future, this can be extended to a complete system of
1181 orthogonal idempotents.
1185 - ``c`` -- an idempotent of this algebra.
1189 A triple (J0, J5, J1) containing two subalgebras and one subspace
1192 - ``J0`` -- the algebra on the eigenspace of ``c.operator()``
1193 corresponding to the eigenvalue zero.
1195 - ``J5`` -- the eigenspace (NOT a subalgebra) of ``c.operator()``
1196 corresponding to the eigenvalue one-half.
1198 - ``J1`` -- the algebra on the eigenspace of ``c.operator()``
1199 corresponding to the eigenvalue one.
1201 These are the only possible eigenspaces for that operator, and this
1202 algebra is a direct sum of them. The spaces ``J0`` and ``J1`` are
1203 orthogonal, and are subalgebras of this algebra with the appropriate
1208 sage: from mjo.eja.eja_algebra import random_eja, RealSymmetricEJA
1212 The canonical example comes from the symmetric matrices, which
1213 decompose into diagonal and off-diagonal parts::
1215 sage: J = RealSymmetricEJA(3)
1216 sage: C = matrix(QQ, [ [1,0,0],
1220 sage: J0,J5,J1 = J.peirce_decomposition(c)
1222 Euclidean Jordan algebra of dimension 1...
1224 Vector space of degree 6 and dimension 2...
1226 Euclidean Jordan algebra of dimension 3...
1227 sage: J0.one().to_matrix()
1231 sage: orig_df = AA.options.display_format
1232 sage: AA.options.display_format = 'radical'
1233 sage: J.from_vector(J5.basis()[0]).to_matrix()
1237 sage: J.from_vector(J5.basis()[1]).to_matrix()
1241 sage: AA.options.display_format = orig_df
1242 sage: J1.one().to_matrix()
1249 Every algebra decomposes trivially with respect to its identity
1252 sage: set_random_seed()
1253 sage: J = random_eja()
1254 sage: J0,J5,J1 = J.peirce_decomposition(J.one())
1255 sage: J0.dimension() == 0 and J5.dimension() == 0
1257 sage: J1.superalgebra() == J and J1.dimension() == J.dimension()
1260 The decomposition is into eigenspaces, and its components are
1261 therefore necessarily orthogonal. Moreover, the identity
1262 elements in the two subalgebras are the projections onto their
1263 respective subspaces of the superalgebra's identity element::
1265 sage: set_random_seed()
1266 sage: J = random_eja()
1267 sage: x = J.random_element()
1268 sage: if not J.is_trivial():
1269 ....: while x.is_nilpotent():
1270 ....: x = J.random_element()
1271 sage: c = x.subalgebra_idempotent()
1272 sage: J0,J5,J1 = J.peirce_decomposition(c)
1274 sage: for (w,y,z) in zip(J0.basis(), J5.basis(), J1.basis()):
1275 ....: w = w.superalgebra_element()
1276 ....: y = J.from_vector(y)
1277 ....: z = z.superalgebra_element()
1278 ....: ipsum += w.inner_product(y).abs()
1279 ....: ipsum += w.inner_product(z).abs()
1280 ....: ipsum += y.inner_product(z).abs()
1283 sage: J1(c) == J1.one()
1285 sage: J0(J.one() - c) == J0.one()
1289 if not c
.is_idempotent():
1290 raise ValueError("element is not idempotent: %s" % c
)
1292 # Default these to what they should be if they turn out to be
1293 # trivial, because eigenspaces_left() won't return eigenvalues
1294 # corresponding to trivial spaces (e.g. it returns only the
1295 # eigenspace corresponding to lambda=1 if you take the
1296 # decomposition relative to the identity element).
1297 trivial
= self
.subalgebra(())
1298 J0
= trivial
# eigenvalue zero
1299 J5
= VectorSpace(self
.base_ring(), 0) # eigenvalue one-half
1300 J1
= trivial
# eigenvalue one
1302 for (eigval
, eigspace
) in c
.operator().matrix().right_eigenspaces():
1303 if eigval
== ~
(self
.base_ring()(2)):
1306 gens
= tuple( self
.from_vector(b
) for b
in eigspace
.basis() )
1307 subalg
= self
.subalgebra(gens
, check_axioms
=False)
1313 raise ValueError("unexpected eigenvalue: %s" % eigval
)
1318 def random_element(self
, thorough
=False):
1320 Return a random element of this algebra.
1322 Our algebra superclass method only returns a linear
1323 combination of at most two basis elements. We instead
1324 want the vector space "random element" method that
1325 returns a more diverse selection.
1329 - ``thorough`` -- (boolean; default False) whether or not we
1330 should generate irrational coefficients for the random
1331 element when our base ring is irrational; this slows the
1332 algebra operations to a crawl, but any truly random method
1336 # For a general base ring... maybe we can trust this to do the
1337 # right thing? Unlikely, but.
1338 V
= self
.vector_space()
1339 v
= V
.random_element()
1341 if self
.base_ring() is AA
:
1342 # The "random element" method of the algebraic reals is
1343 # stupid at the moment, and only returns integers between
1344 # -2 and 2, inclusive:
1346 # https://trac.sagemath.org/ticket/30875
1348 # Instead, we implement our own "random vector" method,
1349 # and then coerce that into the algebra. We use the vector
1350 # space degree here instead of the dimension because a
1351 # subalgebra could (for example) be spanned by only two
1352 # vectors, each with five coordinates. We need to
1353 # generate all five coordinates.
1355 v
*= QQbar
.random_element().real()
1357 v
*= QQ
.random_element()
1359 return self
.from_vector(V
.coordinate_vector(v
))
1361 def random_elements(self
, count
, thorough
=False):
1363 Return ``count`` random elements as a tuple.
1367 - ``thorough`` -- (boolean; default False) whether or not we
1368 should generate irrational coefficients for the random
1369 elements when our base ring is irrational; this slows the
1370 algebra operations to a crawl, but any truly random method
1375 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1379 sage: J = JordanSpinEJA(3)
1380 sage: x,y,z = J.random_elements(3)
1381 sage: all( [ x in J, y in J, z in J ])
1383 sage: len( J.random_elements(10) ) == 10
1387 return tuple( self
.random_element(thorough
)
1388 for idx
in range(count
) )
1392 def _charpoly_coefficients(self
):
1394 The `r` polynomial coefficients of the "characteristic polynomial
1399 sage: from mjo.eja.eja_algebra import random_eja
1403 The theory shows that these are all homogeneous polynomials of
1406 sage: set_random_seed()
1407 sage: J = random_eja()
1408 sage: all(p.is_homogeneous() for p in J._charpoly_coefficients())
1412 n
= self
.dimension()
1413 R
= self
.coordinate_polynomial_ring()
1415 F
= R
.fraction_field()
1418 # From a result in my book, these are the entries of the
1419 # basis representation of L_x.
1420 return sum( vars[k
]*self
.monomial(k
).operator().matrix()[i
,j
]
1423 L_x
= matrix(F
, n
, n
, L_x_i_j
)
1426 if self
.rank
.is_in_cache():
1428 # There's no need to pad the system with redundant
1429 # columns if we *know* they'll be redundant.
1432 # Compute an extra power in case the rank is equal to
1433 # the dimension (otherwise, we would stop at x^(r-1)).
1434 x_powers
= [ (L_x
**k
)*self
.one().to_vector()
1435 for k
in range(n
+1) ]
1436 A
= matrix
.column(F
, x_powers
[:n
])
1437 AE
= A
.extended_echelon_form()
1444 # The theory says that only the first "r" coefficients are
1445 # nonzero, and they actually live in the original polynomial
1446 # ring and not the fraction field. We negate them because in
1447 # the actual characteristic polynomial, they get moved to the
1448 # other side where x^r lives. We don't bother to trim A_rref
1449 # down to a square matrix and solve the resulting system,
1450 # because the upper-left r-by-r portion of A_rref is
1451 # guaranteed to be the identity matrix, so e.g.
1453 # A_rref.solve_right(Y)
1455 # would just be returning Y.
1456 return (-E
*b
)[:r
].change_ring(R
)
1461 Return the rank of this EJA.
1463 This is a cached method because we know the rank a priori for
1464 all of the algebras we can construct. Thus we can avoid the
1465 expensive ``_charpoly_coefficients()`` call unless we truly
1466 need to compute the whole characteristic polynomial.
1470 sage: from mjo.eja.eja_algebra import (HadamardEJA,
1471 ....: JordanSpinEJA,
1472 ....: RealSymmetricEJA,
1473 ....: ComplexHermitianEJA,
1474 ....: QuaternionHermitianEJA,
1479 The rank of the Jordan spin algebra is always two::
1481 sage: JordanSpinEJA(2).rank()
1483 sage: JordanSpinEJA(3).rank()
1485 sage: JordanSpinEJA(4).rank()
1488 The rank of the `n`-by-`n` Hermitian real, complex, or
1489 quaternion matrices is `n`::
1491 sage: RealSymmetricEJA(4).rank()
1493 sage: ComplexHermitianEJA(3).rank()
1495 sage: QuaternionHermitianEJA(2).rank()
1500 Ensure that every EJA that we know how to construct has a
1501 positive integer rank, unless the algebra is trivial in
1502 which case its rank will be zero::
1504 sage: set_random_seed()
1505 sage: J = random_eja()
1509 sage: r > 0 or (r == 0 and J.is_trivial())
1512 Ensure that computing the rank actually works, since the ranks
1513 of all simple algebras are known and will be cached by default::
1515 sage: set_random_seed() # long time
1516 sage: J = random_eja() # long time
1517 sage: cached = J.rank() # long time
1518 sage: J.rank.clear_cache() # long time
1519 sage: J.rank() == cached # long time
1523 return len(self
._charpoly
_coefficients
())
1526 def subalgebra(self
, basis
, **kwargs
):
1528 Create a subalgebra of this algebra from the given basis.
1530 from mjo
.eja
.eja_subalgebra
import FiniteDimensionalEJASubalgebra
1531 return FiniteDimensionalEJASubalgebra(self
, basis
, **kwargs
)
1534 def vector_space(self
):
1536 Return the vector space that underlies this algebra.
1540 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1544 sage: J = RealSymmetricEJA(2)
1545 sage: J.vector_space()
1546 Vector space of dimension 3 over...
1549 return self
.zero().to_vector().parent().ambient_vector_space()
1553 class RationalBasisEJA(FiniteDimensionalEJA
):
1555 Algebras whose supplied basis elements have all rational entries.
1559 sage: from mjo.eja.eja_algebra import BilinearFormEJA
1563 The supplied basis is orthonormalized by default::
1565 sage: B = matrix(QQ, [[1, 0, 0], [0, 25, -32], [0, -32, 41]])
1566 sage: J = BilinearFormEJA(B)
1567 sage: J.matrix_basis()
1584 # Abuse the check_field parameter to check that the entries of
1585 # out basis (in ambient coordinates) are in the field QQ.
1586 # Use _all2list to get the vector coordinates of octonion
1587 # entries and not the octonions themselves (which are not
1589 if not all( all(b_i
in QQ
for b_i
in _all2list(b
))
1591 raise TypeError("basis not rational")
1593 super().__init
__(basis
,
1597 check_field
=check_field
,
1600 self
._rational
_algebra
= None
1602 # There's no point in constructing the extra algebra if this
1603 # one is already rational.
1605 # Note: the same Jordan and inner-products work here,
1606 # because they are necessarily defined with respect to
1607 # ambient coordinates and not any particular basis.
1608 self
._rational
_algebra
= FiniteDimensionalEJA(
1613 associative
=self
.is_associative(),
1614 orthonormalize
=False,
1619 def _charpoly_coefficients(self
):
1623 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
1624 ....: JordanSpinEJA)
1628 The base ring of the resulting polynomial coefficients is what
1629 it should be, and not the rationals (unless the algebra was
1630 already over the rationals)::
1632 sage: J = JordanSpinEJA(3)
1633 sage: J._charpoly_coefficients()
1634 (X1^2 - X2^2 - X3^2, -2*X1)
1635 sage: a0 = J._charpoly_coefficients()[0]
1637 Algebraic Real Field
1638 sage: a0.base_ring()
1639 Algebraic Real Field
1642 if self
._rational
_algebra
is None:
1643 # There's no need to construct *another* algebra over the
1644 # rationals if this one is already over the
1645 # rationals. Likewise, if we never orthonormalized our
1646 # basis, we might as well just use the given one.
1647 return super()._charpoly
_coefficients
()
1649 # Do the computation over the rationals. The answer will be
1650 # the same, because all we've done is a change of basis.
1651 # Then, change back from QQ to our real base ring
1652 a
= ( a_i
.change_ring(self
.base_ring())
1653 for a_i
in self
._rational
_algebra
._charpoly
_coefficients
() )
1655 if self
._deortho
_matrix
is None:
1656 # This can happen if our base ring was, say, AA and we
1657 # chose not to (or didn't need to) orthonormalize. It's
1658 # still faster to do the computations over QQ even if
1659 # the numbers in the boxes stay the same.
1662 # Otherwise, convert the coordinate variables back to the
1663 # deorthonormalized ones.
1664 R
= self
.coordinate_polynomial_ring()
1665 from sage
.modules
.free_module_element
import vector
1666 X
= vector(R
, R
.gens())
1667 BX
= self
._deortho
_matrix
*X
1669 subs_dict
= { X[i]: BX[i] for i in range(len(X)) }
1670 return tuple( a_i
.subs(subs_dict
) for a_i
in a
)
1672 class ConcreteEJA(FiniteDimensionalEJA
):
1674 A class for the Euclidean Jordan algebras that we know by name.
1676 These are the Jordan algebras whose basis, multiplication table,
1677 rank, and so on are known a priori. More to the point, they are
1678 the Euclidean Jordan algebras for which we are able to conjure up
1679 a "random instance."
1683 sage: from mjo.eja.eja_algebra import ConcreteEJA
1687 Our basis is normalized with respect to the algebra's inner
1688 product, unless we specify otherwise::
1690 sage: set_random_seed()
1691 sage: J = ConcreteEJA.random_instance()
1692 sage: all( b.norm() == 1 for b in J.gens() )
1695 Since our basis is orthonormal with respect to the algebra's inner
1696 product, and since we know that this algebra is an EJA, any
1697 left-multiplication operator's matrix will be symmetric because
1698 natural->EJA basis representation is an isometry and within the
1699 EJA the operator is self-adjoint by the Jordan axiom::
1701 sage: set_random_seed()
1702 sage: J = ConcreteEJA.random_instance()
1703 sage: x = J.random_element()
1704 sage: x.operator().is_self_adjoint()
1709 def _max_random_instance_size():
1711 Return an integer "size" that is an upper bound on the size of
1712 this algebra when it is used in a random test
1713 case. Unfortunately, the term "size" is ambiguous -- when
1714 dealing with `R^n` under either the Hadamard or Jordan spin
1715 product, the "size" refers to the dimension `n`. When dealing
1716 with a matrix algebra (real symmetric or complex/quaternion
1717 Hermitian), it refers to the size of the matrix, which is far
1718 less than the dimension of the underlying vector space.
1720 This method must be implemented in each subclass.
1722 raise NotImplementedError
1725 def random_instance(cls
, *args
, **kwargs
):
1727 Return a random instance of this type of algebra.
1729 This method should be implemented in each subclass.
1731 from sage
.misc
.prandom
import choice
1732 eja_class
= choice(cls
.__subclasses
__())
1734 # These all bubble up to the RationalBasisEJA superclass
1735 # constructor, so any (kw)args valid there are also valid
1737 return eja_class
.random_instance(*args
, **kwargs
)
1742 def _denormalized_basis(A
):
1744 Returns a basis for the space of complex Hermitian n-by-n matrices.
1746 Why do we embed these? Basically, because all of numerical linear
1747 algebra assumes that you're working with vectors consisting of `n`
1748 entries from a field and scalars from the same field. There's no way
1749 to tell SageMath that (for example) the vectors contain complex
1750 numbers, while the scalar field is real.
1754 sage: from mjo.hurwitz import (ComplexMatrixAlgebra,
1755 ....: QuaternionMatrixAlgebra,
1756 ....: OctonionMatrixAlgebra)
1757 sage: from mjo.eja.eja_algebra import MatrixEJA
1761 sage: set_random_seed()
1762 sage: n = ZZ.random_element(1,5)
1763 sage: A = MatrixSpace(QQ, n)
1764 sage: B = MatrixEJA._denormalized_basis(A)
1765 sage: all( M.is_hermitian() for M in B)
1770 sage: set_random_seed()
1771 sage: n = ZZ.random_element(1,5)
1772 sage: A = ComplexMatrixAlgebra(n, scalars=QQ)
1773 sage: B = MatrixEJA._denormalized_basis(A)
1774 sage: all( M.is_hermitian() for M in B)
1779 sage: set_random_seed()
1780 sage: n = ZZ.random_element(1,5)
1781 sage: A = QuaternionMatrixAlgebra(n, scalars=QQ)
1782 sage: B = MatrixEJA._denormalized_basis(A)
1783 sage: all( M.is_hermitian() for M in B )
1788 sage: set_random_seed()
1789 sage: n = ZZ.random_element(1,5)
1790 sage: A = OctonionMatrixAlgebra(n, scalars=QQ)
1791 sage: B = MatrixEJA._denormalized_basis(A)
1792 sage: all( M.is_hermitian() for M in B )
1796 # These work for real MatrixSpace, whose monomials only have
1797 # two coordinates (because the last one would always be "1").
1798 es
= A
.base_ring().gens()
1799 gen
= lambda A
,m
: A
.monomial(m
[:2])
1801 if hasattr(A
, 'entry_algebra_gens'):
1802 # We've got a MatrixAlgebra, and its monomials will have
1803 # three coordinates.
1804 es
= A
.entry_algebra_gens()
1805 gen
= lambda A
,m
: A
.monomial(m
)
1808 for i
in range(A
.nrows()):
1809 for j
in range(i
+1):
1811 E_ii
= gen(A
, (i
,j
,es
[0]))
1815 E_ij
= gen(A
, (i
,j
,e
))
1816 E_ij
+= E_ij
.conjugate_transpose()
1819 return tuple( basis
)
1822 def jordan_product(X
,Y
):
1823 return (X
*Y
+ Y
*X
)/2
1826 def trace_inner_product(X
,Y
):
1828 A trace inner-product for matrices that aren't embedded in the
1829 reals. It takes MATRICES as arguments, not EJA elements.
1833 sage: from mjo.eja.eja_algebra import (RealSymmetricEJA,
1834 ....: ComplexHermitianEJA,
1835 ....: QuaternionHermitianEJA,
1836 ....: OctonionHermitianEJA)
1840 sage: J = RealSymmetricEJA(2,field=QQ,orthonormalize=False)
1841 sage: I = J.one().to_matrix()
1842 sage: J.trace_inner_product(I, -I)
1847 sage: J = ComplexHermitianEJA(2,field=QQ,orthonormalize=False)
1848 sage: I = J.one().to_matrix()
1849 sage: J.trace_inner_product(I, -I)
1854 sage: J = QuaternionHermitianEJA(2,field=QQ,orthonormalize=False)
1855 sage: I = J.one().to_matrix()
1856 sage: J.trace_inner_product(I, -I)
1861 sage: J = OctonionHermitianEJA(2,field=QQ,orthonormalize=False)
1862 sage: I = J.one().to_matrix()
1863 sage: J.trace_inner_product(I, -I)
1868 if hasattr(tr
, 'coefficient'):
1869 # Works for octonions, and has to come first because they
1870 # also have a "real()" method that doesn't return an
1871 # element of the scalar ring.
1872 return tr
.coefficient(0)
1873 elif hasattr(tr
, 'coefficient_tuple'):
1874 # Works for quaternions.
1875 return tr
.coefficient_tuple()[0]
1877 # Works for real and complex numbers.
1882 class RealSymmetricEJA(RationalBasisEJA
, ConcreteEJA
, MatrixEJA
):
1884 The rank-n simple EJA consisting of real symmetric n-by-n
1885 matrices, the usual symmetric Jordan product, and the trace inner
1886 product. It has dimension `(n^2 + n)/2` over the reals.
1890 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1894 sage: J = RealSymmetricEJA(2)
1895 sage: b0, b1, b2 = J.gens()
1903 In theory, our "field" can be any subfield of the reals::
1905 sage: RealSymmetricEJA(2, field=RDF, check_axioms=True)
1906 Euclidean Jordan algebra of dimension 3 over Real Double Field
1907 sage: RealSymmetricEJA(2, field=RR, check_axioms=True)
1908 Euclidean Jordan algebra of dimension 3 over Real Field with
1909 53 bits of precision
1913 The dimension of this algebra is `(n^2 + n) / 2`::
1915 sage: set_random_seed()
1916 sage: n_max = RealSymmetricEJA._max_random_instance_size()
1917 sage: n = ZZ.random_element(1, n_max)
1918 sage: J = RealSymmetricEJA(n)
1919 sage: J.dimension() == (n^2 + n)/2
1922 The Jordan multiplication is what we think it is::
1924 sage: set_random_seed()
1925 sage: J = RealSymmetricEJA.random_instance()
1926 sage: x,y = J.random_elements(2)
1927 sage: actual = (x*y).to_matrix()
1928 sage: X = x.to_matrix()
1929 sage: Y = y.to_matrix()
1930 sage: expected = (X*Y + Y*X)/2
1931 sage: actual == expected
1933 sage: J(expected) == x*y
1936 We can change the generator prefix::
1938 sage: RealSymmetricEJA(3, prefix='q').gens()
1939 (q0, q1, q2, q3, q4, q5)
1941 We can construct the (trivial) algebra of rank zero::
1943 sage: RealSymmetricEJA(0)
1944 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1948 def _max_random_instance_size():
1949 return 4 # Dimension 10
1952 def random_instance(cls
, **kwargs
):
1954 Return a random instance of this type of algebra.
1956 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1957 return cls(n
, **kwargs
)
1959 def __init__(self
, n
, field
=AA
, **kwargs
):
1960 # We know this is a valid EJA, but will double-check
1961 # if the user passes check_axioms=True.
1962 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
1964 A
= MatrixSpace(field
, n
)
1965 super().__init
__(self
._denormalized
_basis
(A
),
1966 self
.jordan_product
,
1967 self
.trace_inner_product
,
1971 # TODO: this could be factored out somehow, but is left here
1972 # because the MatrixEJA is not presently a subclass of the
1973 # FDEJA class that defines rank() and one().
1974 self
.rank
.set_cache(n
)
1976 self
.one
.set_cache( self
.zero() )
1978 self
.one
.set_cache(self(A
.one()))
1982 class ComplexHermitianEJA(RationalBasisEJA
, ConcreteEJA
, MatrixEJA
):
1984 The rank-n simple EJA consisting of complex Hermitian n-by-n
1985 matrices over the real numbers, the usual symmetric Jordan product,
1986 and the real-part-of-trace inner product. It has dimension `n^2` over
1991 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1995 In theory, our "field" can be any subfield of the reals::
1997 sage: ComplexHermitianEJA(2, field=RDF, check_axioms=True)
1998 Euclidean Jordan algebra of dimension 4 over Real Double Field
1999 sage: ComplexHermitianEJA(2, field=RR, check_axioms=True)
2000 Euclidean Jordan algebra of dimension 4 over Real Field with
2001 53 bits of precision
2005 The dimension of this algebra is `n^2`::
2007 sage: set_random_seed()
2008 sage: n_max = ComplexHermitianEJA._max_random_instance_size()
2009 sage: n = ZZ.random_element(1, n_max)
2010 sage: J = ComplexHermitianEJA(n)
2011 sage: J.dimension() == n^2
2014 The Jordan multiplication is what we think it is::
2016 sage: set_random_seed()
2017 sage: J = ComplexHermitianEJA.random_instance()
2018 sage: x,y = J.random_elements(2)
2019 sage: actual = (x*y).to_matrix()
2020 sage: X = x.to_matrix()
2021 sage: Y = y.to_matrix()
2022 sage: expected = (X*Y + Y*X)/2
2023 sage: actual == expected
2025 sage: J(expected) == x*y
2028 We can change the generator prefix::
2030 sage: ComplexHermitianEJA(2, prefix='z').gens()
2033 We can construct the (trivial) algebra of rank zero::
2035 sage: ComplexHermitianEJA(0)
2036 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2039 def __init__(self
, n
, field
=AA
, **kwargs
):
2040 # We know this is a valid EJA, but will double-check
2041 # if the user passes check_axioms=True.
2042 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2044 from mjo
.hurwitz
import ComplexMatrixAlgebra
2045 A
= ComplexMatrixAlgebra(n
, scalars
=field
)
2046 super().__init
__(self
._denormalized
_basis
(A
),
2047 self
.jordan_product
,
2048 self
.trace_inner_product
,
2052 # TODO: this could be factored out somehow, but is left here
2053 # because the MatrixEJA is not presently a subclass of the
2054 # FDEJA class that defines rank() and one().
2055 self
.rank
.set_cache(n
)
2057 self
.one
.set_cache( self
.zero() )
2059 self
.one
.set_cache(self(A
.one()))
2062 def _max_random_instance_size():
2063 return 3 # Dimension 9
2066 def random_instance(cls
, **kwargs
):
2068 Return a random instance of this type of algebra.
2070 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2071 return cls(n
, **kwargs
)
2074 class QuaternionHermitianEJA(RationalBasisEJA
, ConcreteEJA
, MatrixEJA
):
2076 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
2077 matrices, the usual symmetric Jordan product, and the
2078 real-part-of-trace inner product. It has dimension `2n^2 - n` over
2083 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2087 In theory, our "field" can be any subfield of the reals::
2089 sage: QuaternionHermitianEJA(2, field=RDF, check_axioms=True)
2090 Euclidean Jordan algebra of dimension 6 over Real Double Field
2091 sage: QuaternionHermitianEJA(2, field=RR, check_axioms=True)
2092 Euclidean Jordan algebra of dimension 6 over Real Field with
2093 53 bits of precision
2097 The dimension of this algebra is `2*n^2 - n`::
2099 sage: set_random_seed()
2100 sage: n_max = QuaternionHermitianEJA._max_random_instance_size()
2101 sage: n = ZZ.random_element(1, n_max)
2102 sage: J = QuaternionHermitianEJA(n)
2103 sage: J.dimension() == 2*(n^2) - n
2106 The Jordan multiplication is what we think it is::
2108 sage: set_random_seed()
2109 sage: J = QuaternionHermitianEJA.random_instance()
2110 sage: x,y = J.random_elements(2)
2111 sage: actual = (x*y).to_matrix()
2112 sage: X = x.to_matrix()
2113 sage: Y = y.to_matrix()
2114 sage: expected = (X*Y + Y*X)/2
2115 sage: actual == expected
2117 sage: J(expected) == x*y
2120 We can change the generator prefix::
2122 sage: QuaternionHermitianEJA(2, prefix='a').gens()
2123 (a0, a1, a2, a3, a4, a5)
2125 We can construct the (trivial) algebra of rank zero::
2127 sage: QuaternionHermitianEJA(0)
2128 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2131 def __init__(self
, n
, field
=AA
, **kwargs
):
2132 # We know this is a valid EJA, but will double-check
2133 # if the user passes check_axioms=True.
2134 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2136 from mjo
.hurwitz
import QuaternionMatrixAlgebra
2137 A
= QuaternionMatrixAlgebra(n
, scalars
=field
)
2138 super().__init
__(self
._denormalized
_basis
(A
),
2139 self
.jordan_product
,
2140 self
.trace_inner_product
,
2144 # TODO: this could be factored out somehow, but is left here
2145 # because the MatrixEJA is not presently a subclass of the
2146 # FDEJA class that defines rank() and one().
2147 self
.rank
.set_cache(n
)
2149 self
.one
.set_cache( self
.zero() )
2151 self
.one
.set_cache(self(A
.one()))
2155 def _max_random_instance_size():
2157 The maximum rank of a random QuaternionHermitianEJA.
2159 return 2 # Dimension 6
2162 def random_instance(cls
, **kwargs
):
2164 Return a random instance of this type of algebra.
2166 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2167 return cls(n
, **kwargs
)
2169 class OctonionHermitianEJA(RationalBasisEJA
, ConcreteEJA
, MatrixEJA
):
2173 sage: from mjo.eja.eja_algebra import (FiniteDimensionalEJA,
2174 ....: OctonionHermitianEJA)
2175 sage: from mjo.hurwitz import Octonions, OctonionMatrixAlgebra
2179 The 3-by-3 algebra satisfies the axioms of an EJA::
2181 sage: OctonionHermitianEJA(3, # long time
2182 ....: field=QQ, # long time
2183 ....: orthonormalize=False, # long time
2184 ....: check_axioms=True) # long time
2185 Euclidean Jordan algebra of dimension 27 over Rational Field
2187 After a change-of-basis, the 2-by-2 algebra has the same
2188 multiplication table as the ten-dimensional Jordan spin algebra::
2190 sage: A = OctonionMatrixAlgebra(2,Octonions(QQ),QQ)
2191 sage: b = OctonionHermitianEJA._denormalized_basis(A)
2192 sage: basis = (b[0] + b[9],) + b[1:9] + (b[0] - b[9],)
2193 sage: jp = OctonionHermitianEJA.jordan_product
2194 sage: ip = OctonionHermitianEJA.trace_inner_product
2195 sage: J = FiniteDimensionalEJA(basis,
2199 ....: orthonormalize=False)
2200 sage: J.multiplication_table()
2201 +----++----+----+----+----+----+----+----+----+----+----+
2202 | * || b0 | b1 | b2 | b3 | b4 | b5 | b6 | b7 | b8 | b9 |
2203 +====++====+====+====+====+====+====+====+====+====+====+
2204 | b0 || b0 | b1 | b2 | b3 | b4 | b5 | b6 | b7 | b8 | b9 |
2205 +----++----+----+----+----+----+----+----+----+----+----+
2206 | b1 || b1 | b0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
2207 +----++----+----+----+----+----+----+----+----+----+----+
2208 | b2 || b2 | 0 | b0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
2209 +----++----+----+----+----+----+----+----+----+----+----+
2210 | b3 || b3 | 0 | 0 | b0 | 0 | 0 | 0 | 0 | 0 | 0 |
2211 +----++----+----+----+----+----+----+----+----+----+----+
2212 | b4 || b4 | 0 | 0 | 0 | b0 | 0 | 0 | 0 | 0 | 0 |
2213 +----++----+----+----+----+----+----+----+----+----+----+
2214 | b5 || b5 | 0 | 0 | 0 | 0 | b0 | 0 | 0 | 0 | 0 |
2215 +----++----+----+----+----+----+----+----+----+----+----+
2216 | b6 || b6 | 0 | 0 | 0 | 0 | 0 | b0 | 0 | 0 | 0 |
2217 +----++----+----+----+----+----+----+----+----+----+----+
2218 | b7 || b7 | 0 | 0 | 0 | 0 | 0 | 0 | b0 | 0 | 0 |
2219 +----++----+----+----+----+----+----+----+----+----+----+
2220 | b8 || b8 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | b0 | 0 |
2221 +----++----+----+----+----+----+----+----+----+----+----+
2222 | b9 || b9 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | b0 |
2223 +----++----+----+----+----+----+----+----+----+----+----+
2227 We can actually construct the 27-dimensional Albert algebra,
2228 and we get the right unit element if we recompute it::
2230 sage: J = OctonionHermitianEJA(3, # long time
2231 ....: field=QQ, # long time
2232 ....: orthonormalize=False) # long time
2233 sage: J.one.clear_cache() # long time
2234 sage: J.one() # long time
2236 sage: J.one().to_matrix() # long time
2245 The 2-by-2 algebra is isomorphic to the ten-dimensional Jordan
2246 spin algebra, but just to be sure, we recompute its rank::
2248 sage: J = OctonionHermitianEJA(2, # long time
2249 ....: field=QQ, # long time
2250 ....: orthonormalize=False) # long time
2251 sage: J.rank.clear_cache() # long time
2252 sage: J.rank() # long time
2257 def _max_random_instance_size():
2259 The maximum rank of a random QuaternionHermitianEJA.
2261 return 1 # Dimension 1
2264 def random_instance(cls
, **kwargs
):
2266 Return a random instance of this type of algebra.
2268 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2269 return cls(n
, **kwargs
)
2271 def __init__(self
, n
, field
=AA
, **kwargs
):
2273 # Otherwise we don't get an EJA.
2274 raise ValueError("n cannot exceed 3")
2276 # We know this is a valid EJA, but will double-check
2277 # if the user passes check_axioms=True.
2278 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2280 from mjo
.hurwitz
import OctonionMatrixAlgebra
2281 A
= OctonionMatrixAlgebra(n
, scalars
=field
)
2282 super().__init
__(self
._denormalized
_basis
(A
),
2283 self
.jordan_product
,
2284 self
.trace_inner_product
,
2288 # TODO: this could be factored out somehow, but is left here
2289 # because the MatrixEJA is not presently a subclass of the
2290 # FDEJA class that defines rank() and one().
2291 self
.rank
.set_cache(n
)
2293 self
.one
.set_cache( self
.zero() )
2295 self
.one
.set_cache(self(A
.one()))
2298 class AlbertEJA(OctonionHermitianEJA
):
2300 The Albert algebra is the algebra of three-by-three Hermitian
2301 matrices whose entries are octonions.
2305 sage: from mjo.eja.eja_algebra import AlbertEJA
2309 sage: AlbertEJA(field=QQ, orthonormalize=False)
2310 Euclidean Jordan algebra of dimension 27 over Rational Field
2311 sage: AlbertEJA() # long time
2312 Euclidean Jordan algebra of dimension 27 over Algebraic Real Field
2315 def __init__(self
, *args
, **kwargs
):
2316 super().__init
__(3, *args
, **kwargs
)
2319 class HadamardEJA(RationalBasisEJA
, ConcreteEJA
):
2321 Return the Euclidean Jordan algebra on `R^n` with the Hadamard
2322 (pointwise real-number multiplication) Jordan product and the
2323 usual inner-product.
2325 This is nothing more than the Cartesian product of ``n`` copies of
2326 the one-dimensional Jordan spin algebra, and is the most common
2327 example of a non-simple Euclidean Jordan algebra.
2331 sage: from mjo.eja.eja_algebra import HadamardEJA
2335 This multiplication table can be verified by hand::
2337 sage: J = HadamardEJA(3)
2338 sage: b0,b1,b2 = J.gens()
2354 We can change the generator prefix::
2356 sage: HadamardEJA(3, prefix='r').gens()
2359 def __init__(self
, n
, field
=AA
, **kwargs
):
2361 jordan_product
= lambda x
,y
: x
2362 inner_product
= lambda x
,y
: x
2364 def jordan_product(x
,y
):
2366 return P( xi
*yi
for (xi
,yi
) in zip(x
,y
) )
2368 def inner_product(x
,y
):
2371 # New defaults for keyword arguments. Don't orthonormalize
2372 # because our basis is already orthonormal with respect to our
2373 # inner-product. Don't check the axioms, because we know this
2374 # is a valid EJA... but do double-check if the user passes
2375 # check_axioms=True. Note: we DON'T override the "check_field"
2376 # default here, because the user can pass in a field!
2377 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2378 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2380 column_basis
= tuple( b
.column()
2381 for b
in FreeModule(field
, n
).basis() )
2382 super().__init
__(column_basis
,
2388 self
.rank
.set_cache(n
)
2391 self
.one
.set_cache( self
.zero() )
2393 self
.one
.set_cache( sum(self
.gens()) )
2396 def _max_random_instance_size():
2398 The maximum dimension of a random HadamardEJA.
2403 def random_instance(cls
, **kwargs
):
2405 Return a random instance of this type of algebra.
2407 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2408 return cls(n
, **kwargs
)
2411 class BilinearFormEJA(RationalBasisEJA
, ConcreteEJA
):
2413 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2414 with the half-trace inner product and jordan product ``x*y =
2415 (<Bx,y>,y_bar>, x0*y_bar + y0*x_bar)`` where `B = 1 \times B22` is
2416 a symmetric positive-definite "bilinear form" matrix. Its
2417 dimension is the size of `B`, and it has rank two in dimensions
2418 larger than two. It reduces to the ``JordanSpinEJA`` when `B` is
2419 the identity matrix of order ``n``.
2421 We insist that the one-by-one upper-left identity block of `B` be
2422 passed in as well so that we can be passed a matrix of size zero
2423 to construct a trivial algebra.
2427 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
2428 ....: JordanSpinEJA)
2432 When no bilinear form is specified, the identity matrix is used,
2433 and the resulting algebra is the Jordan spin algebra::
2435 sage: B = matrix.identity(AA,3)
2436 sage: J0 = BilinearFormEJA(B)
2437 sage: J1 = JordanSpinEJA(3)
2438 sage: J0.multiplication_table() == J0.multiplication_table()
2441 An error is raised if the matrix `B` does not correspond to a
2442 positive-definite bilinear form::
2444 sage: B = matrix.random(QQ,2,3)
2445 sage: J = BilinearFormEJA(B)
2446 Traceback (most recent call last):
2448 ValueError: bilinear form is not positive-definite
2449 sage: B = matrix.zero(QQ,3)
2450 sage: J = BilinearFormEJA(B)
2451 Traceback (most recent call last):
2453 ValueError: bilinear form is not positive-definite
2457 We can create a zero-dimensional algebra::
2459 sage: B = matrix.identity(AA,0)
2460 sage: J = BilinearFormEJA(B)
2464 We can check the multiplication condition given in the Jordan, von
2465 Neumann, and Wigner paper (and also discussed on my "On the
2466 symmetry..." paper). Note that this relies heavily on the standard
2467 choice of basis, as does anything utilizing the bilinear form
2468 matrix. We opt not to orthonormalize the basis, because if we
2469 did, we would have to normalize the `s_{i}` in a similar manner::
2471 sage: set_random_seed()
2472 sage: n = ZZ.random_element(5)
2473 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2474 sage: B11 = matrix.identity(QQ,1)
2475 sage: B22 = M.transpose()*M
2476 sage: B = block_matrix(2,2,[ [B11,0 ],
2478 sage: J = BilinearFormEJA(B, orthonormalize=False)
2479 sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
2480 sage: V = J.vector_space()
2481 sage: sis = [ J( V([0] + (M.inverse()*ei).list()).column() )
2482 ....: for ei in eis ]
2483 sage: actual = [ sis[i]*sis[j]
2484 ....: for i in range(n-1)
2485 ....: for j in range(n-1) ]
2486 sage: expected = [ J.one() if i == j else J.zero()
2487 ....: for i in range(n-1)
2488 ....: for j in range(n-1) ]
2489 sage: actual == expected
2493 def __init__(self
, B
, field
=AA
, **kwargs
):
2494 # The matrix "B" is supplied by the user in most cases,
2495 # so it makes sense to check whether or not its positive-
2496 # definite unless we are specifically asked not to...
2497 if ("check_axioms" not in kwargs
) or kwargs
["check_axioms"]:
2498 if not B
.is_positive_definite():
2499 raise ValueError("bilinear form is not positive-definite")
2501 # However, all of the other data for this EJA is computed
2502 # by us in manner that guarantees the axioms are
2503 # satisfied. So, again, unless we are specifically asked to
2504 # verify things, we'll skip the rest of the checks.
2505 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2507 def inner_product(x
,y
):
2508 return (y
.T
*B
*x
)[0,0]
2510 def jordan_product(x
,y
):
2516 z0
= inner_product(y
,x
)
2517 zbar
= y0
*xbar
+ x0
*ybar
2518 return P([z0
] + zbar
.list())
2521 column_basis
= tuple( b
.column()
2522 for b
in FreeModule(field
, n
).basis() )
2524 # TODO: I haven't actually checked this, but it seems legit.
2529 super().__init
__(column_basis
,
2533 associative
=associative
,
2536 # The rank of this algebra is two, unless we're in a
2537 # one-dimensional ambient space (because the rank is bounded
2538 # by the ambient dimension).
2539 self
.rank
.set_cache(min(n
,2))
2542 self
.one
.set_cache( self
.zero() )
2544 self
.one
.set_cache( self
.monomial(0) )
2547 def _max_random_instance_size():
2549 The maximum dimension of a random BilinearFormEJA.
2554 def random_instance(cls
, **kwargs
):
2556 Return a random instance of this algebra.
2558 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2560 B
= matrix
.identity(ZZ
, n
)
2561 return cls(B
, **kwargs
)
2563 B11
= matrix
.identity(ZZ
, 1)
2564 M
= matrix
.random(ZZ
, n
-1)
2565 I
= matrix
.identity(ZZ
, n
-1)
2567 while alpha
.is_zero():
2568 alpha
= ZZ
.random_element().abs()
2569 B22
= M
.transpose()*M
+ alpha
*I
2571 from sage
.matrix
.special
import block_matrix
2572 B
= block_matrix(2,2, [ [B11
, ZZ(0) ],
2575 return cls(B
, **kwargs
)
2578 class JordanSpinEJA(BilinearFormEJA
):
2580 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2581 with the usual inner product and jordan product ``x*y =
2582 (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2587 sage: from mjo.eja.eja_algebra import JordanSpinEJA
2591 This multiplication table can be verified by hand::
2593 sage: J = JordanSpinEJA(4)
2594 sage: b0,b1,b2,b3 = J.gens()
2610 We can change the generator prefix::
2612 sage: JordanSpinEJA(2, prefix='B').gens()
2617 Ensure that we have the usual inner product on `R^n`::
2619 sage: set_random_seed()
2620 sage: J = JordanSpinEJA.random_instance()
2621 sage: x,y = J.random_elements(2)
2622 sage: actual = x.inner_product(y)
2623 sage: expected = x.to_vector().inner_product(y.to_vector())
2624 sage: actual == expected
2628 def __init__(self
, n
, *args
, **kwargs
):
2629 # This is a special case of the BilinearFormEJA with the
2630 # identity matrix as its bilinear form.
2631 B
= matrix
.identity(ZZ
, n
)
2633 # Don't orthonormalize because our basis is already
2634 # orthonormal with respect to our inner-product.
2635 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2637 # But also don't pass check_field=False here, because the user
2638 # can pass in a field!
2639 super().__init
__(B
, *args
, **kwargs
)
2642 def _max_random_instance_size():
2644 The maximum dimension of a random JordanSpinEJA.
2649 def random_instance(cls
, **kwargs
):
2651 Return a random instance of this type of algebra.
2653 Needed here to override the implementation for ``BilinearFormEJA``.
2655 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2656 return cls(n
, **kwargs
)
2659 class TrivialEJA(RationalBasisEJA
, ConcreteEJA
):
2661 The trivial Euclidean Jordan algebra consisting of only a zero element.
2665 sage: from mjo.eja.eja_algebra import TrivialEJA
2669 sage: J = TrivialEJA()
2676 sage: 7*J.one()*12*J.one()
2678 sage: J.one().inner_product(J.one())
2680 sage: J.one().norm()
2682 sage: J.one().subalgebra_generated_by()
2683 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2688 def __init__(self
, **kwargs
):
2689 jordan_product
= lambda x
,y
: x
2690 inner_product
= lambda x
,y
: 0
2693 # New defaults for keyword arguments
2694 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2695 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2697 super().__init
__(basis
,
2703 # The rank is zero using my definition, namely the dimension of the
2704 # largest subalgebra generated by any element.
2705 self
.rank
.set_cache(0)
2706 self
.one
.set_cache( self
.zero() )
2709 def random_instance(cls
, **kwargs
):
2710 # We don't take a "size" argument so the superclass method is
2711 # inappropriate for us.
2712 return cls(**kwargs
)
2715 class CartesianProductEJA(FiniteDimensionalEJA
):
2717 The external (orthogonal) direct sum of two or more Euclidean
2718 Jordan algebras. Every Euclidean Jordan algebra decomposes into an
2719 orthogonal direct sum of simple Euclidean Jordan algebras which is
2720 then isometric to a Cartesian product, so no generality is lost by
2721 providing only this construction.
2725 sage: from mjo.eja.eja_algebra import (random_eja,
2726 ....: CartesianProductEJA,
2728 ....: JordanSpinEJA,
2729 ....: RealSymmetricEJA)
2733 The Jordan product is inherited from our factors and implemented by
2734 our CombinatorialFreeModule Cartesian product superclass::
2736 sage: set_random_seed()
2737 sage: J1 = HadamardEJA(2)
2738 sage: J2 = RealSymmetricEJA(2)
2739 sage: J = cartesian_product([J1,J2])
2740 sage: x,y = J.random_elements(2)
2744 The ability to retrieve the original factors is implemented by our
2745 CombinatorialFreeModule Cartesian product superclass::
2747 sage: J1 = HadamardEJA(2, field=QQ)
2748 sage: J2 = JordanSpinEJA(3, field=QQ)
2749 sage: J = cartesian_product([J1,J2])
2750 sage: J.cartesian_factors()
2751 (Euclidean Jordan algebra of dimension 2 over Rational Field,
2752 Euclidean Jordan algebra of dimension 3 over Rational Field)
2754 You can provide more than two factors::
2756 sage: J1 = HadamardEJA(2)
2757 sage: J2 = JordanSpinEJA(3)
2758 sage: J3 = RealSymmetricEJA(3)
2759 sage: cartesian_product([J1,J2,J3])
2760 Euclidean Jordan algebra of dimension 2 over Algebraic Real
2761 Field (+) Euclidean Jordan algebra of dimension 3 over Algebraic
2762 Real Field (+) Euclidean Jordan algebra of dimension 6 over
2763 Algebraic Real Field
2765 Rank is additive on a Cartesian product::
2767 sage: J1 = HadamardEJA(1)
2768 sage: J2 = RealSymmetricEJA(2)
2769 sage: J = cartesian_product([J1,J2])
2770 sage: J1.rank.clear_cache()
2771 sage: J2.rank.clear_cache()
2772 sage: J.rank.clear_cache()
2775 sage: J.rank() == J1.rank() + J2.rank()
2778 The same rank computation works over the rationals, with whatever
2781 sage: J1 = HadamardEJA(1, field=QQ, orthonormalize=False)
2782 sage: J2 = RealSymmetricEJA(2, field=QQ, orthonormalize=False)
2783 sage: J = cartesian_product([J1,J2])
2784 sage: J1.rank.clear_cache()
2785 sage: J2.rank.clear_cache()
2786 sage: J.rank.clear_cache()
2789 sage: J.rank() == J1.rank() + J2.rank()
2792 The product algebra will be associative if and only if all of its
2793 components are associative::
2795 sage: J1 = HadamardEJA(2)
2796 sage: J1.is_associative()
2798 sage: J2 = HadamardEJA(3)
2799 sage: J2.is_associative()
2801 sage: J3 = RealSymmetricEJA(3)
2802 sage: J3.is_associative()
2804 sage: CP1 = cartesian_product([J1,J2])
2805 sage: CP1.is_associative()
2807 sage: CP2 = cartesian_product([J1,J3])
2808 sage: CP2.is_associative()
2811 Cartesian products of Cartesian products work::
2813 sage: J1 = JordanSpinEJA(1)
2814 sage: J2 = JordanSpinEJA(1)
2815 sage: J3 = JordanSpinEJA(1)
2816 sage: J = cartesian_product([J1,cartesian_product([J2,J3])])
2817 sage: J.multiplication_table()
2818 +----++----+----+----+
2819 | * || b0 | b1 | b2 |
2820 +====++====+====+====+
2821 | b0 || b0 | 0 | 0 |
2822 +----++----+----+----+
2823 | b1 || 0 | b1 | 0 |
2824 +----++----+----+----+
2825 | b2 || 0 | 0 | b2 |
2826 +----++----+----+----+
2827 sage: HadamardEJA(3).multiplication_table()
2828 +----++----+----+----+
2829 | * || b0 | b1 | b2 |
2830 +====++====+====+====+
2831 | b0 || b0 | 0 | 0 |
2832 +----++----+----+----+
2833 | b1 || 0 | b1 | 0 |
2834 +----++----+----+----+
2835 | b2 || 0 | 0 | b2 |
2836 +----++----+----+----+
2840 All factors must share the same base field::
2842 sage: J1 = HadamardEJA(2, field=QQ)
2843 sage: J2 = RealSymmetricEJA(2)
2844 sage: CartesianProductEJA((J1,J2))
2845 Traceback (most recent call last):
2847 ValueError: all factors must share the same base field
2849 The cached unit element is the same one that would be computed::
2851 sage: set_random_seed() # long time
2852 sage: J1 = random_eja() # long time
2853 sage: J2 = random_eja() # long time
2854 sage: J = cartesian_product([J1,J2]) # long time
2855 sage: actual = J.one() # long time
2856 sage: J.one.clear_cache() # long time
2857 sage: expected = J.one() # long time
2858 sage: actual == expected # long time
2862 Element
= FiniteDimensionalEJAElement
2865 def __init__(self
, factors
, **kwargs
):
2870 self
._sets
= factors
2872 field
= factors
[0].base_ring()
2873 if not all( J
.base_ring() == field
for J
in factors
):
2874 raise ValueError("all factors must share the same base field")
2876 associative
= all( f
.is_associative() for f
in factors
)
2878 MS
= self
.matrix_space()
2882 for b
in factors
[i
].matrix_basis():
2887 basis
= tuple( MS(b
) for b
in basis
)
2889 # Define jordan/inner products that operate on that matrix_basis.
2890 def jordan_product(x
,y
):
2892 (factors
[i
](x
[i
])*factors
[i
](y
[i
])).to_matrix()
2896 def inner_product(x
, y
):
2898 factors
[i
](x
[i
]).inner_product(factors
[i
](y
[i
]))
2902 # There's no need to check the field since it already came
2903 # from an EJA. Likewise the axioms are guaranteed to be
2904 # satisfied, unless the guy writing this class sucks.
2906 # If you want the basis to be orthonormalized, orthonormalize
2908 FiniteDimensionalEJA
.__init
__(self
,
2913 orthonormalize
=False,
2914 associative
=associative
,
2915 cartesian_product
=True,
2919 ones
= tuple(J
.one().to_matrix() for J
in factors
)
2920 self
.one
.set_cache(self(ones
))
2921 self
.rank
.set_cache(sum(J
.rank() for J
in factors
))
2923 def cartesian_factors(self
):
2924 # Copy/pasted from CombinatorialFreeModule_CartesianProduct.
2927 def cartesian_factor(self
, i
):
2929 Return the ``i``th factor of this algebra.
2931 return self
._sets
[i
]
2934 # Copy/pasted from CombinatorialFreeModule_CartesianProduct.
2935 from sage
.categories
.cartesian_product
import cartesian_product
2936 return cartesian_product
.symbol
.join("%s" % factor
2937 for factor
in self
._sets
)
2939 def matrix_space(self
):
2941 Return the space that our matrix basis lives in as a Cartesian
2944 We don't simply use the ``cartesian_product()`` functor here
2945 because it acts differently on SageMath MatrixSpaces and our
2946 custom MatrixAlgebras, which are CombinatorialFreeModules. We
2947 always want the result to be represented (and indexed) as
2952 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
2954 ....: OctonionHermitianEJA,
2955 ....: RealSymmetricEJA)
2959 sage: J1 = HadamardEJA(1)
2960 sage: J2 = RealSymmetricEJA(2)
2961 sage: J = cartesian_product([J1,J2])
2962 sage: J.matrix_space()
2963 The Cartesian product of (Full MatrixSpace of 1 by 1 dense
2964 matrices over Algebraic Real Field, Full MatrixSpace of 2
2965 by 2 dense matrices over Algebraic Real Field)
2969 sage: J1 = ComplexHermitianEJA(1)
2970 sage: J2 = ComplexHermitianEJA(1)
2971 sage: J = cartesian_product([J1,J2])
2972 sage: J.one().to_matrix()[0]
2976 sage: J.one().to_matrix()[1]
2983 sage: J1 = OctonionHermitianEJA(1)
2984 sage: J2 = OctonionHermitianEJA(1)
2985 sage: J = cartesian_product([J1,J2])
2986 sage: J.one().to_matrix()[0]
2990 sage: J.one().to_matrix()[1]
2996 scalars
= self
.cartesian_factor(0).base_ring()
2998 # This category isn't perfect, but is good enough for what we
3000 cat
= MagmaticAlgebras(scalars
).FiniteDimensional().WithBasis()
3001 cat
= cat
.Unital().CartesianProducts()
3002 factors
= tuple( J
.matrix_space() for J
in self
.cartesian_factors() )
3004 from sage
.sets
.cartesian_product
import CartesianProduct
3005 return CartesianProduct(factors
, cat
)
3009 def cartesian_projection(self
, i
):
3013 sage: from mjo.eja.eja_algebra import (random_eja,
3014 ....: JordanSpinEJA,
3016 ....: RealSymmetricEJA,
3017 ....: ComplexHermitianEJA)
3021 The projection morphisms are Euclidean Jordan algebra
3024 sage: J1 = HadamardEJA(2)
3025 sage: J2 = RealSymmetricEJA(2)
3026 sage: J = cartesian_product([J1,J2])
3027 sage: J.cartesian_projection(0)
3028 Linear operator between finite-dimensional Euclidean Jordan
3029 algebras represented by the matrix:
3032 Domain: Euclidean Jordan algebra of dimension 2 over Algebraic
3033 Real Field (+) Euclidean Jordan algebra of dimension 3 over
3034 Algebraic Real Field
3035 Codomain: Euclidean Jordan algebra of dimension 2 over Algebraic
3037 sage: J.cartesian_projection(1)
3038 Linear operator between finite-dimensional Euclidean Jordan
3039 algebras represented by the matrix:
3043 Domain: Euclidean Jordan algebra of dimension 2 over Algebraic
3044 Real Field (+) Euclidean Jordan algebra of dimension 3 over
3045 Algebraic Real Field
3046 Codomain: Euclidean Jordan algebra of dimension 3 over Algebraic
3049 The projections work the way you'd expect on the vector
3050 representation of an element::
3052 sage: J1 = JordanSpinEJA(2)
3053 sage: J2 = ComplexHermitianEJA(2)
3054 sage: J = cartesian_product([J1,J2])
3055 sage: pi_left = J.cartesian_projection(0)
3056 sage: pi_right = J.cartesian_projection(1)
3057 sage: pi_left(J.one()).to_vector()
3059 sage: pi_right(J.one()).to_vector()
3061 sage: J.one().to_vector()
3066 The answer never changes::
3068 sage: set_random_seed()
3069 sage: J1 = random_eja()
3070 sage: J2 = random_eja()
3071 sage: J = cartesian_product([J1,J2])
3072 sage: P0 = J.cartesian_projection(0)
3073 sage: P1 = J.cartesian_projection(0)
3078 offset
= sum( self
.cartesian_factor(k
).dimension()
3080 Ji
= self
.cartesian_factor(i
)
3081 Pi
= self
._module
_morphism
(lambda j
: Ji
.monomial(j
- offset
),
3084 return FiniteDimensionalEJAOperator(self
,Ji
,Pi
.matrix())
3087 def cartesian_embedding(self
, i
):
3091 sage: from mjo.eja.eja_algebra import (random_eja,
3092 ....: JordanSpinEJA,
3094 ....: RealSymmetricEJA)
3098 The embedding morphisms are Euclidean Jordan algebra
3101 sage: J1 = HadamardEJA(2)
3102 sage: J2 = RealSymmetricEJA(2)
3103 sage: J = cartesian_product([J1,J2])
3104 sage: J.cartesian_embedding(0)
3105 Linear operator between finite-dimensional Euclidean Jordan
3106 algebras represented by the matrix:
3112 Domain: Euclidean Jordan algebra of dimension 2 over
3113 Algebraic Real Field
3114 Codomain: Euclidean Jordan algebra of dimension 2 over
3115 Algebraic Real Field (+) Euclidean Jordan algebra of
3116 dimension 3 over Algebraic Real Field
3117 sage: J.cartesian_embedding(1)
3118 Linear operator between finite-dimensional Euclidean Jordan
3119 algebras represented by the matrix:
3125 Domain: Euclidean Jordan algebra of dimension 3 over
3126 Algebraic Real Field
3127 Codomain: Euclidean Jordan algebra of dimension 2 over
3128 Algebraic Real Field (+) Euclidean Jordan algebra of
3129 dimension 3 over Algebraic Real Field
3131 The embeddings work the way you'd expect on the vector
3132 representation of an element::
3134 sage: J1 = JordanSpinEJA(3)
3135 sage: J2 = RealSymmetricEJA(2)
3136 sage: J = cartesian_product([J1,J2])
3137 sage: iota_left = J.cartesian_embedding(0)
3138 sage: iota_right = J.cartesian_embedding(1)
3139 sage: iota_left(J1.zero()) == J.zero()
3141 sage: iota_right(J2.zero()) == J.zero()
3143 sage: J1.one().to_vector()
3145 sage: iota_left(J1.one()).to_vector()
3147 sage: J2.one().to_vector()
3149 sage: iota_right(J2.one()).to_vector()
3151 sage: J.one().to_vector()
3156 The answer never changes::
3158 sage: set_random_seed()
3159 sage: J1 = random_eja()
3160 sage: J2 = random_eja()
3161 sage: J = cartesian_product([J1,J2])
3162 sage: E0 = J.cartesian_embedding(0)
3163 sage: E1 = J.cartesian_embedding(0)
3167 Composing a projection with the corresponding inclusion should
3168 produce the identity map, and mismatching them should produce
3171 sage: set_random_seed()
3172 sage: J1 = random_eja()
3173 sage: J2 = random_eja()
3174 sage: J = cartesian_product([J1,J2])
3175 sage: iota_left = J.cartesian_embedding(0)
3176 sage: iota_right = J.cartesian_embedding(1)
3177 sage: pi_left = J.cartesian_projection(0)
3178 sage: pi_right = J.cartesian_projection(1)
3179 sage: pi_left*iota_left == J1.one().operator()
3181 sage: pi_right*iota_right == J2.one().operator()
3183 sage: (pi_left*iota_right).is_zero()
3185 sage: (pi_right*iota_left).is_zero()
3189 offset
= sum( self
.cartesian_factor(k
).dimension()
3191 Ji
= self
.cartesian_factor(i
)
3192 Ei
= Ji
._module
_morphism
(lambda j
: self
.monomial(j
+ offset
),
3194 return FiniteDimensionalEJAOperator(Ji
,self
,Ei
.matrix())
3198 FiniteDimensionalEJA
.CartesianProduct
= CartesianProductEJA
3200 class RationalBasisCartesianProductEJA(CartesianProductEJA
,
3203 A separate class for products of algebras for which we know a
3208 sage: from mjo.eja.eja_algebra import (HadamardEJA,
3209 ....: JordanSpinEJA,
3210 ....: OctonionHermitianEJA,
3211 ....: RealSymmetricEJA)
3215 This gives us fast characteristic polynomial computations in
3216 product algebras, too::
3219 sage: J1 = JordanSpinEJA(2)
3220 sage: J2 = RealSymmetricEJA(3)
3221 sage: J = cartesian_product([J1,J2])
3222 sage: J.characteristic_polynomial_of().degree()
3229 The ``cartesian_product()`` function only uses the first factor to
3230 decide where the result will live; thus we have to be careful to
3231 check that all factors do indeed have a `_rational_algebra` member
3232 before we try to access it::
3234 sage: J1 = OctonionHermitianEJA(1) # no rational basis
3235 sage: J2 = HadamardEJA(2)
3236 sage: cartesian_product([J1,J2])
3237 Euclidean Jordan algebra of dimension 1 over Algebraic Real Field
3238 (+) Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
3239 sage: cartesian_product([J2,J1])
3240 Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
3241 (+) Euclidean Jordan algebra of dimension 1 over Algebraic Real Field
3244 def __init__(self
, algebras
, **kwargs
):
3245 CartesianProductEJA
.__init
__(self
, algebras
, **kwargs
)
3247 self
._rational
_algebra
= None
3248 if self
.vector_space().base_field() is not QQ
:
3249 if all( hasattr(r
, "_rational_algebra") for r
in algebras
):
3250 self
._rational
_algebra
= cartesian_product([
3251 r
._rational
_algebra
for r
in algebras
3255 RationalBasisEJA
.CartesianProduct
= RationalBasisCartesianProductEJA
3257 def random_eja(*args
, **kwargs
):
3258 J1
= ConcreteEJA
.random_instance(*args
, **kwargs
)
3260 # This might make Cartesian products appear roughly as often as
3261 # any other ConcreteEJA.
3262 if ZZ
.random_element(len(ConcreteEJA
.__subclasses
__()) + 1) == 0:
3263 # Use random_eja() again so we can get more than two factors.
3264 J2
= random_eja(*args
, **kwargs
)
3265 J
= cartesian_product([J1
,J2
])