2 Representations and constructions for Euclidean Jordan algebras.
4 A Euclidean Jordan algebra is a Jordan algebra that has some
7 1. It is finite-dimensional.
8 2. Its scalar field is the real numbers.
9 3a. An inner product is defined on it, and...
10 3b. That inner product is compatible with the Jordan product
11 in the sense that `<x*y,z> = <y,x*z>` for all elements
12 `x,y,z` in the algebra.
14 Every Euclidean Jordan algebra is formally-real: for any two elements
15 `x` and `y` in the algebra, `x^{2} + y^{2} = 0` implies that `x = y =
16 0`. Conversely, every finite-dimensional formally-real Jordan algebra
17 can be made into a Euclidean Jordan algebra with an appropriate choice
20 Formally-real Jordan algebras were originally studied as a framework
21 for quantum mechanics. Today, Euclidean Jordan algebras are crucial in
22 symmetric cone optimization, since every symmetric cone arises as the
23 cone of squares in some Euclidean Jordan algebra.
25 It is known that every Euclidean Jordan algebra decomposes into an
26 orthogonal direct sum (essentially, a Cartesian product) of simple
27 algebras, and that moreover, up to Jordan-algebra isomorphism, there
28 are only five families of simple algebras. We provide constructions
29 for these simple algebras:
31 * :class:`BilinearFormEJA`
32 * :class:`RealSymmetricEJA`
33 * :class:`ComplexHermitianEJA`
34 * :class:`QuaternionHermitianEJA`
36 Missing from this list is the algebra of three-by-three octononion
37 Hermitian matrices, as there is (as of yet) no implementation of the
38 octonions in SageMath. In addition to these, we provide two other
39 example constructions,
41 * :class:`HadamardEJA`
44 The Jordan spin algebra is a bilinear form algebra where the bilinear
45 form is the identity. The Hadamard EJA is simply a Cartesian product
46 of one-dimensional spin algebras. And last but not least, the trivial
47 EJA is exactly what you think. Cartesian products of these are also
48 supported using the usual ``cartesian_product()`` function; as a
49 result, we support (up to isomorphism) all Euclidean Jordan algebras
50 that don't involve octonions.
54 sage: from mjo.eja.eja_algebra import random_eja
59 Euclidean Jordan algebra of dimension...
62 from itertools
import repeat
64 from sage
.algebras
.quatalg
.quaternion_algebra
import QuaternionAlgebra
65 from sage
.categories
.magmatic_algebras
import MagmaticAlgebras
66 from sage
.categories
.sets_cat
import cartesian_product
67 from sage
.combinat
.free_module
import (CombinatorialFreeModule
,
68 CombinatorialFreeModule_CartesianProduct
)
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.
116 Element
= FiniteDimensionalEJAElement
125 cartesian_product
=False,
130 # Keep track of whether or not the matrix basis consists of
131 # tuples, since we need special cases for them damned near
132 # everywhere. This is INDEPENDENT of whether or not the
133 # algebra is a cartesian product, since a subalgebra of a
134 # cartesian product will have a basis of tuples, but will not
135 # in general itself be a cartesian product algebra.
136 self
._matrix
_basis
_is
_cartesian
= False
139 if hasattr(basis
[0], 'cartesian_factors'):
140 self
._matrix
_basis
_is
_cartesian
= True
143 if not field
.is_subring(RR
):
144 # Note: this does return true for the real algebraic
145 # field, the rationals, and any quadratic field where
146 # we've specified a real embedding.
147 raise ValueError("scalar field is not real")
149 # If the basis given to us wasn't over the field that it's
150 # supposed to be over, fix that. Or, you know, crash.
151 if not cartesian_product
:
152 # The field for a cartesian product algebra comes from one
153 # of its factors and is the same for all factors, so
154 # there's no need to "reapply" it on product algebras.
155 if self
._matrix
_basis
_is
_cartesian
:
156 # OK since if n == 0, the basis does not consist of tuples.
157 P
= basis
[0].parent()
158 basis
= tuple( P(tuple(b_i
.change_ring(field
) for b_i
in b
))
161 basis
= tuple( b
.change_ring(field
) for b
in basis
)
165 # Check commutativity of the Jordan and inner-products.
166 # This has to be done before we build the multiplication
167 # and inner-product tables/matrices, because we take
168 # advantage of symmetry in the process.
169 if not all( jordan_product(bi
,bj
) == jordan_product(bj
,bi
)
172 raise ValueError("Jordan product is not commutative")
174 if not all( inner_product(bi
,bj
) == inner_product(bj
,bi
)
177 raise ValueError("inner-product is not commutative")
180 category
= MagmaticAlgebras(field
).FiniteDimensional()
181 category
= category
.WithBasis().Unital()
183 # Element subalgebras can take advantage of this.
184 category
= category
.Associative()
185 if cartesian_product
:
186 category
= category
.CartesianProducts()
188 # Call the superclass constructor so that we can use its from_vector()
189 # method to build our multiplication table.
190 CombinatorialFreeModule
.__init
__(self
,
197 # Now comes all of the hard work. We'll be constructing an
198 # ambient vector space V that our (vectorized) basis lives in,
199 # as well as a subspace W of V spanned by those (vectorized)
200 # basis elements. The W-coordinates are the coefficients that
201 # we see in things like x = 1*e1 + 2*e2.
206 degree
= len(_all2list(basis
[0]))
208 # Build an ambient space that fits our matrix basis when
209 # written out as "long vectors."
210 V
= VectorSpace(field
, degree
)
212 # The matrix that will hole the orthonormal -> unorthonormal
213 # coordinate transformation.
214 self
._deortho
_matrix
= None
217 # Save a copy of the un-orthonormalized basis for later.
218 # Convert it to ambient V (vector) coordinates while we're
219 # at it, because we'd have to do it later anyway.
220 deortho_vector_basis
= tuple( V(_all2list(b
)) for b
in basis
)
222 from mjo
.eja
.eja_utils
import gram_schmidt
223 basis
= tuple(gram_schmidt(basis
, inner_product
))
225 # Save the (possibly orthonormalized) matrix basis for
227 self
._matrix
_basis
= basis
229 # Now create the vector space for the algebra, which will have
230 # its own set of non-ambient coordinates (in terms of the
232 vector_basis
= tuple( V(_all2list(b
)) for b
in basis
)
233 W
= V
.span_of_basis( vector_basis
, check
=check_axioms
)
236 # Now "W" is the vector space of our algebra coordinates. The
237 # variables "X1", "X2",... refer to the entries of vectors in
238 # W. Thus to convert back and forth between the orthonormal
239 # coordinates and the given ones, we need to stick the original
241 U
= V
.span_of_basis( deortho_vector_basis
, check
=check_axioms
)
242 self
._deortho
_matrix
= matrix( U
.coordinate_vector(q
)
243 for q
in vector_basis
)
246 # Now we actually compute the multiplication and inner-product
247 # tables/matrices using the possibly-orthonormalized basis.
248 self
._inner
_product
_matrix
= matrix
.identity(field
, n
)
249 self
._multiplication
_table
= [ [0 for j
in range(i
+1)]
252 # Note: the Jordan and inner-products are defined in terms
253 # of the ambient basis. It's important that their arguments
254 # are in ambient coordinates as well.
257 # ortho basis w.r.t. ambient coords
261 # The jordan product returns a matrixy answer, so we
262 # have to convert it to the algebra coordinates.
263 elt
= jordan_product(q_i
, q_j
)
264 elt
= W
.coordinate_vector(V(_all2list(elt
)))
265 self
._multiplication
_table
[i
][j
] = self
.from_vector(elt
)
267 if not orthonormalize
:
268 # If we're orthonormalizing the basis with respect
269 # to an inner-product, then the inner-product
270 # matrix with respect to the resulting basis is
271 # just going to be the identity.
272 ip
= inner_product(q_i
, q_j
)
273 self
._inner
_product
_matrix
[i
,j
] = ip
274 self
._inner
_product
_matrix
[j
,i
] = ip
276 self
._inner
_product
_matrix
._cache
= {'hermitian': True}
277 self
._inner
_product
_matrix
.set_immutable()
280 if not self
._is
_jordanian
():
281 raise ValueError("Jordan identity does not hold")
282 if not self
._inner
_product
_is
_associative
():
283 raise ValueError("inner product is not associative")
286 def _coerce_map_from_base_ring(self
):
288 Disable the map from the base ring into the algebra.
290 Performing a nonsense conversion like this automatically
291 is counterpedagogical. The fallback is to try the usual
292 element constructor, which should also fail.
296 sage: from mjo.eja.eja_algebra import random_eja
300 sage: set_random_seed()
301 sage: J = random_eja()
303 Traceback (most recent call last):
305 ValueError: not an element of this algebra
311 def product_on_basis(self
, i
, j
):
313 Returns the Jordan product of the `i` and `j`th basis elements.
315 This completely defines the Jordan product on the algebra, and
316 is used direclty by our superclass machinery to implement
321 sage: from mjo.eja.eja_algebra import random_eja
325 sage: set_random_seed()
326 sage: J = random_eja()
327 sage: n = J.dimension()
330 sage: ei_ej = J.zero()*J.zero()
332 ....: i = ZZ.random_element(n)
333 ....: j = ZZ.random_element(n)
334 ....: ei = J.gens()[i]
335 ....: ej = J.gens()[j]
336 ....: ei_ej = J.product_on_basis(i,j)
341 # We only stored the lower-triangular portion of the
342 # multiplication table.
344 return self
._multiplication
_table
[i
][j
]
346 return self
._multiplication
_table
[j
][i
]
348 def inner_product(self
, x
, y
):
350 The inner product associated with this Euclidean Jordan algebra.
352 Defaults to the trace inner product, but can be overridden by
353 subclasses if they are sure that the necessary properties are
358 sage: from mjo.eja.eja_algebra import (random_eja,
360 ....: BilinearFormEJA)
364 Our inner product is "associative," which means the following for
365 a symmetric bilinear form::
367 sage: set_random_seed()
368 sage: J = random_eja()
369 sage: x,y,z = J.random_elements(3)
370 sage: (x*y).inner_product(z) == y.inner_product(x*z)
375 Ensure that this is the usual inner product for the algebras
378 sage: set_random_seed()
379 sage: J = HadamardEJA.random_instance()
380 sage: x,y = J.random_elements(2)
381 sage: actual = x.inner_product(y)
382 sage: expected = x.to_vector().inner_product(y.to_vector())
383 sage: actual == expected
386 Ensure that this is one-half of the trace inner-product in a
387 BilinearFormEJA that isn't just the reals (when ``n`` isn't
388 one). This is in Faraut and Koranyi, and also my "On the
391 sage: set_random_seed()
392 sage: J = BilinearFormEJA.random_instance()
393 sage: n = J.dimension()
394 sage: x = J.random_element()
395 sage: y = J.random_element()
396 sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
400 B
= self
._inner
_product
_matrix
401 return (B
*x
.to_vector()).inner_product(y
.to_vector())
404 def is_associative(self
):
406 Return whether or not this algebra's Jordan product is associative.
410 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
414 sage: J = ComplexHermitianEJA(3, field=QQ, orthonormalize=False)
415 sage: J.is_associative()
417 sage: x = sum(J.gens())
418 sage: A = x.subalgebra_generated_by(orthonormalize=False)
419 sage: A.is_associative()
423 return "Associative" in self
.category().axioms()
425 def _is_jordanian(self
):
427 Whether or not this algebra's multiplication table respects the
428 Jordan identity `(x^{2})(xy) = x(x^{2}y)`.
430 We only check one arrangement of `x` and `y`, so for a
431 ``True`` result to be truly true, you should also check
432 :meth:`is_commutative`. This method should of course always
433 return ``True``, unless this algebra was constructed with
434 ``check_axioms=False`` and passed an invalid multiplication table.
436 return all( (self
.gens()[i
]**2)*(self
.gens()[i
]*self
.gens()[j
])
438 (self
.gens()[i
])*((self
.gens()[i
]**2)*self
.gens()[j
])
439 for i
in range(self
.dimension())
440 for j
in range(self
.dimension()) )
442 def _inner_product_is_associative(self
):
444 Return whether or not this algebra's inner product `B` is
445 associative; that is, whether or not `B(xy,z) = B(x,yz)`.
447 This method should of course always return ``True``, unless
448 this algebra was constructed with ``check_axioms=False`` and
449 passed an invalid Jordan or inner-product.
453 # Used to check whether or not something is zero.
456 # This choice is sufficient to allow the construction of
457 # QuaternionHermitianEJA(2, field=RDF) with check_axioms=True.
460 for i
in range(self
.dimension()):
461 for j
in range(self
.dimension()):
462 for k
in range(self
.dimension()):
466 diff
= (x
*y
).inner_product(z
) - x
.inner_product(y
*z
)
468 if diff
.abs() > epsilon
:
473 def _element_constructor_(self
, elt
):
475 Construct an element of this algebra from its vector or matrix
478 This gets called only after the parent element _call_ method
479 fails to find a coercion for the argument.
483 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
485 ....: RealSymmetricEJA)
489 The identity in `S^n` is converted to the identity in the EJA::
491 sage: J = RealSymmetricEJA(3)
492 sage: I = matrix.identity(QQ,3)
493 sage: J(I) == J.one()
496 This skew-symmetric matrix can't be represented in the EJA::
498 sage: J = RealSymmetricEJA(3)
499 sage: A = matrix(QQ,3, lambda i,j: i-j)
501 Traceback (most recent call last):
503 ValueError: not an element of this algebra
505 Tuples work as well, provided that the matrix basis for the
506 algebra consists of them::
508 sage: J1 = HadamardEJA(3)
509 sage: J2 = RealSymmetricEJA(2)
510 sage: J = cartesian_product([J1,J2])
511 sage: J( (J1.matrix_basis()[1], J2.matrix_basis()[2]) )
516 Ensure that we can convert any element of the two non-matrix
517 simple algebras (whose matrix representations are columns)
518 back and forth faithfully::
520 sage: set_random_seed()
521 sage: J = HadamardEJA.random_instance()
522 sage: x = J.random_element()
523 sage: J(x.to_vector().column()) == x
525 sage: J = JordanSpinEJA.random_instance()
526 sage: x = J.random_element()
527 sage: J(x.to_vector().column()) == x
530 We cannot coerce elements between algebras just because their
531 matrix representations are compatible::
533 sage: J1 = HadamardEJA(3)
534 sage: J2 = JordanSpinEJA(3)
536 Traceback (most recent call last):
538 ValueError: not an element of this algebra
540 Traceback (most recent call last):
542 ValueError: not an element of this algebra
545 msg
= "not an element of this algebra"
546 if elt
in self
.base_ring():
547 # Ensure that no base ring -> algebra coercion is performed
548 # by this method. There's some stupidity in sage that would
549 # otherwise propagate to this method; for example, sage thinks
550 # that the integer 3 belongs to the space of 2-by-2 matrices.
551 raise ValueError(msg
)
554 # Try to convert a vector into a column-matrix...
556 except (AttributeError, TypeError):
557 # and ignore failure, because we weren't really expecting
558 # a vector as an argument anyway.
561 if elt
not in self
.matrix_space():
562 raise ValueError(msg
)
564 # Thanks for nothing! Matrix spaces aren't vector spaces in
565 # Sage, so we have to figure out its matrix-basis coordinates
566 # ourselves. We use the basis space's ring instead of the
567 # element's ring because the basis space might be an algebraic
568 # closure whereas the base ring of the 3-by-3 identity matrix
569 # could be QQ instead of QQbar.
571 # And, we also have to handle Cartesian product bases (when
572 # the matrix basis consists of tuples) here. The "good news"
573 # is that we're already converting everything to long vectors,
574 # and that strategy works for tuples as well.
576 # We pass check=False because the matrix basis is "guaranteed"
577 # to be linearly independent... right? Ha ha.
579 V
= VectorSpace(self
.base_ring(), len(elt
))
580 W
= V
.span_of_basis( (V(_all2list(s
)) for s
in self
.matrix_basis()),
584 coords
= W
.coordinate_vector(V(elt
))
585 except ArithmeticError: # vector is not in free module
586 raise ValueError(msg
)
588 return self
.from_vector(coords
)
592 Return a string representation of ``self``.
596 sage: from mjo.eja.eja_algebra import JordanSpinEJA
600 Ensure that it says what we think it says::
602 sage: JordanSpinEJA(2, field=AA)
603 Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
604 sage: JordanSpinEJA(3, field=RDF)
605 Euclidean Jordan algebra of dimension 3 over Real Double Field
608 fmt
= "Euclidean Jordan algebra of dimension {} over {}"
609 return fmt
.format(self
.dimension(), self
.base_ring())
613 def characteristic_polynomial_of(self
):
615 Return the algebra's "characteristic polynomial of" function,
616 which is itself a multivariate polynomial that, when evaluated
617 at the coordinates of some algebra element, returns that
618 element's characteristic polynomial.
620 The resulting polynomial has `n+1` variables, where `n` is the
621 dimension of this algebra. The first `n` variables correspond to
622 the coordinates of an algebra element: when evaluated at the
623 coordinates of an algebra element with respect to a certain
624 basis, the result is a univariate polynomial (in the one
625 remaining variable ``t``), namely the characteristic polynomial
630 sage: from mjo.eja.eja_algebra import JordanSpinEJA, TrivialEJA
634 The characteristic polynomial in the spin algebra is given in
635 Alizadeh, Example 11.11::
637 sage: J = JordanSpinEJA(3)
638 sage: p = J.characteristic_polynomial_of(); p
639 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
640 sage: xvec = J.one().to_vector()
644 By definition, the characteristic polynomial is a monic
645 degree-zero polynomial in a rank-zero algebra. Note that
646 Cayley-Hamilton is indeed satisfied since the polynomial
647 ``1`` evaluates to the identity element of the algebra on
650 sage: J = TrivialEJA()
651 sage: J.characteristic_polynomial_of()
658 # The list of coefficient polynomials a_0, a_1, a_2, ..., a_(r-1).
659 a
= self
._charpoly
_coefficients
()
661 # We go to a bit of trouble here to reorder the
662 # indeterminates, so that it's easier to evaluate the
663 # characteristic polynomial at x's coordinates and get back
664 # something in terms of t, which is what we want.
665 S
= PolynomialRing(self
.base_ring(),'t')
669 S
= PolynomialRing(S
, R
.variable_names())
672 return (t
**r
+ sum( a
[k
]*(t
**k
) for k
in range(r
) ))
674 def coordinate_polynomial_ring(self
):
676 The multivariate polynomial ring in which this algebra's
677 :meth:`characteristic_polynomial_of` lives.
681 sage: from mjo.eja.eja_algebra import (HadamardEJA,
682 ....: RealSymmetricEJA)
686 sage: J = HadamardEJA(2)
687 sage: J.coordinate_polynomial_ring()
688 Multivariate Polynomial Ring in X1, X2...
689 sage: J = RealSymmetricEJA(3,field=QQ,orthonormalize=False)
690 sage: J.coordinate_polynomial_ring()
691 Multivariate Polynomial Ring in X1, X2, X3, X4, X5, X6...
694 var_names
= tuple( "X%d" % z
for z
in range(1, self
.dimension()+1) )
695 return PolynomialRing(self
.base_ring(), var_names
)
697 def inner_product(self
, x
, y
):
699 The inner product associated with this Euclidean Jordan algebra.
701 Defaults to the trace inner product, but can be overridden by
702 subclasses if they are sure that the necessary properties are
707 sage: from mjo.eja.eja_algebra import (random_eja,
709 ....: BilinearFormEJA)
713 Our inner product is "associative," which means the following for
714 a symmetric bilinear form::
716 sage: set_random_seed()
717 sage: J = random_eja()
718 sage: x,y,z = J.random_elements(3)
719 sage: (x*y).inner_product(z) == y.inner_product(x*z)
724 Ensure that this is the usual inner product for the algebras
727 sage: set_random_seed()
728 sage: J = HadamardEJA.random_instance()
729 sage: x,y = J.random_elements(2)
730 sage: actual = x.inner_product(y)
731 sage: expected = x.to_vector().inner_product(y.to_vector())
732 sage: actual == expected
735 Ensure that this is one-half of the trace inner-product in a
736 BilinearFormEJA that isn't just the reals (when ``n`` isn't
737 one). This is in Faraut and Koranyi, and also my "On the
740 sage: set_random_seed()
741 sage: J = BilinearFormEJA.random_instance()
742 sage: n = J.dimension()
743 sage: x = J.random_element()
744 sage: y = J.random_element()
745 sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
748 B
= self
._inner
_product
_matrix
749 return (B
*x
.to_vector()).inner_product(y
.to_vector())
752 def is_trivial(self
):
754 Return whether or not this algebra is trivial.
756 A trivial algebra contains only the zero element.
760 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
765 sage: J = ComplexHermitianEJA(3)
771 sage: J = TrivialEJA()
776 return self
.dimension() == 0
779 def multiplication_table(self
):
781 Return a visual representation of this algebra's multiplication
782 table (on basis elements).
786 sage: from mjo.eja.eja_algebra import JordanSpinEJA
790 sage: J = JordanSpinEJA(4)
791 sage: J.multiplication_table()
792 +----++----+----+----+----+
793 | * || e0 | e1 | e2 | e3 |
794 +====++====+====+====+====+
795 | e0 || e0 | e1 | e2 | e3 |
796 +----++----+----+----+----+
797 | e1 || e1 | e0 | 0 | 0 |
798 +----++----+----+----+----+
799 | e2 || e2 | 0 | e0 | 0 |
800 +----++----+----+----+----+
801 | e3 || e3 | 0 | 0 | e0 |
802 +----++----+----+----+----+
806 # Prepend the header row.
807 M
= [["*"] + list(self
.gens())]
809 # And to each subsequent row, prepend an entry that belongs to
810 # the left-side "header column."
811 M
+= [ [self
.gens()[i
]] + [ self
.product_on_basis(i
,j
)
815 return table(M
, header_row
=True, header_column
=True, frame
=True)
818 def matrix_basis(self
):
820 Return an (often more natural) representation of this algebras
821 basis as an ordered tuple of matrices.
823 Every finite-dimensional Euclidean Jordan Algebra is a, up to
824 Jordan isomorphism, a direct sum of five simple
825 algebras---four of which comprise Hermitian matrices. And the
826 last type of algebra can of course be thought of as `n`-by-`1`
827 column matrices (ambiguusly called column vectors) to avoid
828 special cases. As a result, matrices (and column vectors) are
829 a natural representation format for Euclidean Jordan algebra
832 But, when we construct an algebra from a basis of matrices,
833 those matrix representations are lost in favor of coordinate
834 vectors *with respect to* that basis. We could eventually
835 convert back if we tried hard enough, but having the original
836 representations handy is valuable enough that we simply store
837 them and return them from this method.
839 Why implement this for non-matrix algebras? Avoiding special
840 cases for the :class:`BilinearFormEJA` pays with simplicity in
841 its own right. But mainly, we would like to be able to assume
842 that elements of a :class:`CartesianProductEJA` can be displayed
843 nicely, without having to have special classes for direct sums
844 one of whose components was a matrix algebra.
848 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
849 ....: RealSymmetricEJA)
853 sage: J = RealSymmetricEJA(2)
855 Finite family {0: e0, 1: e1, 2: e2}
856 sage: J.matrix_basis()
858 [1 0] [ 0 0.7071067811865475?] [0 0]
859 [0 0], [0.7071067811865475? 0], [0 1]
864 sage: J = JordanSpinEJA(2)
866 Finite family {0: e0, 1: e1}
867 sage: J.matrix_basis()
873 return self
._matrix
_basis
876 def matrix_space(self
):
878 Return the matrix space in which this algebra's elements live, if
879 we think of them as matrices (including column vectors of the
882 "By default" this will be an `n`-by-`1` column-matrix space,
883 except when the algebra is trivial. There it's `n`-by-`n`
884 (where `n` is zero), to ensure that two elements of the matrix
885 space (empty matrices) can be multiplied. For algebras of
886 matrices, this returns the space in which their
887 real embeddings live.
891 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
893 ....: QuaternionHermitianEJA,
898 By default, the matrix representation is just a column-matrix
899 equivalent to the vector representation::
901 sage: J = JordanSpinEJA(3)
902 sage: J.matrix_space()
903 Full MatrixSpace of 3 by 1 dense matrices over Algebraic
906 The matrix representation in the trivial algebra is
907 zero-by-zero instead of the usual `n`-by-one::
909 sage: J = TrivialEJA()
910 sage: J.matrix_space()
911 Full MatrixSpace of 0 by 0 dense matrices over Algebraic
914 The matrix space for complex/quaternion Hermitian matrix EJA
915 is the space in which their real-embeddings live, not the
916 original complex/quaternion matrix space::
918 sage: J = ComplexHermitianEJA(2,field=QQ,orthonormalize=False)
919 sage: J.matrix_space()
920 Full MatrixSpace of 4 by 4 dense matrices over Rational Field
921 sage: J = QuaternionHermitianEJA(1,field=QQ,orthonormalize=False)
922 sage: J.matrix_space()
923 Full MatrixSpace of 4 by 4 dense matrices over Rational Field
926 if self
.is_trivial():
927 return MatrixSpace(self
.base_ring(), 0)
929 return self
.matrix_basis()[0].parent()
935 Return the unit element of this algebra.
939 sage: from mjo.eja.eja_algebra import (HadamardEJA,
944 We can compute unit element in the Hadamard EJA::
946 sage: J = HadamardEJA(5)
948 e0 + e1 + e2 + e3 + e4
950 The unit element in the Hadamard EJA is inherited in the
951 subalgebras generated by its elements::
953 sage: J = HadamardEJA(5)
955 e0 + e1 + e2 + e3 + e4
956 sage: x = sum(J.gens())
957 sage: A = x.subalgebra_generated_by(orthonormalize=False)
960 sage: A.one().superalgebra_element()
961 e0 + e1 + e2 + e3 + e4
965 The identity element acts like the identity, regardless of
966 whether or not we orthonormalize::
968 sage: set_random_seed()
969 sage: J = random_eja()
970 sage: x = J.random_element()
971 sage: J.one()*x == x and x*J.one() == x
973 sage: A = x.subalgebra_generated_by()
974 sage: y = A.random_element()
975 sage: A.one()*y == y and y*A.one() == y
980 sage: set_random_seed()
981 sage: J = random_eja(field=QQ, orthonormalize=False)
982 sage: x = J.random_element()
983 sage: J.one()*x == x and x*J.one() == x
985 sage: A = x.subalgebra_generated_by(orthonormalize=False)
986 sage: y = A.random_element()
987 sage: A.one()*y == y and y*A.one() == y
990 The matrix of the unit element's operator is the identity,
991 regardless of the base field and whether or not we
994 sage: set_random_seed()
995 sage: J = random_eja()
996 sage: actual = J.one().operator().matrix()
997 sage: expected = matrix.identity(J.base_ring(), J.dimension())
998 sage: actual == expected
1000 sage: x = J.random_element()
1001 sage: A = x.subalgebra_generated_by()
1002 sage: actual = A.one().operator().matrix()
1003 sage: expected = matrix.identity(A.base_ring(), A.dimension())
1004 sage: actual == expected
1009 sage: set_random_seed()
1010 sage: J = random_eja(field=QQ, orthonormalize=False)
1011 sage: actual = J.one().operator().matrix()
1012 sage: expected = matrix.identity(J.base_ring(), J.dimension())
1013 sage: actual == expected
1015 sage: x = J.random_element()
1016 sage: A = x.subalgebra_generated_by(orthonormalize=False)
1017 sage: actual = A.one().operator().matrix()
1018 sage: expected = matrix.identity(A.base_ring(), A.dimension())
1019 sage: actual == expected
1022 Ensure that the cached unit element (often precomputed by
1023 hand) agrees with the computed one::
1025 sage: set_random_seed()
1026 sage: J = random_eja()
1027 sage: cached = J.one()
1028 sage: J.one.clear_cache()
1029 sage: J.one() == cached
1034 sage: set_random_seed()
1035 sage: J = random_eja(field=QQ, orthonormalize=False)
1036 sage: cached = J.one()
1037 sage: J.one.clear_cache()
1038 sage: J.one() == cached
1042 # We can brute-force compute the matrices of the operators
1043 # that correspond to the basis elements of this algebra.
1044 # If some linear combination of those basis elements is the
1045 # algebra identity, then the same linear combination of
1046 # their matrices has to be the identity matrix.
1048 # Of course, matrices aren't vectors in sage, so we have to
1049 # appeal to the "long vectors" isometry.
1050 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
1052 # Now we use basic linear algebra to find the coefficients,
1053 # of the matrices-as-vectors-linear-combination, which should
1054 # work for the original algebra basis too.
1055 A
= matrix(self
.base_ring(), oper_vecs
)
1057 # We used the isometry on the left-hand side already, but we
1058 # still need to do it for the right-hand side. Recall that we
1059 # wanted something that summed to the identity matrix.
1060 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
1062 # Now if there's an identity element in the algebra, this
1063 # should work. We solve on the left to avoid having to
1064 # transpose the matrix "A".
1065 return self
.from_vector(A
.solve_left(b
))
1068 def peirce_decomposition(self
, c
):
1070 The Peirce decomposition of this algebra relative to the
1073 In the future, this can be extended to a complete system of
1074 orthogonal idempotents.
1078 - ``c`` -- an idempotent of this algebra.
1082 A triple (J0, J5, J1) containing two subalgebras and one subspace
1085 - ``J0`` -- the algebra on the eigenspace of ``c.operator()``
1086 corresponding to the eigenvalue zero.
1088 - ``J5`` -- the eigenspace (NOT a subalgebra) of ``c.operator()``
1089 corresponding to the eigenvalue one-half.
1091 - ``J1`` -- the algebra on the eigenspace of ``c.operator()``
1092 corresponding to the eigenvalue one.
1094 These are the only possible eigenspaces for that operator, and this
1095 algebra is a direct sum of them. The spaces ``J0`` and ``J1`` are
1096 orthogonal, and are subalgebras of this algebra with the appropriate
1101 sage: from mjo.eja.eja_algebra import random_eja, RealSymmetricEJA
1105 The canonical example comes from the symmetric matrices, which
1106 decompose into diagonal and off-diagonal parts::
1108 sage: J = RealSymmetricEJA(3)
1109 sage: C = matrix(QQ, [ [1,0,0],
1113 sage: J0,J5,J1 = J.peirce_decomposition(c)
1115 Euclidean Jordan algebra of dimension 1...
1117 Vector space of degree 6 and dimension 2...
1119 Euclidean Jordan algebra of dimension 3...
1120 sage: J0.one().to_matrix()
1124 sage: orig_df = AA.options.display_format
1125 sage: AA.options.display_format = 'radical'
1126 sage: J.from_vector(J5.basis()[0]).to_matrix()
1130 sage: J.from_vector(J5.basis()[1]).to_matrix()
1134 sage: AA.options.display_format = orig_df
1135 sage: J1.one().to_matrix()
1142 Every algebra decomposes trivially with respect to its identity
1145 sage: set_random_seed()
1146 sage: J = random_eja()
1147 sage: J0,J5,J1 = J.peirce_decomposition(J.one())
1148 sage: J0.dimension() == 0 and J5.dimension() == 0
1150 sage: J1.superalgebra() == J and J1.dimension() == J.dimension()
1153 The decomposition is into eigenspaces, and its components are
1154 therefore necessarily orthogonal. Moreover, the identity
1155 elements in the two subalgebras are the projections onto their
1156 respective subspaces of the superalgebra's identity element::
1158 sage: set_random_seed()
1159 sage: J = random_eja()
1160 sage: x = J.random_element()
1161 sage: if not J.is_trivial():
1162 ....: while x.is_nilpotent():
1163 ....: x = J.random_element()
1164 sage: c = x.subalgebra_idempotent()
1165 sage: J0,J5,J1 = J.peirce_decomposition(c)
1167 sage: for (w,y,z) in zip(J0.basis(), J5.basis(), J1.basis()):
1168 ....: w = w.superalgebra_element()
1169 ....: y = J.from_vector(y)
1170 ....: z = z.superalgebra_element()
1171 ....: ipsum += w.inner_product(y).abs()
1172 ....: ipsum += w.inner_product(z).abs()
1173 ....: ipsum += y.inner_product(z).abs()
1176 sage: J1(c) == J1.one()
1178 sage: J0(J.one() - c) == J0.one()
1182 if not c
.is_idempotent():
1183 raise ValueError("element is not idempotent: %s" % c
)
1185 # Default these to what they should be if they turn out to be
1186 # trivial, because eigenspaces_left() won't return eigenvalues
1187 # corresponding to trivial spaces (e.g. it returns only the
1188 # eigenspace corresponding to lambda=1 if you take the
1189 # decomposition relative to the identity element).
1190 trivial
= self
.subalgebra(())
1191 J0
= trivial
# eigenvalue zero
1192 J5
= VectorSpace(self
.base_ring(), 0) # eigenvalue one-half
1193 J1
= trivial
# eigenvalue one
1195 for (eigval
, eigspace
) in c
.operator().matrix().right_eigenspaces():
1196 if eigval
== ~
(self
.base_ring()(2)):
1199 gens
= tuple( self
.from_vector(b
) for b
in eigspace
.basis() )
1200 subalg
= self
.subalgebra(gens
, check_axioms
=False)
1206 raise ValueError("unexpected eigenvalue: %s" % eigval
)
1211 def random_element(self
, thorough
=False):
1213 Return a random element of this algebra.
1215 Our algebra superclass method only returns a linear
1216 combination of at most two basis elements. We instead
1217 want the vector space "random element" method that
1218 returns a more diverse selection.
1222 - ``thorough`` -- (boolean; default False) whether or not we
1223 should generate irrational coefficients for the random
1224 element when our base ring is irrational; this slows the
1225 algebra operations to a crawl, but any truly random method
1229 # For a general base ring... maybe we can trust this to do the
1230 # right thing? Unlikely, but.
1231 V
= self
.vector_space()
1232 v
= V
.random_element()
1234 if self
.base_ring() is AA
:
1235 # The "random element" method of the algebraic reals is
1236 # stupid at the moment, and only returns integers between
1237 # -2 and 2, inclusive:
1239 # https://trac.sagemath.org/ticket/30875
1241 # Instead, we implement our own "random vector" method,
1242 # and then coerce that into the algebra. We use the vector
1243 # space degree here instead of the dimension because a
1244 # subalgebra could (for example) be spanned by only two
1245 # vectors, each with five coordinates. We need to
1246 # generate all five coordinates.
1248 v
*= QQbar
.random_element().real()
1250 v
*= QQ
.random_element()
1252 return self
.from_vector(V
.coordinate_vector(v
))
1254 def random_elements(self
, count
, thorough
=False):
1256 Return ``count`` random elements as a tuple.
1260 - ``thorough`` -- (boolean; default False) whether or not we
1261 should generate irrational coefficients for the random
1262 elements when our base ring is irrational; this slows the
1263 algebra operations to a crawl, but any truly random method
1268 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1272 sage: J = JordanSpinEJA(3)
1273 sage: x,y,z = J.random_elements(3)
1274 sage: all( [ x in J, y in J, z in J ])
1276 sage: len( J.random_elements(10) ) == 10
1280 return tuple( self
.random_element(thorough
)
1281 for idx
in range(count
) )
1285 def _charpoly_coefficients(self
):
1287 The `r` polynomial coefficients of the "characteristic polynomial
1292 sage: from mjo.eja.eja_algebra import random_eja
1296 The theory shows that these are all homogeneous polynomials of
1299 sage: set_random_seed()
1300 sage: J = random_eja()
1301 sage: all(p.is_homogeneous() for p in J._charpoly_coefficients())
1305 n
= self
.dimension()
1306 R
= self
.coordinate_polynomial_ring()
1308 F
= R
.fraction_field()
1311 # From a result in my book, these are the entries of the
1312 # basis representation of L_x.
1313 return sum( vars[k
]*self
.gens()[k
].operator().matrix()[i
,j
]
1316 L_x
= matrix(F
, n
, n
, L_x_i_j
)
1319 if self
.rank
.is_in_cache():
1321 # There's no need to pad the system with redundant
1322 # columns if we *know* they'll be redundant.
1325 # Compute an extra power in case the rank is equal to
1326 # the dimension (otherwise, we would stop at x^(r-1)).
1327 x_powers
= [ (L_x
**k
)*self
.one().to_vector()
1328 for k
in range(n
+1) ]
1329 A
= matrix
.column(F
, x_powers
[:n
])
1330 AE
= A
.extended_echelon_form()
1337 # The theory says that only the first "r" coefficients are
1338 # nonzero, and they actually live in the original polynomial
1339 # ring and not the fraction field. We negate them because in
1340 # the actual characteristic polynomial, they get moved to the
1341 # other side where x^r lives. We don't bother to trim A_rref
1342 # down to a square matrix and solve the resulting system,
1343 # because the upper-left r-by-r portion of A_rref is
1344 # guaranteed to be the identity matrix, so e.g.
1346 # A_rref.solve_right(Y)
1348 # would just be returning Y.
1349 return (-E
*b
)[:r
].change_ring(R
)
1354 Return the rank of this EJA.
1356 This is a cached method because we know the rank a priori for
1357 all of the algebras we can construct. Thus we can avoid the
1358 expensive ``_charpoly_coefficients()`` call unless we truly
1359 need to compute the whole characteristic polynomial.
1363 sage: from mjo.eja.eja_algebra import (HadamardEJA,
1364 ....: JordanSpinEJA,
1365 ....: RealSymmetricEJA,
1366 ....: ComplexHermitianEJA,
1367 ....: QuaternionHermitianEJA,
1372 The rank of the Jordan spin algebra is always two::
1374 sage: JordanSpinEJA(2).rank()
1376 sage: JordanSpinEJA(3).rank()
1378 sage: JordanSpinEJA(4).rank()
1381 The rank of the `n`-by-`n` Hermitian real, complex, or
1382 quaternion matrices is `n`::
1384 sage: RealSymmetricEJA(4).rank()
1386 sage: ComplexHermitianEJA(3).rank()
1388 sage: QuaternionHermitianEJA(2).rank()
1393 Ensure that every EJA that we know how to construct has a
1394 positive integer rank, unless the algebra is trivial in
1395 which case its rank will be zero::
1397 sage: set_random_seed()
1398 sage: J = random_eja()
1402 sage: r > 0 or (r == 0 and J.is_trivial())
1405 Ensure that computing the rank actually works, since the ranks
1406 of all simple algebras are known and will be cached by default::
1408 sage: set_random_seed() # long time
1409 sage: J = random_eja() # long time
1410 sage: cached = J.rank() # long time
1411 sage: J.rank.clear_cache() # long time
1412 sage: J.rank() == cached # long time
1416 return len(self
._charpoly
_coefficients
())
1419 def subalgebra(self
, basis
, **kwargs
):
1421 Create a subalgebra of this algebra from the given basis.
1423 from mjo
.eja
.eja_subalgebra
import FiniteDimensionalEJASubalgebra
1424 return FiniteDimensionalEJASubalgebra(self
, basis
, **kwargs
)
1427 def vector_space(self
):
1429 Return the vector space that underlies this algebra.
1433 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1437 sage: J = RealSymmetricEJA(2)
1438 sage: J.vector_space()
1439 Vector space of dimension 3 over...
1442 return self
.zero().to_vector().parent().ambient_vector_space()
1446 class RationalBasisEJA(FiniteDimensionalEJA
):
1448 New class for algebras whose supplied basis elements have all rational entries.
1452 sage: from mjo.eja.eja_algebra import BilinearFormEJA
1456 The supplied basis is orthonormalized by default::
1458 sage: B = matrix(QQ, [[1, 0, 0], [0, 25, -32], [0, -32, 41]])
1459 sage: J = BilinearFormEJA(B)
1460 sage: J.matrix_basis()
1477 # Abuse the check_field parameter to check that the entries of
1478 # out basis (in ambient coordinates) are in the field QQ.
1479 if not all( all(b_i
in QQ
for b_i
in b
.list()) for b
in basis
):
1480 raise TypeError("basis not rational")
1482 self
._rational
_algebra
= None
1484 # There's no point in constructing the extra algebra if this
1485 # one is already rational.
1487 # Note: the same Jordan and inner-products work here,
1488 # because they are necessarily defined with respect to
1489 # ambient coordinates and not any particular basis.
1490 self
._rational
_algebra
= FiniteDimensionalEJA(
1495 orthonormalize
=False,
1499 super().__init
__(basis
,
1503 check_field
=check_field
,
1507 def _charpoly_coefficients(self
):
1511 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
1512 ....: JordanSpinEJA)
1516 The base ring of the resulting polynomial coefficients is what
1517 it should be, and not the rationals (unless the algebra was
1518 already over the rationals)::
1520 sage: J = JordanSpinEJA(3)
1521 sage: J._charpoly_coefficients()
1522 (X1^2 - X2^2 - X3^2, -2*X1)
1523 sage: a0 = J._charpoly_coefficients()[0]
1525 Algebraic Real Field
1526 sage: a0.base_ring()
1527 Algebraic Real Field
1530 if self
._rational
_algebra
is None:
1531 # There's no need to construct *another* algebra over the
1532 # rationals if this one is already over the
1533 # rationals. Likewise, if we never orthonormalized our
1534 # basis, we might as well just use the given one.
1535 return super()._charpoly
_coefficients
()
1537 # Do the computation over the rationals. The answer will be
1538 # the same, because all we've done is a change of basis.
1539 # Then, change back from QQ to our real base ring
1540 a
= ( a_i
.change_ring(self
.base_ring())
1541 for a_i
in self
._rational
_algebra
._charpoly
_coefficients
() )
1543 if self
._deortho
_matrix
is None:
1544 # This can happen if our base ring was, say, AA and we
1545 # chose not to (or didn't need to) orthonormalize. It's
1546 # still faster to do the computations over QQ even if
1547 # the numbers in the boxes stay the same.
1550 # Otherwise, convert the coordinate variables back to the
1551 # deorthonormalized ones.
1552 R
= self
.coordinate_polynomial_ring()
1553 from sage
.modules
.free_module_element
import vector
1554 X
= vector(R
, R
.gens())
1555 BX
= self
._deortho
_matrix
*X
1557 subs_dict
= { X[i]: BX[i] for i in range(len(X)) }
1558 return tuple( a_i
.subs(subs_dict
) for a_i
in a
)
1560 class ConcreteEJA(RationalBasisEJA
):
1562 A class for the Euclidean Jordan algebras that we know by name.
1564 These are the Jordan algebras whose basis, multiplication table,
1565 rank, and so on are known a priori. More to the point, they are
1566 the Euclidean Jordan algebras for which we are able to conjure up
1567 a "random instance."
1571 sage: from mjo.eja.eja_algebra import ConcreteEJA
1575 Our basis is normalized with respect to the algebra's inner
1576 product, unless we specify otherwise::
1578 sage: set_random_seed()
1579 sage: J = ConcreteEJA.random_instance()
1580 sage: all( b.norm() == 1 for b in J.gens() )
1583 Since our basis is orthonormal with respect to the algebra's inner
1584 product, and since we know that this algebra is an EJA, any
1585 left-multiplication operator's matrix will be symmetric because
1586 natural->EJA basis representation is an isometry and within the
1587 EJA the operator is self-adjoint by the Jordan axiom::
1589 sage: set_random_seed()
1590 sage: J = ConcreteEJA.random_instance()
1591 sage: x = J.random_element()
1592 sage: x.operator().is_self_adjoint()
1597 def _max_random_instance_size():
1599 Return an integer "size" that is an upper bound on the size of
1600 this algebra when it is used in a random test
1601 case. Unfortunately, the term "size" is ambiguous -- when
1602 dealing with `R^n` under either the Hadamard or Jordan spin
1603 product, the "size" refers to the dimension `n`. When dealing
1604 with a matrix algebra (real symmetric or complex/quaternion
1605 Hermitian), it refers to the size of the matrix, which is far
1606 less than the dimension of the underlying vector space.
1608 This method must be implemented in each subclass.
1610 raise NotImplementedError
1613 def random_instance(cls
, *args
, **kwargs
):
1615 Return a random instance of this type of algebra.
1617 This method should be implemented in each subclass.
1619 from sage
.misc
.prandom
import choice
1620 eja_class
= choice(cls
.__subclasses
__())
1622 # These all bubble up to the RationalBasisEJA superclass
1623 # constructor, so any (kw)args valid there are also valid
1625 return eja_class
.random_instance(*args
, **kwargs
)
1630 def dimension_over_reals():
1632 The dimension of this matrix's base ring over the reals.
1634 The reals are dimension one over themselves, obviously; that's
1635 just `\mathbb{R}^{1}`. Likewise, the complex numbers `a + bi`
1636 have dimension two. Finally, the quaternions have dimension
1637 four over the reals.
1639 This is used to determine the size of the matrix returned from
1640 :meth:`real_embed`, among other things.
1642 raise NotImplementedError
1645 def real_embed(cls
,M
):
1647 Embed the matrix ``M`` into a space of real matrices.
1649 The matrix ``M`` can have entries in any field at the moment:
1650 the real numbers, complex numbers, or quaternions. And although
1651 they are not a field, we can probably support octonions at some
1652 point, too. This function returns a real matrix that "acts like"
1653 the original with respect to matrix multiplication; i.e.
1655 real_embed(M*N) = real_embed(M)*real_embed(N)
1658 if M
.ncols() != M
.nrows():
1659 raise ValueError("the matrix 'M' must be square")
1664 def real_unembed(cls
,M
):
1666 The inverse of :meth:`real_embed`.
1668 if M
.ncols() != M
.nrows():
1669 raise ValueError("the matrix 'M' must be square")
1670 if not ZZ(M
.nrows()).mod(cls
.dimension_over_reals()).is_zero():
1671 raise ValueError("the matrix 'M' must be a real embedding")
1675 def jordan_product(X
,Y
):
1676 return (X
*Y
+ Y
*X
)/2
1679 def trace_inner_product(cls
,X
,Y
):
1681 Compute the trace inner-product of two real-embeddings.
1685 sage: from mjo.eja.eja_algebra import (RealSymmetricEJA,
1686 ....: ComplexHermitianEJA,
1687 ....: QuaternionHermitianEJA)
1691 This gives the same answer as it would if we computed the trace
1692 from the unembedded (original) matrices::
1694 sage: set_random_seed()
1695 sage: J = RealSymmetricEJA.random_instance()
1696 sage: x,y = J.random_elements(2)
1697 sage: Xe = x.to_matrix()
1698 sage: Ye = y.to_matrix()
1699 sage: X = J.real_unembed(Xe)
1700 sage: Y = J.real_unembed(Ye)
1701 sage: expected = (X*Y).trace()
1702 sage: actual = J.trace_inner_product(Xe,Ye)
1703 sage: actual == expected
1708 sage: set_random_seed()
1709 sage: J = ComplexHermitianEJA.random_instance()
1710 sage: x,y = J.random_elements(2)
1711 sage: Xe = x.to_matrix()
1712 sage: Ye = y.to_matrix()
1713 sage: X = J.real_unembed(Xe)
1714 sage: Y = J.real_unembed(Ye)
1715 sage: expected = (X*Y).trace().real()
1716 sage: actual = J.trace_inner_product(Xe,Ye)
1717 sage: actual == expected
1722 sage: set_random_seed()
1723 sage: J = QuaternionHermitianEJA.random_instance()
1724 sage: x,y = J.random_elements(2)
1725 sage: Xe = x.to_matrix()
1726 sage: Ye = y.to_matrix()
1727 sage: X = J.real_unembed(Xe)
1728 sage: Y = J.real_unembed(Ye)
1729 sage: expected = (X*Y).trace().coefficient_tuple()[0]
1730 sage: actual = J.trace_inner_product(Xe,Ye)
1731 sage: actual == expected
1735 Xu
= cls
.real_unembed(X
)
1736 Yu
= cls
.real_unembed(Y
)
1737 tr
= (Xu
*Yu
).trace()
1740 # Works in QQ, AA, RDF, et cetera.
1742 except AttributeError:
1743 # A quaternion doesn't have a real() method, but does
1744 # have coefficient_tuple() method that returns the
1745 # coefficients of 1, i, j, and k -- in that order.
1746 return tr
.coefficient_tuple()[0]
1749 class RealMatrixEJA(MatrixEJA
):
1751 def dimension_over_reals():
1755 class RealSymmetricEJA(ConcreteEJA
, RealMatrixEJA
):
1757 The rank-n simple EJA consisting of real symmetric n-by-n
1758 matrices, the usual symmetric Jordan product, and the trace inner
1759 product. It has dimension `(n^2 + n)/2` over the reals.
1763 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1767 sage: J = RealSymmetricEJA(2)
1768 sage: e0, e1, e2 = J.gens()
1776 In theory, our "field" can be any subfield of the reals::
1778 sage: RealSymmetricEJA(2, field=RDF)
1779 Euclidean Jordan algebra of dimension 3 over Real Double Field
1780 sage: RealSymmetricEJA(2, field=RR)
1781 Euclidean Jordan algebra of dimension 3 over Real Field with
1782 53 bits of precision
1786 The dimension of this algebra is `(n^2 + n) / 2`::
1788 sage: set_random_seed()
1789 sage: n_max = RealSymmetricEJA._max_random_instance_size()
1790 sage: n = ZZ.random_element(1, n_max)
1791 sage: J = RealSymmetricEJA(n)
1792 sage: J.dimension() == (n^2 + n)/2
1795 The Jordan multiplication is what we think it is::
1797 sage: set_random_seed()
1798 sage: J = RealSymmetricEJA.random_instance()
1799 sage: x,y = J.random_elements(2)
1800 sage: actual = (x*y).to_matrix()
1801 sage: X = x.to_matrix()
1802 sage: Y = y.to_matrix()
1803 sage: expected = (X*Y + Y*X)/2
1804 sage: actual == expected
1806 sage: J(expected) == x*y
1809 We can change the generator prefix::
1811 sage: RealSymmetricEJA(3, prefix='q').gens()
1812 (q0, q1, q2, q3, q4, q5)
1814 We can construct the (trivial) algebra of rank zero::
1816 sage: RealSymmetricEJA(0)
1817 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1821 def _denormalized_basis(cls
, n
):
1823 Return a basis for the space of real symmetric n-by-n matrices.
1827 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1831 sage: set_random_seed()
1832 sage: n = ZZ.random_element(1,5)
1833 sage: B = RealSymmetricEJA._denormalized_basis(n)
1834 sage: all( M.is_symmetric() for M in B)
1838 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1842 for j
in range(i
+1):
1843 Eij
= matrix(ZZ
, n
, lambda k
,l
: k
==i
and l
==j
)
1847 Sij
= Eij
+ Eij
.transpose()
1853 def _max_random_instance_size():
1854 return 4 # Dimension 10
1857 def random_instance(cls
, **kwargs
):
1859 Return a random instance of this type of algebra.
1861 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1862 return cls(n
, **kwargs
)
1864 def __init__(self
, n
, **kwargs
):
1865 # We know this is a valid EJA, but will double-check
1866 # if the user passes check_axioms=True.
1867 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
1869 super(RealSymmetricEJA
, self
).__init
__(self
._denormalized
_basis
(n
),
1870 self
.jordan_product
,
1871 self
.trace_inner_product
,
1874 # TODO: this could be factored out somehow, but is left here
1875 # because the MatrixEJA is not presently a subclass of the
1876 # FDEJA class that defines rank() and one().
1877 self
.rank
.set_cache(n
)
1878 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
1879 self
.one
.set_cache(self(idV
))
1883 class ComplexMatrixEJA(MatrixEJA
):
1884 # A manual dictionary-cache for the complex_extension() method,
1885 # since apparently @classmethods can't also be @cached_methods.
1886 _complex_extension
= {}
1889 def complex_extension(cls
,field
):
1891 The complex field that we embed/unembed, as an extension
1892 of the given ``field``.
1894 if field
in cls
._complex
_extension
:
1895 return cls
._complex
_extension
[field
]
1897 # Sage doesn't know how to adjoin the complex "i" (the root of
1898 # x^2 + 1) to a field in a general way. Here, we just enumerate
1899 # all of the cases that I have cared to support so far.
1901 # Sage doesn't know how to embed AA into QQbar, i.e. how
1902 # to adjoin sqrt(-1) to AA.
1904 elif not field
.is_exact():
1906 F
= field
.complex_field()
1908 # Works for QQ and... maybe some other fields.
1909 R
= PolynomialRing(field
, 'z')
1911 F
= field
.extension(z
**2 + 1, 'I', embedding
=CLF(-1).sqrt())
1913 cls
._complex
_extension
[field
] = F
1917 def dimension_over_reals():
1921 def real_embed(cls
,M
):
1923 Embed the n-by-n complex matrix ``M`` into the space of real
1924 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1925 bi` to the block matrix ``[[a,b],[-b,a]]``.
1929 sage: from mjo.eja.eja_algebra import ComplexMatrixEJA
1933 sage: F = QuadraticField(-1, 'I')
1934 sage: x1 = F(4 - 2*i)
1935 sage: x2 = F(1 + 2*i)
1938 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1939 sage: ComplexMatrixEJA.real_embed(M)
1948 Embedding is a homomorphism (isomorphism, in fact)::
1950 sage: set_random_seed()
1951 sage: n = ZZ.random_element(3)
1952 sage: F = QuadraticField(-1, 'I')
1953 sage: X = random_matrix(F, n)
1954 sage: Y = random_matrix(F, n)
1955 sage: Xe = ComplexMatrixEJA.real_embed(X)
1956 sage: Ye = ComplexMatrixEJA.real_embed(Y)
1957 sage: XYe = ComplexMatrixEJA.real_embed(X*Y)
1962 super(ComplexMatrixEJA
,cls
).real_embed(M
)
1965 # We don't need any adjoined elements...
1966 field
= M
.base_ring().base_ring()
1972 blocks
.append(matrix(field
, 2, [ [ a
, b
],
1975 return matrix
.block(field
, n
, blocks
)
1979 def real_unembed(cls
,M
):
1981 The inverse of _embed_complex_matrix().
1985 sage: from mjo.eja.eja_algebra import ComplexMatrixEJA
1989 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1990 ....: [-2, 1, -4, 3],
1991 ....: [ 9, 10, 11, 12],
1992 ....: [-10, 9, -12, 11] ])
1993 sage: ComplexMatrixEJA.real_unembed(A)
1995 [ 10*I + 9 12*I + 11]
1999 Unembedding is the inverse of embedding::
2001 sage: set_random_seed()
2002 sage: F = QuadraticField(-1, 'I')
2003 sage: M = random_matrix(F, 3)
2004 sage: Me = ComplexMatrixEJA.real_embed(M)
2005 sage: ComplexMatrixEJA.real_unembed(Me) == M
2009 super(ComplexMatrixEJA
,cls
).real_unembed(M
)
2011 d
= cls
.dimension_over_reals()
2012 F
= cls
.complex_extension(M
.base_ring())
2015 # Go top-left to bottom-right (reading order), converting every
2016 # 2-by-2 block we see to a single complex element.
2018 for k
in range(n
/d
):
2019 for j
in range(n
/d
):
2020 submat
= M
[d
*k
:d
*k
+d
,d
*j
:d
*j
+d
]
2021 if submat
[0,0] != submat
[1,1]:
2022 raise ValueError('bad on-diagonal submatrix')
2023 if submat
[0,1] != -submat
[1,0]:
2024 raise ValueError('bad off-diagonal submatrix')
2025 z
= submat
[0,0] + submat
[0,1]*i
2028 return matrix(F
, n
/d
, elements
)
2031 class ComplexHermitianEJA(ConcreteEJA
, ComplexMatrixEJA
):
2033 The rank-n simple EJA consisting of complex Hermitian n-by-n
2034 matrices over the real numbers, the usual symmetric Jordan product,
2035 and the real-part-of-trace inner product. It has dimension `n^2` over
2040 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
2044 In theory, our "field" can be any subfield of the reals::
2046 sage: ComplexHermitianEJA(2, field=RDF)
2047 Euclidean Jordan algebra of dimension 4 over Real Double Field
2048 sage: ComplexHermitianEJA(2, field=RR)
2049 Euclidean Jordan algebra of dimension 4 over Real Field with
2050 53 bits of precision
2054 The dimension of this algebra is `n^2`::
2056 sage: set_random_seed()
2057 sage: n_max = ComplexHermitianEJA._max_random_instance_size()
2058 sage: n = ZZ.random_element(1, n_max)
2059 sage: J = ComplexHermitianEJA(n)
2060 sage: J.dimension() == n^2
2063 The Jordan multiplication is what we think it is::
2065 sage: set_random_seed()
2066 sage: J = ComplexHermitianEJA.random_instance()
2067 sage: x,y = J.random_elements(2)
2068 sage: actual = (x*y).to_matrix()
2069 sage: X = x.to_matrix()
2070 sage: Y = y.to_matrix()
2071 sage: expected = (X*Y + Y*X)/2
2072 sage: actual == expected
2074 sage: J(expected) == x*y
2077 We can change the generator prefix::
2079 sage: ComplexHermitianEJA(2, prefix='z').gens()
2082 We can construct the (trivial) algebra of rank zero::
2084 sage: ComplexHermitianEJA(0)
2085 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2090 def _denormalized_basis(cls
, n
):
2092 Returns a basis for the space of complex Hermitian n-by-n matrices.
2094 Why do we embed these? Basically, because all of numerical linear
2095 algebra assumes that you're working with vectors consisting of `n`
2096 entries from a field and scalars from the same field. There's no way
2097 to tell SageMath that (for example) the vectors contain complex
2098 numbers, while the scalar field is real.
2102 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
2106 sage: set_random_seed()
2107 sage: n = ZZ.random_element(1,5)
2108 sage: B = ComplexHermitianEJA._denormalized_basis(n)
2109 sage: all( M.is_symmetric() for M in B)
2114 R
= PolynomialRing(field
, 'z')
2116 F
= field
.extension(z
**2 + 1, 'I')
2119 # This is like the symmetric case, but we need to be careful:
2121 # * We want conjugate-symmetry, not just symmetry.
2122 # * The diagonal will (as a result) be real.
2125 Eij
= matrix
.zero(F
,n
)
2127 for j
in range(i
+1):
2131 Sij
= cls
.real_embed(Eij
)
2134 # The second one has a minus because it's conjugated.
2135 Eij
[j
,i
] = 1 # Eij = Eij + Eij.transpose()
2136 Sij_real
= cls
.real_embed(Eij
)
2138 # Eij = I*Eij - I*Eij.transpose()
2141 Sij_imag
= cls
.real_embed(Eij
)
2147 # Since we embedded these, we can drop back to the "field" that we
2148 # started with instead of the complex extension "F".
2149 return tuple( s
.change_ring(field
) for s
in S
)
2152 def __init__(self
, n
, **kwargs
):
2153 # We know this is a valid EJA, but will double-check
2154 # if the user passes check_axioms=True.
2155 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2157 super(ComplexHermitianEJA
, self
).__init
__(self
._denormalized
_basis
(n
),
2158 self
.jordan_product
,
2159 self
.trace_inner_product
,
2161 # TODO: this could be factored out somehow, but is left here
2162 # because the MatrixEJA is not presently a subclass of the
2163 # FDEJA class that defines rank() and one().
2164 self
.rank
.set_cache(n
)
2165 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
2166 self
.one
.set_cache(self(idV
))
2169 def _max_random_instance_size():
2170 return 3 # Dimension 9
2173 def random_instance(cls
, **kwargs
):
2175 Return a random instance of this type of algebra.
2177 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2178 return cls(n
, **kwargs
)
2180 class QuaternionMatrixEJA(MatrixEJA
):
2182 # A manual dictionary-cache for the quaternion_extension() method,
2183 # since apparently @classmethods can't also be @cached_methods.
2184 _quaternion_extension
= {}
2187 def quaternion_extension(cls
,field
):
2189 The quaternion field that we embed/unembed, as an extension
2190 of the given ``field``.
2192 if field
in cls
._quaternion
_extension
:
2193 return cls
._quaternion
_extension
[field
]
2195 Q
= QuaternionAlgebra(field
,-1,-1)
2197 cls
._quaternion
_extension
[field
] = Q
2201 def dimension_over_reals():
2205 def real_embed(cls
,M
):
2207 Embed the n-by-n quaternion matrix ``M`` into the space of real
2208 matrices of size 4n-by-4n by first sending each quaternion entry `z
2209 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
2210 c+di],[-c + di, a-bi]]`, and then embedding those into a real
2215 sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA
2219 sage: Q = QuaternionAlgebra(QQ,-1,-1)
2220 sage: i,j,k = Q.gens()
2221 sage: x = 1 + 2*i + 3*j + 4*k
2222 sage: M = matrix(Q, 1, [[x]])
2223 sage: QuaternionMatrixEJA.real_embed(M)
2229 Embedding is a homomorphism (isomorphism, in fact)::
2231 sage: set_random_seed()
2232 sage: n = ZZ.random_element(2)
2233 sage: Q = QuaternionAlgebra(QQ,-1,-1)
2234 sage: X = random_matrix(Q, n)
2235 sage: Y = random_matrix(Q, n)
2236 sage: Xe = QuaternionMatrixEJA.real_embed(X)
2237 sage: Ye = QuaternionMatrixEJA.real_embed(Y)
2238 sage: XYe = QuaternionMatrixEJA.real_embed(X*Y)
2243 super(QuaternionMatrixEJA
,cls
).real_embed(M
)
2244 quaternions
= M
.base_ring()
2247 F
= QuadraticField(-1, 'I')
2252 t
= z
.coefficient_tuple()
2257 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
2258 [-c
+ d
*i
, a
- b
*i
]])
2259 realM
= ComplexMatrixEJA
.real_embed(cplxM
)
2260 blocks
.append(realM
)
2262 # We should have real entries by now, so use the realest field
2263 # we've got for the return value.
2264 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
2269 def real_unembed(cls
,M
):
2271 The inverse of _embed_quaternion_matrix().
2275 sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA
2279 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
2280 ....: [-2, 1, -4, 3],
2281 ....: [-3, 4, 1, -2],
2282 ....: [-4, -3, 2, 1]])
2283 sage: QuaternionMatrixEJA.real_unembed(M)
2284 [1 + 2*i + 3*j + 4*k]
2288 Unembedding is the inverse of embedding::
2290 sage: set_random_seed()
2291 sage: Q = QuaternionAlgebra(QQ, -1, -1)
2292 sage: M = random_matrix(Q, 3)
2293 sage: Me = QuaternionMatrixEJA.real_embed(M)
2294 sage: QuaternionMatrixEJA.real_unembed(Me) == M
2298 super(QuaternionMatrixEJA
,cls
).real_unembed(M
)
2300 d
= cls
.dimension_over_reals()
2302 # Use the base ring of the matrix to ensure that its entries can be
2303 # multiplied by elements of the quaternion algebra.
2304 Q
= cls
.quaternion_extension(M
.base_ring())
2307 # Go top-left to bottom-right (reading order), converting every
2308 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
2311 for l
in range(n
/d
):
2312 for m
in range(n
/d
):
2313 submat
= ComplexMatrixEJA
.real_unembed(
2314 M
[d
*l
:d
*l
+d
,d
*m
:d
*m
+d
] )
2315 if submat
[0,0] != submat
[1,1].conjugate():
2316 raise ValueError('bad on-diagonal submatrix')
2317 if submat
[0,1] != -submat
[1,0].conjugate():
2318 raise ValueError('bad off-diagonal submatrix')
2319 z
= submat
[0,0].real()
2320 z
+= submat
[0,0].imag()*i
2321 z
+= submat
[0,1].real()*j
2322 z
+= submat
[0,1].imag()*k
2325 return matrix(Q
, n
/d
, elements
)
2328 class QuaternionHermitianEJA(ConcreteEJA
, QuaternionMatrixEJA
):
2330 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
2331 matrices, the usual symmetric Jordan product, and the
2332 real-part-of-trace inner product. It has dimension `2n^2 - n` over
2337 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2341 In theory, our "field" can be any subfield of the reals::
2343 sage: QuaternionHermitianEJA(2, field=RDF)
2344 Euclidean Jordan algebra of dimension 6 over Real Double Field
2345 sage: QuaternionHermitianEJA(2, field=RR)
2346 Euclidean Jordan algebra of dimension 6 over Real Field with
2347 53 bits of precision
2351 The dimension of this algebra is `2*n^2 - n`::
2353 sage: set_random_seed()
2354 sage: n_max = QuaternionHermitianEJA._max_random_instance_size()
2355 sage: n = ZZ.random_element(1, n_max)
2356 sage: J = QuaternionHermitianEJA(n)
2357 sage: J.dimension() == 2*(n^2) - n
2360 The Jordan multiplication is what we think it is::
2362 sage: set_random_seed()
2363 sage: J = QuaternionHermitianEJA.random_instance()
2364 sage: x,y = J.random_elements(2)
2365 sage: actual = (x*y).to_matrix()
2366 sage: X = x.to_matrix()
2367 sage: Y = y.to_matrix()
2368 sage: expected = (X*Y + Y*X)/2
2369 sage: actual == expected
2371 sage: J(expected) == x*y
2374 We can change the generator prefix::
2376 sage: QuaternionHermitianEJA(2, prefix='a').gens()
2377 (a0, a1, a2, a3, a4, a5)
2379 We can construct the (trivial) algebra of rank zero::
2381 sage: QuaternionHermitianEJA(0)
2382 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2386 def _denormalized_basis(cls
, n
):
2388 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
2390 Why do we embed these? Basically, because all of numerical
2391 linear algebra assumes that you're working with vectors consisting
2392 of `n` entries from a field and scalars from the same field. There's
2393 no way to tell SageMath that (for example) the vectors contain
2394 complex numbers, while the scalar field is real.
2398 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2402 sage: set_random_seed()
2403 sage: n = ZZ.random_element(1,5)
2404 sage: B = QuaternionHermitianEJA._denormalized_basis(n)
2405 sage: all( M.is_symmetric() for M in B )
2410 Q
= QuaternionAlgebra(QQ
,-1,-1)
2413 # This is like the symmetric case, but we need to be careful:
2415 # * We want conjugate-symmetry, not just symmetry.
2416 # * The diagonal will (as a result) be real.
2419 Eij
= matrix
.zero(Q
,n
)
2421 for j
in range(i
+1):
2425 Sij
= cls
.real_embed(Eij
)
2428 # The second, third, and fourth ones have a minus
2429 # because they're conjugated.
2430 # Eij = Eij + Eij.transpose()
2432 Sij_real
= cls
.real_embed(Eij
)
2434 # Eij = I*(Eij - Eij.transpose())
2437 Sij_I
= cls
.real_embed(Eij
)
2439 # Eij = J*(Eij - Eij.transpose())
2442 Sij_J
= cls
.real_embed(Eij
)
2444 # Eij = K*(Eij - Eij.transpose())
2447 Sij_K
= cls
.real_embed(Eij
)
2453 # Since we embedded these, we can drop back to the "field" that we
2454 # started with instead of the quaternion algebra "Q".
2455 return tuple( s
.change_ring(field
) for s
in S
)
2458 def __init__(self
, n
, **kwargs
):
2459 # We know this is a valid EJA, but will double-check
2460 # if the user passes check_axioms=True.
2461 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2463 super(QuaternionHermitianEJA
, self
).__init
__(self
._denormalized
_basis
(n
),
2464 self
.jordan_product
,
2465 self
.trace_inner_product
,
2467 # TODO: this could be factored out somehow, but is left here
2468 # because the MatrixEJA is not presently a subclass of the
2469 # FDEJA class that defines rank() and one().
2470 self
.rank
.set_cache(n
)
2471 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
2472 self
.one
.set_cache(self(idV
))
2476 def _max_random_instance_size():
2478 The maximum rank of a random QuaternionHermitianEJA.
2480 return 2 # Dimension 6
2483 def random_instance(cls
, **kwargs
):
2485 Return a random instance of this type of algebra.
2487 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2488 return cls(n
, **kwargs
)
2491 class HadamardEJA(ConcreteEJA
):
2493 Return the Euclidean Jordan Algebra corresponding to the set
2494 `R^n` under the Hadamard product.
2496 Note: this is nothing more than the Cartesian product of ``n``
2497 copies of the spin algebra. Once Cartesian product algebras
2498 are implemented, this can go.
2502 sage: from mjo.eja.eja_algebra import HadamardEJA
2506 This multiplication table can be verified by hand::
2508 sage: J = HadamardEJA(3)
2509 sage: e0,e1,e2 = J.gens()
2525 We can change the generator prefix::
2527 sage: HadamardEJA(3, prefix='r').gens()
2531 def __init__(self
, n
, **kwargs
):
2533 jordan_product
= lambda x
,y
: x
2534 inner_product
= lambda x
,y
: x
2536 def jordan_product(x
,y
):
2538 return P( xi
*yi
for (xi
,yi
) in zip(x
,y
) )
2540 def inner_product(x
,y
):
2543 # New defaults for keyword arguments. Don't orthonormalize
2544 # because our basis is already orthonormal with respect to our
2545 # inner-product. Don't check the axioms, because we know this
2546 # is a valid EJA... but do double-check if the user passes
2547 # check_axioms=True. Note: we DON'T override the "check_field"
2548 # default here, because the user can pass in a field!
2549 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2550 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2552 column_basis
= tuple( b
.column() for b
in FreeModule(ZZ
, n
).basis() )
2553 super().__init
__(column_basis
,
2558 self
.rank
.set_cache(n
)
2561 self
.one
.set_cache( self
.zero() )
2563 self
.one
.set_cache( sum(self
.gens()) )
2566 def _max_random_instance_size():
2568 The maximum dimension of a random HadamardEJA.
2573 def random_instance(cls
, **kwargs
):
2575 Return a random instance of this type of algebra.
2577 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2578 return cls(n
, **kwargs
)
2581 class BilinearFormEJA(ConcreteEJA
):
2583 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2584 with the half-trace inner product and jordan product ``x*y =
2585 (<Bx,y>,y_bar>, x0*y_bar + y0*x_bar)`` where `B = 1 \times B22` is
2586 a symmetric positive-definite "bilinear form" matrix. Its
2587 dimension is the size of `B`, and it has rank two in dimensions
2588 larger than two. It reduces to the ``JordanSpinEJA`` when `B` is
2589 the identity matrix of order ``n``.
2591 We insist that the one-by-one upper-left identity block of `B` be
2592 passed in as well so that we can be passed a matrix of size zero
2593 to construct a trivial algebra.
2597 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
2598 ....: JordanSpinEJA)
2602 When no bilinear form is specified, the identity matrix is used,
2603 and the resulting algebra is the Jordan spin algebra::
2605 sage: B = matrix.identity(AA,3)
2606 sage: J0 = BilinearFormEJA(B)
2607 sage: J1 = JordanSpinEJA(3)
2608 sage: J0.multiplication_table() == J0.multiplication_table()
2611 An error is raised if the matrix `B` does not correspond to a
2612 positive-definite bilinear form::
2614 sage: B = matrix.random(QQ,2,3)
2615 sage: J = BilinearFormEJA(B)
2616 Traceback (most recent call last):
2618 ValueError: bilinear form is not positive-definite
2619 sage: B = matrix.zero(QQ,3)
2620 sage: J = BilinearFormEJA(B)
2621 Traceback (most recent call last):
2623 ValueError: bilinear form is not positive-definite
2627 We can create a zero-dimensional algebra::
2629 sage: B = matrix.identity(AA,0)
2630 sage: J = BilinearFormEJA(B)
2634 We can check the multiplication condition given in the Jordan, von
2635 Neumann, and Wigner paper (and also discussed on my "On the
2636 symmetry..." paper). Note that this relies heavily on the standard
2637 choice of basis, as does anything utilizing the bilinear form
2638 matrix. We opt not to orthonormalize the basis, because if we
2639 did, we would have to normalize the `s_{i}` in a similar manner::
2641 sage: set_random_seed()
2642 sage: n = ZZ.random_element(5)
2643 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2644 sage: B11 = matrix.identity(QQ,1)
2645 sage: B22 = M.transpose()*M
2646 sage: B = block_matrix(2,2,[ [B11,0 ],
2648 sage: J = BilinearFormEJA(B, orthonormalize=False)
2649 sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
2650 sage: V = J.vector_space()
2651 sage: sis = [ J( V([0] + (M.inverse()*ei).list()).column() )
2652 ....: for ei in eis ]
2653 sage: actual = [ sis[i]*sis[j]
2654 ....: for i in range(n-1)
2655 ....: for j in range(n-1) ]
2656 sage: expected = [ J.one() if i == j else J.zero()
2657 ....: for i in range(n-1)
2658 ....: for j in range(n-1) ]
2659 sage: actual == expected
2663 def __init__(self
, B
, **kwargs
):
2664 # The matrix "B" is supplied by the user in most cases,
2665 # so it makes sense to check whether or not its positive-
2666 # definite unless we are specifically asked not to...
2667 if ("check_axioms" not in kwargs
) or kwargs
["check_axioms"]:
2668 if not B
.is_positive_definite():
2669 raise ValueError("bilinear form is not positive-definite")
2671 # However, all of the other data for this EJA is computed
2672 # by us in manner that guarantees the axioms are
2673 # satisfied. So, again, unless we are specifically asked to
2674 # verify things, we'll skip the rest of the checks.
2675 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2677 def inner_product(x
,y
):
2678 return (y
.T
*B
*x
)[0,0]
2680 def jordan_product(x
,y
):
2686 z0
= inner_product(y
,x
)
2687 zbar
= y0
*xbar
+ x0
*ybar
2688 return P([z0
] + zbar
.list())
2691 column_basis
= tuple( b
.column() for b
in FreeModule(ZZ
, n
).basis() )
2692 super(BilinearFormEJA
, self
).__init
__(column_basis
,
2697 # The rank of this algebra is two, unless we're in a
2698 # one-dimensional ambient space (because the rank is bounded
2699 # by the ambient dimension).
2700 self
.rank
.set_cache(min(n
,2))
2703 self
.one
.set_cache( self
.zero() )
2705 self
.one
.set_cache( self
.monomial(0) )
2708 def _max_random_instance_size():
2710 The maximum dimension of a random BilinearFormEJA.
2715 def random_instance(cls
, **kwargs
):
2717 Return a random instance of this algebra.
2719 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2721 B
= matrix
.identity(ZZ
, n
)
2722 return cls(B
, **kwargs
)
2724 B11
= matrix
.identity(ZZ
, 1)
2725 M
= matrix
.random(ZZ
, n
-1)
2726 I
= matrix
.identity(ZZ
, n
-1)
2728 while alpha
.is_zero():
2729 alpha
= ZZ
.random_element().abs()
2730 B22
= M
.transpose()*M
+ alpha
*I
2732 from sage
.matrix
.special
import block_matrix
2733 B
= block_matrix(2,2, [ [B11
, ZZ(0) ],
2736 return cls(B
, **kwargs
)
2739 class JordanSpinEJA(BilinearFormEJA
):
2741 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2742 with the usual inner product and jordan product ``x*y =
2743 (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2748 sage: from mjo.eja.eja_algebra import JordanSpinEJA
2752 This multiplication table can be verified by hand::
2754 sage: J = JordanSpinEJA(4)
2755 sage: e0,e1,e2,e3 = J.gens()
2771 We can change the generator prefix::
2773 sage: JordanSpinEJA(2, prefix='B').gens()
2778 Ensure that we have the usual inner product on `R^n`::
2780 sage: set_random_seed()
2781 sage: J = JordanSpinEJA.random_instance()
2782 sage: x,y = J.random_elements(2)
2783 sage: actual = x.inner_product(y)
2784 sage: expected = x.to_vector().inner_product(y.to_vector())
2785 sage: actual == expected
2789 def __init__(self
, n
, **kwargs
):
2790 # This is a special case of the BilinearFormEJA with the
2791 # identity matrix as its bilinear form.
2792 B
= matrix
.identity(ZZ
, n
)
2794 # Don't orthonormalize because our basis is already
2795 # orthonormal with respect to our inner-product.
2796 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2798 # But also don't pass check_field=False here, because the user
2799 # can pass in a field!
2800 super(JordanSpinEJA
, self
).__init
__(B
, **kwargs
)
2803 def _max_random_instance_size():
2805 The maximum dimension of a random JordanSpinEJA.
2810 def random_instance(cls
, **kwargs
):
2812 Return a random instance of this type of algebra.
2814 Needed here to override the implementation for ``BilinearFormEJA``.
2816 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2817 return cls(n
, **kwargs
)
2820 class TrivialEJA(ConcreteEJA
):
2822 The trivial Euclidean Jordan algebra consisting of only a zero element.
2826 sage: from mjo.eja.eja_algebra import TrivialEJA
2830 sage: J = TrivialEJA()
2837 sage: 7*J.one()*12*J.one()
2839 sage: J.one().inner_product(J.one())
2841 sage: J.one().norm()
2843 sage: J.one().subalgebra_generated_by()
2844 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2849 def __init__(self
, **kwargs
):
2850 jordan_product
= lambda x
,y
: x
2851 inner_product
= lambda x
,y
: 0
2854 # New defaults for keyword arguments
2855 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2856 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2858 super(TrivialEJA
, self
).__init
__(basis
,
2862 # The rank is zero using my definition, namely the dimension of the
2863 # largest subalgebra generated by any element.
2864 self
.rank
.set_cache(0)
2865 self
.one
.set_cache( self
.zero() )
2868 def random_instance(cls
, **kwargs
):
2869 # We don't take a "size" argument so the superclass method is
2870 # inappropriate for us.
2871 return cls(**kwargs
)
2874 class CartesianProductEJA(CombinatorialFreeModule_CartesianProduct
,
2875 FiniteDimensionalEJA
):
2877 The external (orthogonal) direct sum of two or more Euclidean
2878 Jordan algebras. Every Euclidean Jordan algebra decomposes into an
2879 orthogonal direct sum of simple Euclidean Jordan algebras which is
2880 then isometric to a Cartesian product, so no generality is lost by
2881 providing only this construction.
2885 sage: from mjo.eja.eja_algebra import (random_eja,
2886 ....: CartesianProductEJA,
2888 ....: JordanSpinEJA,
2889 ....: RealSymmetricEJA)
2893 The Jordan product is inherited from our factors and implemented by
2894 our CombinatorialFreeModule Cartesian product superclass::
2896 sage: set_random_seed()
2897 sage: J1 = HadamardEJA(2)
2898 sage: J2 = RealSymmetricEJA(2)
2899 sage: J = cartesian_product([J1,J2])
2900 sage: x,y = J.random_elements(2)
2904 The ability to retrieve the original factors is implemented by our
2905 CombinatorialFreeModule Cartesian product superclass::
2907 sage: J1 = HadamardEJA(2, field=QQ)
2908 sage: J2 = JordanSpinEJA(3, field=QQ)
2909 sage: J = cartesian_product([J1,J2])
2910 sage: J.cartesian_factors()
2911 (Euclidean Jordan algebra of dimension 2 over Rational Field,
2912 Euclidean Jordan algebra of dimension 3 over Rational Field)
2914 You can provide more than two factors::
2916 sage: J1 = HadamardEJA(2)
2917 sage: J2 = JordanSpinEJA(3)
2918 sage: J3 = RealSymmetricEJA(3)
2919 sage: cartesian_product([J1,J2,J3])
2920 Euclidean Jordan algebra of dimension 2 over Algebraic Real
2921 Field (+) Euclidean Jordan algebra of dimension 3 over Algebraic
2922 Real Field (+) Euclidean Jordan algebra of dimension 6 over
2923 Algebraic Real Field
2925 Rank is additive on a Cartesian product::
2927 sage: J1 = HadamardEJA(1)
2928 sage: J2 = RealSymmetricEJA(2)
2929 sage: J = cartesian_product([J1,J2])
2930 sage: J1.rank.clear_cache()
2931 sage: J2.rank.clear_cache()
2932 sage: J.rank.clear_cache()
2935 sage: J.rank() == J1.rank() + J2.rank()
2938 The same rank computation works over the rationals, with whatever
2941 sage: J1 = HadamardEJA(1, field=QQ, orthonormalize=False)
2942 sage: J2 = RealSymmetricEJA(2, field=QQ, orthonormalize=False)
2943 sage: J = cartesian_product([J1,J2])
2944 sage: J1.rank.clear_cache()
2945 sage: J2.rank.clear_cache()
2946 sage: J.rank.clear_cache()
2949 sage: J.rank() == J1.rank() + J2.rank()
2952 The product algebra will be associative if and only if all of its
2953 components are associative::
2955 sage: J1 = HadamardEJA(2)
2956 sage: J1.is_associative()
2958 sage: J2 = HadamardEJA(3)
2959 sage: J2.is_associative()
2961 sage: J3 = RealSymmetricEJA(3)
2962 sage: J3.is_associative()
2964 sage: CP1 = cartesian_product([J1,J2])
2965 sage: CP1.is_associative()
2967 sage: CP2 = cartesian_product([J1,J3])
2968 sage: CP2.is_associative()
2973 All factors must share the same base field::
2975 sage: J1 = HadamardEJA(2, field=QQ)
2976 sage: J2 = RealSymmetricEJA(2)
2977 sage: CartesianProductEJA((J1,J2))
2978 Traceback (most recent call last):
2980 ValueError: all factors must share the same base field
2982 The cached unit element is the same one that would be computed::
2984 sage: set_random_seed() # long time
2985 sage: J1 = random_eja() # long time
2986 sage: J2 = random_eja() # long time
2987 sage: J = cartesian_product([J1,J2]) # long time
2988 sage: actual = J.one() # long time
2989 sage: J.one.clear_cache() # long time
2990 sage: expected = J.one() # long time
2991 sage: actual == expected # long time
2995 Element
= FiniteDimensionalEJAElement
2998 def __init__(self
, algebras
, **kwargs
):
2999 CombinatorialFreeModule_CartesianProduct
.__init
__(self
,
3002 field
= algebras
[0].base_ring()
3003 if not all( J
.base_ring() == field
for J
in algebras
):
3004 raise ValueError("all factors must share the same base field")
3006 associative
= all( m
.is_associative() for m
in algebras
)
3008 # The definition of matrix_space() and self.basis() relies
3009 # only on the stuff in the CFM_CartesianProduct class, which
3010 # we've already initialized.
3011 Js
= self
.cartesian_factors()
3013 MS
= self
.matrix_space()
3015 MS(tuple( self
.cartesian_projection(i
)(b
).to_matrix()
3016 for i
in range(m
) ))
3017 for b
in self
.basis()
3020 # Define jordan/inner products that operate on that matrix_basis.
3021 def jordan_product(x
,y
):
3023 (Js
[i
](x
[i
])*Js
[i
](y
[i
])).to_matrix() for i
in range(m
)
3026 def inner_product(x
, y
):
3028 Js
[i
](x
[i
]).inner_product(Js
[i
](y
[i
])) for i
in range(m
)
3031 # There's no need to check the field since it already came
3032 # from an EJA. Likewise the axioms are guaranteed to be
3033 # satisfied, unless the guy writing this class sucks.
3035 # If you want the basis to be orthonormalized, orthonormalize
3037 FiniteDimensionalEJA
.__init
__(self
,
3042 orthonormalize
=False,
3043 associative
=associative
,
3044 cartesian_product
=True,
3048 ones
= tuple(J
.one() for J
in algebras
)
3049 self
.one
.set_cache(self
._cartesian
_product
_of
_elements
(ones
))
3050 self
.rank
.set_cache(sum(J
.rank() for J
in algebras
))
3052 def matrix_space(self
):
3054 Return the space that our matrix basis lives in as a Cartesian
3059 sage: from mjo.eja.eja_algebra import (HadamardEJA,
3060 ....: RealSymmetricEJA)
3064 sage: J1 = HadamardEJA(1)
3065 sage: J2 = RealSymmetricEJA(2)
3066 sage: J = cartesian_product([J1,J2])
3067 sage: J.matrix_space()
3068 The Cartesian product of (Full MatrixSpace of 1 by 1 dense
3069 matrices over Algebraic Real Field, Full MatrixSpace of 2
3070 by 2 dense matrices over Algebraic Real Field)
3073 from sage
.categories
.cartesian_product
import cartesian_product
3074 return cartesian_product( [J
.matrix_space()
3075 for J
in self
.cartesian_factors()] )
3078 def cartesian_projection(self
, i
):
3082 sage: from mjo.eja.eja_algebra import (random_eja,
3083 ....: JordanSpinEJA,
3085 ....: RealSymmetricEJA,
3086 ....: ComplexHermitianEJA)
3090 The projection morphisms are Euclidean Jordan algebra
3093 sage: J1 = HadamardEJA(2)
3094 sage: J2 = RealSymmetricEJA(2)
3095 sage: J = cartesian_product([J1,J2])
3096 sage: J.cartesian_projection(0)
3097 Linear operator between finite-dimensional Euclidean Jordan
3098 algebras represented by the matrix:
3101 Domain: Euclidean Jordan algebra of dimension 2 over Algebraic
3102 Real Field (+) Euclidean Jordan algebra of dimension 3 over
3103 Algebraic Real Field
3104 Codomain: Euclidean Jordan algebra of dimension 2 over Algebraic
3106 sage: J.cartesian_projection(1)
3107 Linear operator between finite-dimensional Euclidean Jordan
3108 algebras represented by the matrix:
3112 Domain: Euclidean Jordan algebra of dimension 2 over Algebraic
3113 Real Field (+) Euclidean Jordan algebra of dimension 3 over
3114 Algebraic Real Field
3115 Codomain: Euclidean Jordan algebra of dimension 3 over Algebraic
3118 The projections work the way you'd expect on the vector
3119 representation of an element::
3121 sage: J1 = JordanSpinEJA(2)
3122 sage: J2 = ComplexHermitianEJA(2)
3123 sage: J = cartesian_product([J1,J2])
3124 sage: pi_left = J.cartesian_projection(0)
3125 sage: pi_right = J.cartesian_projection(1)
3126 sage: pi_left(J.one()).to_vector()
3128 sage: pi_right(J.one()).to_vector()
3130 sage: J.one().to_vector()
3135 The answer never changes::
3137 sage: set_random_seed()
3138 sage: J1 = random_eja()
3139 sage: J2 = random_eja()
3140 sage: J = cartesian_product([J1,J2])
3141 sage: P0 = J.cartesian_projection(0)
3142 sage: P1 = J.cartesian_projection(0)
3147 Ji
= self
.cartesian_factors()[i
]
3148 # Requires the fix on Trac 31421/31422 to work!
3149 Pi
= super().cartesian_projection(i
)
3150 return FiniteDimensionalEJAOperator(self
,Ji
,Pi
.matrix())
3153 def cartesian_embedding(self
, i
):
3157 sage: from mjo.eja.eja_algebra import (random_eja,
3158 ....: JordanSpinEJA,
3160 ....: RealSymmetricEJA)
3164 The embedding morphisms are Euclidean Jordan algebra
3167 sage: J1 = HadamardEJA(2)
3168 sage: J2 = RealSymmetricEJA(2)
3169 sage: J = cartesian_product([J1,J2])
3170 sage: J.cartesian_embedding(0)
3171 Linear operator between finite-dimensional Euclidean Jordan
3172 algebras represented by the matrix:
3178 Domain: Euclidean Jordan algebra of dimension 2 over
3179 Algebraic Real Field
3180 Codomain: Euclidean Jordan algebra of dimension 2 over
3181 Algebraic Real Field (+) Euclidean Jordan algebra of
3182 dimension 3 over Algebraic Real Field
3183 sage: J.cartesian_embedding(1)
3184 Linear operator between finite-dimensional Euclidean Jordan
3185 algebras represented by the matrix:
3191 Domain: Euclidean Jordan algebra of dimension 3 over
3192 Algebraic Real Field
3193 Codomain: Euclidean Jordan algebra of dimension 2 over
3194 Algebraic Real Field (+) Euclidean Jordan algebra of
3195 dimension 3 over Algebraic Real Field
3197 The embeddings work the way you'd expect on the vector
3198 representation of an element::
3200 sage: J1 = JordanSpinEJA(3)
3201 sage: J2 = RealSymmetricEJA(2)
3202 sage: J = cartesian_product([J1,J2])
3203 sage: iota_left = J.cartesian_embedding(0)
3204 sage: iota_right = J.cartesian_embedding(1)
3205 sage: iota_left(J1.zero()) == J.zero()
3207 sage: iota_right(J2.zero()) == J.zero()
3209 sage: J1.one().to_vector()
3211 sage: iota_left(J1.one()).to_vector()
3213 sage: J2.one().to_vector()
3215 sage: iota_right(J2.one()).to_vector()
3217 sage: J.one().to_vector()
3222 The answer never changes::
3224 sage: set_random_seed()
3225 sage: J1 = random_eja()
3226 sage: J2 = random_eja()
3227 sage: J = cartesian_product([J1,J2])
3228 sage: E0 = J.cartesian_embedding(0)
3229 sage: E1 = J.cartesian_embedding(0)
3233 Composing a projection with the corresponding inclusion should
3234 produce the identity map, and mismatching them should produce
3237 sage: set_random_seed()
3238 sage: J1 = random_eja()
3239 sage: J2 = random_eja()
3240 sage: J = cartesian_product([J1,J2])
3241 sage: iota_left = J.cartesian_embedding(0)
3242 sage: iota_right = J.cartesian_embedding(1)
3243 sage: pi_left = J.cartesian_projection(0)
3244 sage: pi_right = J.cartesian_projection(1)
3245 sage: pi_left*iota_left == J1.one().operator()
3247 sage: pi_right*iota_right == J2.one().operator()
3249 sage: (pi_left*iota_right).is_zero()
3251 sage: (pi_right*iota_left).is_zero()
3255 Ji
= self
.cartesian_factors()[i
]
3256 # Requires the fix on Trac 31421/31422 to work!
3257 Ei
= super().cartesian_embedding(i
)
3258 return FiniteDimensionalEJAOperator(Ji
,self
,Ei
.matrix())
3262 FiniteDimensionalEJA
.CartesianProduct
= CartesianProductEJA
3264 class RationalBasisCartesianProductEJA(CartesianProductEJA
,
3267 A separate class for products of algebras for which we know a
3272 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
3273 ....: RealSymmetricEJA)
3277 This gives us fast characteristic polynomial computations in
3278 product algebras, too::
3281 sage: J1 = JordanSpinEJA(2)
3282 sage: J2 = RealSymmetricEJA(3)
3283 sage: J = cartesian_product([J1,J2])
3284 sage: J.characteristic_polynomial_of().degree()
3290 def __init__(self
, algebras
, **kwargs
):
3291 CartesianProductEJA
.__init
__(self
, algebras
, **kwargs
)
3293 self
._rational
_algebra
= None
3294 if self
.vector_space().base_field() is not QQ
:
3295 self
._rational
_algebra
= cartesian_product([
3296 r
._rational
_algebra
for r
in algebras
3300 RationalBasisEJA
.CartesianProduct
= RationalBasisCartesianProductEJA
3302 random_eja
= ConcreteEJA
.random_instance