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 of basis elements in "matrix form," which
88 must be the same form as the arguments to ``jordan_product``
89 and ``inner_product``. In reality, "matrix form" can be either
90 vectors, matrices, or a Cartesian product (ordered tuple)
91 of vectors or matrices. All of these would ideally be vector
92 spaces in sage with no special-casing needed; but in reality
93 we turn vectors into column-matrices and Cartesian products
94 `(a,b)` into column matrices `(a,b)^{T}` after converting
95 `a` and `b` themselves.
97 - jordan_product -- function of two ``basis`` elements (in
98 matrix form) that returns their jordan product, also in matrix
99 form; this will be applied to ``basis`` to compute a
100 multiplication table for the algebra.
102 - inner_product -- function of two ``basis`` elements (in matrix
103 form) that returns their inner product. This will be applied
104 to ``basis`` to compute an inner-product table (basically a
105 matrix) for this algebra.
107 Element
= FiniteDimensionalEJAElement
116 cartesian_product
=False,
121 # Keep track of whether or not the matrix basis consists of
122 # tuples, since we need special cases for them damned near
123 # everywhere. This is INDEPENDENT of whether or not the
124 # algebra is a cartesian product, since a subalgebra of a
125 # cartesian product will have a basis of tuples, but will not
126 # in general itself be a cartesian product algebra.
127 self
._matrix
_basis
_is
_cartesian
= False
130 if hasattr(basis
[0], 'cartesian_factors'):
131 self
._matrix
_basis
_is
_cartesian
= True
134 if not field
.is_subring(RR
):
135 # Note: this does return true for the real algebraic
136 # field, the rationals, and any quadratic field where
137 # we've specified a real embedding.
138 raise ValueError("scalar field is not real")
140 # If the basis given to us wasn't over the field that it's
141 # supposed to be over, fix that. Or, you know, crash.
142 if not cartesian_product
:
143 # The field for a cartesian product algebra comes from one
144 # of its factors and is the same for all factors, so
145 # there's no need to "reapply" it on product algebras.
146 if self
._matrix
_basis
_is
_cartesian
:
147 # OK since if n == 0, the basis does not consist of tuples.
148 P
= basis
[0].parent()
149 basis
= tuple( P(tuple(b_i
.change_ring(field
) for b_i
in b
))
152 basis
= tuple( b
.change_ring(field
) for b
in basis
)
156 # Check commutativity of the Jordan and inner-products.
157 # This has to be done before we build the multiplication
158 # and inner-product tables/matrices, because we take
159 # advantage of symmetry in the process.
160 if not all( jordan_product(bi
,bj
) == jordan_product(bj
,bi
)
163 raise ValueError("Jordan product is not commutative")
165 if not all( inner_product(bi
,bj
) == inner_product(bj
,bi
)
168 raise ValueError("inner-product is not commutative")
171 category
= MagmaticAlgebras(field
).FiniteDimensional()
172 category
= category
.WithBasis().Unital()
174 # Element subalgebras can take advantage of this.
175 category
= category
.Associative()
176 if cartesian_product
:
177 category
= category
.CartesianProducts()
179 # Call the superclass constructor so that we can use its from_vector()
180 # method to build our multiplication table.
181 CombinatorialFreeModule
.__init
__(self
,
188 # Now comes all of the hard work. We'll be constructing an
189 # ambient vector space V that our (vectorized) basis lives in,
190 # as well as a subspace W of V spanned by those (vectorized)
191 # basis elements. The W-coordinates are the coefficients that
192 # we see in things like x = 1*e1 + 2*e2.
197 degree
= len(_all2list(basis
[0]))
199 # Build an ambient space that fits our matrix basis when
200 # written out as "long vectors."
201 V
= VectorSpace(field
, degree
)
203 # The matrix that will hole the orthonormal -> unorthonormal
204 # coordinate transformation.
205 self
._deortho
_matrix
= None
208 # Save a copy of the un-orthonormalized basis for later.
209 # Convert it to ambient V (vector) coordinates while we're
210 # at it, because we'd have to do it later anyway.
211 deortho_vector_basis
= tuple( V(_all2list(b
)) for b
in basis
)
213 from mjo
.eja
.eja_utils
import gram_schmidt
214 basis
= tuple(gram_schmidt(basis
, inner_product
))
216 # Save the (possibly orthonormalized) matrix basis for
218 self
._matrix
_basis
= basis
220 # Now create the vector space for the algebra, which will have
221 # its own set of non-ambient coordinates (in terms of the
223 vector_basis
= tuple( V(_all2list(b
)) for b
in basis
)
224 W
= V
.span_of_basis( vector_basis
, check
=check_axioms
)
227 # Now "W" is the vector space of our algebra coordinates. The
228 # variables "X1", "X2",... refer to the entries of vectors in
229 # W. Thus to convert back and forth between the orthonormal
230 # coordinates and the given ones, we need to stick the original
232 U
= V
.span_of_basis( deortho_vector_basis
, check
=check_axioms
)
233 self
._deortho
_matrix
= matrix( U
.coordinate_vector(q
)
234 for q
in vector_basis
)
237 # Now we actually compute the multiplication and inner-product
238 # tables/matrices using the possibly-orthonormalized basis.
239 self
._inner
_product
_matrix
= matrix
.identity(field
, n
)
240 self
._multiplication
_table
= [ [0 for j
in range(i
+1)]
243 # Note: the Jordan and inner-products are defined in terms
244 # of the ambient basis. It's important that their arguments
245 # are in ambient coordinates as well.
248 # ortho basis w.r.t. ambient coords
252 # The jordan product returns a matrixy answer, so we
253 # have to convert it to the algebra coordinates.
254 elt
= jordan_product(q_i
, q_j
)
255 elt
= W
.coordinate_vector(V(_all2list(elt
)))
256 self
._multiplication
_table
[i
][j
] = self
.from_vector(elt
)
258 if not orthonormalize
:
259 # If we're orthonormalizing the basis with respect
260 # to an inner-product, then the inner-product
261 # matrix with respect to the resulting basis is
262 # just going to be the identity.
263 ip
= inner_product(q_i
, q_j
)
264 self
._inner
_product
_matrix
[i
,j
] = ip
265 self
._inner
_product
_matrix
[j
,i
] = ip
267 self
._inner
_product
_matrix
._cache
= {'hermitian': True}
268 self
._inner
_product
_matrix
.set_immutable()
271 if not self
._is
_jordanian
():
272 raise ValueError("Jordan identity does not hold")
273 if not self
._inner
_product
_is
_associative
():
274 raise ValueError("inner product is not associative")
277 def _coerce_map_from_base_ring(self
):
279 Disable the map from the base ring into the algebra.
281 Performing a nonsense conversion like this automatically
282 is counterpedagogical. The fallback is to try the usual
283 element constructor, which should also fail.
287 sage: from mjo.eja.eja_algebra import random_eja
291 sage: set_random_seed()
292 sage: J = random_eja()
294 Traceback (most recent call last):
296 ValueError: not an element of this algebra
302 def product_on_basis(self
, i
, j
):
304 Returns the Jordan product of the `i` and `j`th basis elements.
306 This completely defines the Jordan product on the algebra, and
307 is used direclty by our superclass machinery to implement
312 sage: from mjo.eja.eja_algebra import random_eja
316 sage: set_random_seed()
317 sage: J = random_eja()
318 sage: n = J.dimension()
321 sage: ei_ej = J.zero()*J.zero()
323 ....: i = ZZ.random_element(n)
324 ....: j = ZZ.random_element(n)
325 ....: ei = J.gens()[i]
326 ....: ej = J.gens()[j]
327 ....: ei_ej = J.product_on_basis(i,j)
332 # We only stored the lower-triangular portion of the
333 # multiplication table.
335 return self
._multiplication
_table
[i
][j
]
337 return self
._multiplication
_table
[j
][i
]
339 def inner_product(self
, x
, y
):
341 The inner product associated with this Euclidean Jordan algebra.
343 Defaults to the trace inner product, but can be overridden by
344 subclasses if they are sure that the necessary properties are
349 sage: from mjo.eja.eja_algebra import (random_eja,
351 ....: BilinearFormEJA)
355 Our inner product is "associative," which means the following for
356 a symmetric bilinear form::
358 sage: set_random_seed()
359 sage: J = random_eja()
360 sage: x,y,z = J.random_elements(3)
361 sage: (x*y).inner_product(z) == y.inner_product(x*z)
366 Ensure that this is the usual inner product for the algebras
369 sage: set_random_seed()
370 sage: J = HadamardEJA.random_instance()
371 sage: x,y = J.random_elements(2)
372 sage: actual = x.inner_product(y)
373 sage: expected = x.to_vector().inner_product(y.to_vector())
374 sage: actual == expected
377 Ensure that this is one-half of the trace inner-product in a
378 BilinearFormEJA that isn't just the reals (when ``n`` isn't
379 one). This is in Faraut and Koranyi, and also my "On the
382 sage: set_random_seed()
383 sage: J = BilinearFormEJA.random_instance()
384 sage: n = J.dimension()
385 sage: x = J.random_element()
386 sage: y = J.random_element()
387 sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
391 B
= self
._inner
_product
_matrix
392 return (B
*x
.to_vector()).inner_product(y
.to_vector())
395 def is_associative(self
):
397 Return whether or not this algebra's Jordan product is associative.
401 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
405 sage: J = ComplexHermitianEJA(3, field=QQ, orthonormalize=False)
406 sage: J.is_associative()
408 sage: x = sum(J.gens())
409 sage: A = x.subalgebra_generated_by(orthonormalize=False)
410 sage: A.is_associative()
414 return "Associative" in self
.category().axioms()
416 def _is_jordanian(self
):
418 Whether or not this algebra's multiplication table respects the
419 Jordan identity `(x^{2})(xy) = x(x^{2}y)`.
421 We only check one arrangement of `x` and `y`, so for a
422 ``True`` result to be truly true, you should also check
423 :meth:`is_commutative`. This method should of course always
424 return ``True``, unless this algebra was constructed with
425 ``check_axioms=False`` and passed an invalid multiplication table.
427 return all( (self
.gens()[i
]**2)*(self
.gens()[i
]*self
.gens()[j
])
429 (self
.gens()[i
])*((self
.gens()[i
]**2)*self
.gens()[j
])
430 for i
in range(self
.dimension())
431 for j
in range(self
.dimension()) )
433 def _inner_product_is_associative(self
):
435 Return whether or not this algebra's inner product `B` is
436 associative; that is, whether or not `B(xy,z) = B(x,yz)`.
438 This method should of course always return ``True``, unless
439 this algebra was constructed with ``check_axioms=False`` and
440 passed an invalid Jordan or inner-product.
443 # Used to check whether or not something is zero in an inexact
444 # ring. This number is sufficient to allow the construction of
445 # QuaternionHermitianEJA(2, field=RDF) with check_axioms=True.
448 for i
in range(self
.dimension()):
449 for j
in range(self
.dimension()):
450 for k
in range(self
.dimension()):
454 diff
= (x
*y
).inner_product(z
) - x
.inner_product(y
*z
)
456 if self
.base_ring().is_exact():
460 if diff
.abs() > epsilon
:
465 def _element_constructor_(self
, elt
):
467 Construct an element of this algebra from its vector or matrix
470 This gets called only after the parent element _call_ method
471 fails to find a coercion for the argument.
475 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
477 ....: RealSymmetricEJA)
481 The identity in `S^n` is converted to the identity in the EJA::
483 sage: J = RealSymmetricEJA(3)
484 sage: I = matrix.identity(QQ,3)
485 sage: J(I) == J.one()
488 This skew-symmetric matrix can't be represented in the EJA::
490 sage: J = RealSymmetricEJA(3)
491 sage: A = matrix(QQ,3, lambda i,j: i-j)
493 Traceback (most recent call last):
495 ValueError: not an element of this algebra
497 Tuples work as well, provided that the matrix basis for the
498 algebra consists of them::
500 sage: J1 = HadamardEJA(3)
501 sage: J2 = RealSymmetricEJA(2)
502 sage: J = cartesian_product([J1,J2])
503 sage: J( (J1.matrix_basis()[1], J2.matrix_basis()[2]) )
508 Ensure that we can convert any element of the two non-matrix
509 simple algebras (whose matrix representations are columns)
510 back and forth faithfully::
512 sage: set_random_seed()
513 sage: J = HadamardEJA.random_instance()
514 sage: x = J.random_element()
515 sage: J(x.to_vector().column()) == x
517 sage: J = JordanSpinEJA.random_instance()
518 sage: x = J.random_element()
519 sage: J(x.to_vector().column()) == x
522 We cannot coerce elements between algebras just because their
523 matrix representations are compatible::
525 sage: J1 = HadamardEJA(3)
526 sage: J2 = JordanSpinEJA(3)
528 Traceback (most recent call last):
530 ValueError: not an element of this algebra
532 Traceback (most recent call last):
534 ValueError: not an element of this algebra
537 msg
= "not an element of this algebra"
538 if elt
in self
.base_ring():
539 # Ensure that no base ring -> algebra coercion is performed
540 # by this method. There's some stupidity in sage that would
541 # otherwise propagate to this method; for example, sage thinks
542 # that the integer 3 belongs to the space of 2-by-2 matrices.
543 raise ValueError(msg
)
546 # Try to convert a vector into a column-matrix...
548 except (AttributeError, TypeError):
549 # and ignore failure, because we weren't really expecting
550 # a vector as an argument anyway.
553 if elt
not in self
.matrix_space():
554 raise ValueError(msg
)
556 # Thanks for nothing! Matrix spaces aren't vector spaces in
557 # Sage, so we have to figure out its matrix-basis coordinates
558 # ourselves. We use the basis space's ring instead of the
559 # element's ring because the basis space might be an algebraic
560 # closure whereas the base ring of the 3-by-3 identity matrix
561 # could be QQ instead of QQbar.
563 # And, we also have to handle Cartesian product bases (when
564 # the matrix basis consists of tuples) here. The "good news"
565 # is that we're already converting everything to long vectors,
566 # and that strategy works for tuples as well.
568 # We pass check=False because the matrix basis is "guaranteed"
569 # to be linearly independent... right? Ha ha.
571 V
= VectorSpace(self
.base_ring(), len(elt
))
572 W
= V
.span_of_basis( (V(_all2list(s
)) for s
in self
.matrix_basis()),
576 coords
= W
.coordinate_vector(V(elt
))
577 except ArithmeticError: # vector is not in free module
578 raise ValueError(msg
)
580 return self
.from_vector(coords
)
584 Return a string representation of ``self``.
588 sage: from mjo.eja.eja_algebra import JordanSpinEJA
592 Ensure that it says what we think it says::
594 sage: JordanSpinEJA(2, field=AA)
595 Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
596 sage: JordanSpinEJA(3, field=RDF)
597 Euclidean Jordan algebra of dimension 3 over Real Double Field
600 fmt
= "Euclidean Jordan algebra of dimension {} over {}"
601 return fmt
.format(self
.dimension(), self
.base_ring())
605 def characteristic_polynomial_of(self
):
607 Return the algebra's "characteristic polynomial of" function,
608 which is itself a multivariate polynomial that, when evaluated
609 at the coordinates of some algebra element, returns that
610 element's characteristic polynomial.
612 The resulting polynomial has `n+1` variables, where `n` is the
613 dimension of this algebra. The first `n` variables correspond to
614 the coordinates of an algebra element: when evaluated at the
615 coordinates of an algebra element with respect to a certain
616 basis, the result is a univariate polynomial (in the one
617 remaining variable ``t``), namely the characteristic polynomial
622 sage: from mjo.eja.eja_algebra import JordanSpinEJA, TrivialEJA
626 The characteristic polynomial in the spin algebra is given in
627 Alizadeh, Example 11.11::
629 sage: J = JordanSpinEJA(3)
630 sage: p = J.characteristic_polynomial_of(); p
631 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
632 sage: xvec = J.one().to_vector()
636 By definition, the characteristic polynomial is a monic
637 degree-zero polynomial in a rank-zero algebra. Note that
638 Cayley-Hamilton is indeed satisfied since the polynomial
639 ``1`` evaluates to the identity element of the algebra on
642 sage: J = TrivialEJA()
643 sage: J.characteristic_polynomial_of()
650 # The list of coefficient polynomials a_0, a_1, a_2, ..., a_(r-1).
651 a
= self
._charpoly
_coefficients
()
653 # We go to a bit of trouble here to reorder the
654 # indeterminates, so that it's easier to evaluate the
655 # characteristic polynomial at x's coordinates and get back
656 # something in terms of t, which is what we want.
657 S
= PolynomialRing(self
.base_ring(),'t')
661 S
= PolynomialRing(S
, R
.variable_names())
664 return (t
**r
+ sum( a
[k
]*(t
**k
) for k
in range(r
) ))
666 def coordinate_polynomial_ring(self
):
668 The multivariate polynomial ring in which this algebra's
669 :meth:`characteristic_polynomial_of` lives.
673 sage: from mjo.eja.eja_algebra import (HadamardEJA,
674 ....: RealSymmetricEJA)
678 sage: J = HadamardEJA(2)
679 sage: J.coordinate_polynomial_ring()
680 Multivariate Polynomial Ring in X1, X2...
681 sage: J = RealSymmetricEJA(3,field=QQ,orthonormalize=False)
682 sage: J.coordinate_polynomial_ring()
683 Multivariate Polynomial Ring in X1, X2, X3, X4, X5, X6...
686 var_names
= tuple( "X%d" % z
for z
in range(1, self
.dimension()+1) )
687 return PolynomialRing(self
.base_ring(), var_names
)
689 def inner_product(self
, x
, y
):
691 The inner product associated with this Euclidean Jordan algebra.
693 Defaults to the trace inner product, but can be overridden by
694 subclasses if they are sure that the necessary properties are
699 sage: from mjo.eja.eja_algebra import (random_eja,
701 ....: BilinearFormEJA)
705 Our inner product is "associative," which means the following for
706 a symmetric bilinear form::
708 sage: set_random_seed()
709 sage: J = random_eja()
710 sage: x,y,z = J.random_elements(3)
711 sage: (x*y).inner_product(z) == y.inner_product(x*z)
716 Ensure that this is the usual inner product for the algebras
719 sage: set_random_seed()
720 sage: J = HadamardEJA.random_instance()
721 sage: x,y = J.random_elements(2)
722 sage: actual = x.inner_product(y)
723 sage: expected = x.to_vector().inner_product(y.to_vector())
724 sage: actual == expected
727 Ensure that this is one-half of the trace inner-product in a
728 BilinearFormEJA that isn't just the reals (when ``n`` isn't
729 one). This is in Faraut and Koranyi, and also my "On the
732 sage: set_random_seed()
733 sage: J = BilinearFormEJA.random_instance()
734 sage: n = J.dimension()
735 sage: x = J.random_element()
736 sage: y = J.random_element()
737 sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
740 B
= self
._inner
_product
_matrix
741 return (B
*x
.to_vector()).inner_product(y
.to_vector())
744 def is_trivial(self
):
746 Return whether or not this algebra is trivial.
748 A trivial algebra contains only the zero element.
752 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
757 sage: J = ComplexHermitianEJA(3)
763 sage: J = TrivialEJA()
768 return self
.dimension() == 0
771 def multiplication_table(self
):
773 Return a visual representation of this algebra's multiplication
774 table (on basis elements).
778 sage: from mjo.eja.eja_algebra import JordanSpinEJA
782 sage: J = JordanSpinEJA(4)
783 sage: J.multiplication_table()
784 +----++----+----+----+----+
785 | * || e0 | e1 | e2 | e3 |
786 +====++====+====+====+====+
787 | e0 || e0 | e1 | e2 | e3 |
788 +----++----+----+----+----+
789 | e1 || e1 | e0 | 0 | 0 |
790 +----++----+----+----+----+
791 | e2 || e2 | 0 | e0 | 0 |
792 +----++----+----+----+----+
793 | e3 || e3 | 0 | 0 | e0 |
794 +----++----+----+----+----+
798 # Prepend the header row.
799 M
= [["*"] + list(self
.gens())]
801 # And to each subsequent row, prepend an entry that belongs to
802 # the left-side "header column."
803 M
+= [ [self
.gens()[i
]] + [ self
.product_on_basis(i
,j
)
807 return table(M
, header_row
=True, header_column
=True, frame
=True)
810 def matrix_basis(self
):
812 Return an (often more natural) representation of this algebras
813 basis as an ordered tuple of matrices.
815 Every finite-dimensional Euclidean Jordan Algebra is a, up to
816 Jordan isomorphism, a direct sum of five simple
817 algebras---four of which comprise Hermitian matrices. And the
818 last type of algebra can of course be thought of as `n`-by-`1`
819 column matrices (ambiguusly called column vectors) to avoid
820 special cases. As a result, matrices (and column vectors) are
821 a natural representation format for Euclidean Jordan algebra
824 But, when we construct an algebra from a basis of matrices,
825 those matrix representations are lost in favor of coordinate
826 vectors *with respect to* that basis. We could eventually
827 convert back if we tried hard enough, but having the original
828 representations handy is valuable enough that we simply store
829 them and return them from this method.
831 Why implement this for non-matrix algebras? Avoiding special
832 cases for the :class:`BilinearFormEJA` pays with simplicity in
833 its own right. But mainly, we would like to be able to assume
834 that elements of a :class:`CartesianProductEJA` can be displayed
835 nicely, without having to have special classes for direct sums
836 one of whose components was a matrix algebra.
840 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
841 ....: RealSymmetricEJA)
845 sage: J = RealSymmetricEJA(2)
847 Finite family {0: e0, 1: e1, 2: e2}
848 sage: J.matrix_basis()
850 [1 0] [ 0 0.7071067811865475?] [0 0]
851 [0 0], [0.7071067811865475? 0], [0 1]
856 sage: J = JordanSpinEJA(2)
858 Finite family {0: e0, 1: e1}
859 sage: J.matrix_basis()
865 return self
._matrix
_basis
868 def matrix_space(self
):
870 Return the matrix space in which this algebra's elements live, if
871 we think of them as matrices (including column vectors of the
874 "By default" this will be an `n`-by-`1` column-matrix space,
875 except when the algebra is trivial. There it's `n`-by-`n`
876 (where `n` is zero), to ensure that two elements of the matrix
877 space (empty matrices) can be multiplied. For algebras of
878 matrices, this returns the space in which their
879 real embeddings live.
883 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
885 ....: QuaternionHermitianEJA,
890 By default, the matrix representation is just a column-matrix
891 equivalent to the vector representation::
893 sage: J = JordanSpinEJA(3)
894 sage: J.matrix_space()
895 Full MatrixSpace of 3 by 1 dense matrices over Algebraic
898 The matrix representation in the trivial algebra is
899 zero-by-zero instead of the usual `n`-by-one::
901 sage: J = TrivialEJA()
902 sage: J.matrix_space()
903 Full MatrixSpace of 0 by 0 dense matrices over Algebraic
906 The matrix space for complex/quaternion Hermitian matrix EJA
907 is the space in which their real-embeddings live, not the
908 original complex/quaternion matrix space::
910 sage: J = ComplexHermitianEJA(2,field=QQ,orthonormalize=False)
911 sage: J.matrix_space()
912 Full MatrixSpace of 4 by 4 dense matrices over Rational Field
913 sage: J = QuaternionHermitianEJA(1,field=QQ,orthonormalize=False)
914 sage: J.matrix_space()
915 Full MatrixSpace of 4 by 4 dense matrices over Rational Field
918 if self
.is_trivial():
919 return MatrixSpace(self
.base_ring(), 0)
921 return self
.matrix_basis()[0].parent()
927 Return the unit element of this algebra.
931 sage: from mjo.eja.eja_algebra import (HadamardEJA,
936 We can compute unit element in the Hadamard EJA::
938 sage: J = HadamardEJA(5)
940 e0 + e1 + e2 + e3 + e4
942 The unit element in the Hadamard EJA is inherited in the
943 subalgebras generated by its elements::
945 sage: J = HadamardEJA(5)
947 e0 + e1 + e2 + e3 + e4
948 sage: x = sum(J.gens())
949 sage: A = x.subalgebra_generated_by(orthonormalize=False)
952 sage: A.one().superalgebra_element()
953 e0 + e1 + e2 + e3 + e4
957 The identity element acts like the identity, regardless of
958 whether or not we orthonormalize::
960 sage: set_random_seed()
961 sage: J = random_eja()
962 sage: x = J.random_element()
963 sage: J.one()*x == x and x*J.one() == x
965 sage: A = x.subalgebra_generated_by()
966 sage: y = A.random_element()
967 sage: A.one()*y == y and y*A.one() == y
972 sage: set_random_seed()
973 sage: J = random_eja(field=QQ, orthonormalize=False)
974 sage: x = J.random_element()
975 sage: J.one()*x == x and x*J.one() == x
977 sage: A = x.subalgebra_generated_by(orthonormalize=False)
978 sage: y = A.random_element()
979 sage: A.one()*y == y and y*A.one() == y
982 The matrix of the unit element's operator is the identity,
983 regardless of the base field and whether or not we
986 sage: set_random_seed()
987 sage: J = random_eja()
988 sage: actual = J.one().operator().matrix()
989 sage: expected = matrix.identity(J.base_ring(), J.dimension())
990 sage: actual == expected
992 sage: x = J.random_element()
993 sage: A = x.subalgebra_generated_by()
994 sage: actual = A.one().operator().matrix()
995 sage: expected = matrix.identity(A.base_ring(), A.dimension())
996 sage: actual == expected
1001 sage: set_random_seed()
1002 sage: J = random_eja(field=QQ, orthonormalize=False)
1003 sage: actual = J.one().operator().matrix()
1004 sage: expected = matrix.identity(J.base_ring(), J.dimension())
1005 sage: actual == expected
1007 sage: x = J.random_element()
1008 sage: A = x.subalgebra_generated_by(orthonormalize=False)
1009 sage: actual = A.one().operator().matrix()
1010 sage: expected = matrix.identity(A.base_ring(), A.dimension())
1011 sage: actual == expected
1014 Ensure that the cached unit element (often precomputed by
1015 hand) agrees with the computed one::
1017 sage: set_random_seed()
1018 sage: J = random_eja()
1019 sage: cached = J.one()
1020 sage: J.one.clear_cache()
1021 sage: J.one() == cached
1026 sage: set_random_seed()
1027 sage: J = random_eja(field=QQ, orthonormalize=False)
1028 sage: cached = J.one()
1029 sage: J.one.clear_cache()
1030 sage: J.one() == cached
1034 # We can brute-force compute the matrices of the operators
1035 # that correspond to the basis elements of this algebra.
1036 # If some linear combination of those basis elements is the
1037 # algebra identity, then the same linear combination of
1038 # their matrices has to be the identity matrix.
1040 # Of course, matrices aren't vectors in sage, so we have to
1041 # appeal to the "long vectors" isometry.
1042 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
1044 # Now we use basic linear algebra to find the coefficients,
1045 # of the matrices-as-vectors-linear-combination, which should
1046 # work for the original algebra basis too.
1047 A
= matrix(self
.base_ring(), oper_vecs
)
1049 # We used the isometry on the left-hand side already, but we
1050 # still need to do it for the right-hand side. Recall that we
1051 # wanted something that summed to the identity matrix.
1052 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
1054 # Now if there's an identity element in the algebra, this
1055 # should work. We solve on the left to avoid having to
1056 # transpose the matrix "A".
1057 return self
.from_vector(A
.solve_left(b
))
1060 def peirce_decomposition(self
, c
):
1062 The Peirce decomposition of this algebra relative to the
1065 In the future, this can be extended to a complete system of
1066 orthogonal idempotents.
1070 - ``c`` -- an idempotent of this algebra.
1074 A triple (J0, J5, J1) containing two subalgebras and one subspace
1077 - ``J0`` -- the algebra on the eigenspace of ``c.operator()``
1078 corresponding to the eigenvalue zero.
1080 - ``J5`` -- the eigenspace (NOT a subalgebra) of ``c.operator()``
1081 corresponding to the eigenvalue one-half.
1083 - ``J1`` -- the algebra on the eigenspace of ``c.operator()``
1084 corresponding to the eigenvalue one.
1086 These are the only possible eigenspaces for that operator, and this
1087 algebra is a direct sum of them. The spaces ``J0`` and ``J1`` are
1088 orthogonal, and are subalgebras of this algebra with the appropriate
1093 sage: from mjo.eja.eja_algebra import random_eja, RealSymmetricEJA
1097 The canonical example comes from the symmetric matrices, which
1098 decompose into diagonal and off-diagonal parts::
1100 sage: J = RealSymmetricEJA(3)
1101 sage: C = matrix(QQ, [ [1,0,0],
1105 sage: J0,J5,J1 = J.peirce_decomposition(c)
1107 Euclidean Jordan algebra of dimension 1...
1109 Vector space of degree 6 and dimension 2...
1111 Euclidean Jordan algebra of dimension 3...
1112 sage: J0.one().to_matrix()
1116 sage: orig_df = AA.options.display_format
1117 sage: AA.options.display_format = 'radical'
1118 sage: J.from_vector(J5.basis()[0]).to_matrix()
1122 sage: J.from_vector(J5.basis()[1]).to_matrix()
1126 sage: AA.options.display_format = orig_df
1127 sage: J1.one().to_matrix()
1134 Every algebra decomposes trivially with respect to its identity
1137 sage: set_random_seed()
1138 sage: J = random_eja()
1139 sage: J0,J5,J1 = J.peirce_decomposition(J.one())
1140 sage: J0.dimension() == 0 and J5.dimension() == 0
1142 sage: J1.superalgebra() == J and J1.dimension() == J.dimension()
1145 The decomposition is into eigenspaces, and its components are
1146 therefore necessarily orthogonal. Moreover, the identity
1147 elements in the two subalgebras are the projections onto their
1148 respective subspaces of the superalgebra's identity element::
1150 sage: set_random_seed()
1151 sage: J = random_eja()
1152 sage: x = J.random_element()
1153 sage: if not J.is_trivial():
1154 ....: while x.is_nilpotent():
1155 ....: x = J.random_element()
1156 sage: c = x.subalgebra_idempotent()
1157 sage: J0,J5,J1 = J.peirce_decomposition(c)
1159 sage: for (w,y,z) in zip(J0.basis(), J5.basis(), J1.basis()):
1160 ....: w = w.superalgebra_element()
1161 ....: y = J.from_vector(y)
1162 ....: z = z.superalgebra_element()
1163 ....: ipsum += w.inner_product(y).abs()
1164 ....: ipsum += w.inner_product(z).abs()
1165 ....: ipsum += y.inner_product(z).abs()
1168 sage: J1(c) == J1.one()
1170 sage: J0(J.one() - c) == J0.one()
1174 if not c
.is_idempotent():
1175 raise ValueError("element is not idempotent: %s" % c
)
1177 # Default these to what they should be if they turn out to be
1178 # trivial, because eigenspaces_left() won't return eigenvalues
1179 # corresponding to trivial spaces (e.g. it returns only the
1180 # eigenspace corresponding to lambda=1 if you take the
1181 # decomposition relative to the identity element).
1182 trivial
= self
.subalgebra(())
1183 J0
= trivial
# eigenvalue zero
1184 J5
= VectorSpace(self
.base_ring(), 0) # eigenvalue one-half
1185 J1
= trivial
# eigenvalue one
1187 for (eigval
, eigspace
) in c
.operator().matrix().right_eigenspaces():
1188 if eigval
== ~
(self
.base_ring()(2)):
1191 gens
= tuple( self
.from_vector(b
) for b
in eigspace
.basis() )
1192 subalg
= self
.subalgebra(gens
, check_axioms
=False)
1198 raise ValueError("unexpected eigenvalue: %s" % eigval
)
1203 def random_element(self
, thorough
=False):
1205 Return a random element of this algebra.
1207 Our algebra superclass method only returns a linear
1208 combination of at most two basis elements. We instead
1209 want the vector space "random element" method that
1210 returns a more diverse selection.
1214 - ``thorough`` -- (boolean; default False) whether or not we
1215 should generate irrational coefficients for the random
1216 element when our base ring is irrational; this slows the
1217 algebra operations to a crawl, but any truly random method
1221 # For a general base ring... maybe we can trust this to do the
1222 # right thing? Unlikely, but.
1223 V
= self
.vector_space()
1224 v
= V
.random_element()
1226 if self
.base_ring() is AA
:
1227 # The "random element" method of the algebraic reals is
1228 # stupid at the moment, and only returns integers between
1229 # -2 and 2, inclusive:
1231 # https://trac.sagemath.org/ticket/30875
1233 # Instead, we implement our own "random vector" method,
1234 # and then coerce that into the algebra. We use the vector
1235 # space degree here instead of the dimension because a
1236 # subalgebra could (for example) be spanned by only two
1237 # vectors, each with five coordinates. We need to
1238 # generate all five coordinates.
1240 v
*= QQbar
.random_element().real()
1242 v
*= QQ
.random_element()
1244 return self
.from_vector(V
.coordinate_vector(v
))
1246 def random_elements(self
, count
, thorough
=False):
1248 Return ``count`` random elements as a tuple.
1252 - ``thorough`` -- (boolean; default False) whether or not we
1253 should generate irrational coefficients for the random
1254 elements when our base ring is irrational; this slows the
1255 algebra operations to a crawl, but any truly random method
1260 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1264 sage: J = JordanSpinEJA(3)
1265 sage: x,y,z = J.random_elements(3)
1266 sage: all( [ x in J, y in J, z in J ])
1268 sage: len( J.random_elements(10) ) == 10
1272 return tuple( self
.random_element(thorough
)
1273 for idx
in range(count
) )
1277 def _charpoly_coefficients(self
):
1279 The `r` polynomial coefficients of the "characteristic polynomial
1284 sage: from mjo.eja.eja_algebra import random_eja
1288 The theory shows that these are all homogeneous polynomials of
1291 sage: set_random_seed()
1292 sage: J = random_eja()
1293 sage: all(p.is_homogeneous() for p in J._charpoly_coefficients())
1297 n
= self
.dimension()
1298 R
= self
.coordinate_polynomial_ring()
1300 F
= R
.fraction_field()
1303 # From a result in my book, these are the entries of the
1304 # basis representation of L_x.
1305 return sum( vars[k
]*self
.gens()[k
].operator().matrix()[i
,j
]
1308 L_x
= matrix(F
, n
, n
, L_x_i_j
)
1311 if self
.rank
.is_in_cache():
1313 # There's no need to pad the system with redundant
1314 # columns if we *know* they'll be redundant.
1317 # Compute an extra power in case the rank is equal to
1318 # the dimension (otherwise, we would stop at x^(r-1)).
1319 x_powers
= [ (L_x
**k
)*self
.one().to_vector()
1320 for k
in range(n
+1) ]
1321 A
= matrix
.column(F
, x_powers
[:n
])
1322 AE
= A
.extended_echelon_form()
1329 # The theory says that only the first "r" coefficients are
1330 # nonzero, and they actually live in the original polynomial
1331 # ring and not the fraction field. We negate them because in
1332 # the actual characteristic polynomial, they get moved to the
1333 # other side where x^r lives. We don't bother to trim A_rref
1334 # down to a square matrix and solve the resulting system,
1335 # because the upper-left r-by-r portion of A_rref is
1336 # guaranteed to be the identity matrix, so e.g.
1338 # A_rref.solve_right(Y)
1340 # would just be returning Y.
1341 return (-E
*b
)[:r
].change_ring(R
)
1346 Return the rank of this EJA.
1348 This is a cached method because we know the rank a priori for
1349 all of the algebras we can construct. Thus we can avoid the
1350 expensive ``_charpoly_coefficients()`` call unless we truly
1351 need to compute the whole characteristic polynomial.
1355 sage: from mjo.eja.eja_algebra import (HadamardEJA,
1356 ....: JordanSpinEJA,
1357 ....: RealSymmetricEJA,
1358 ....: ComplexHermitianEJA,
1359 ....: QuaternionHermitianEJA,
1364 The rank of the Jordan spin algebra is always two::
1366 sage: JordanSpinEJA(2).rank()
1368 sage: JordanSpinEJA(3).rank()
1370 sage: JordanSpinEJA(4).rank()
1373 The rank of the `n`-by-`n` Hermitian real, complex, or
1374 quaternion matrices is `n`::
1376 sage: RealSymmetricEJA(4).rank()
1378 sage: ComplexHermitianEJA(3).rank()
1380 sage: QuaternionHermitianEJA(2).rank()
1385 Ensure that every EJA that we know how to construct has a
1386 positive integer rank, unless the algebra is trivial in
1387 which case its rank will be zero::
1389 sage: set_random_seed()
1390 sage: J = random_eja()
1394 sage: r > 0 or (r == 0 and J.is_trivial())
1397 Ensure that computing the rank actually works, since the ranks
1398 of all simple algebras are known and will be cached by default::
1400 sage: set_random_seed() # long time
1401 sage: J = random_eja() # long time
1402 sage: cached = J.rank() # long time
1403 sage: J.rank.clear_cache() # long time
1404 sage: J.rank() == cached # long time
1408 return len(self
._charpoly
_coefficients
())
1411 def subalgebra(self
, basis
, **kwargs
):
1413 Create a subalgebra of this algebra from the given basis.
1415 from mjo
.eja
.eja_subalgebra
import FiniteDimensionalEJASubalgebra
1416 return FiniteDimensionalEJASubalgebra(self
, basis
, **kwargs
)
1419 def vector_space(self
):
1421 Return the vector space that underlies this algebra.
1425 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1429 sage: J = RealSymmetricEJA(2)
1430 sage: J.vector_space()
1431 Vector space of dimension 3 over...
1434 return self
.zero().to_vector().parent().ambient_vector_space()
1438 class RationalBasisEJA(FiniteDimensionalEJA
):
1440 New class for algebras whose supplied basis elements have all rational entries.
1444 sage: from mjo.eja.eja_algebra import BilinearFormEJA
1448 The supplied basis is orthonormalized by default::
1450 sage: B = matrix(QQ, [[1, 0, 0], [0, 25, -32], [0, -32, 41]])
1451 sage: J = BilinearFormEJA(B)
1452 sage: J.matrix_basis()
1469 # Abuse the check_field parameter to check that the entries of
1470 # out basis (in ambient coordinates) are in the field QQ.
1471 if not all( all(b_i
in QQ
for b_i
in b
.list()) for b
in basis
):
1472 raise TypeError("basis not rational")
1474 self
._rational
_algebra
= None
1476 # There's no point in constructing the extra algebra if this
1477 # one is already rational.
1479 # Note: the same Jordan and inner-products work here,
1480 # because they are necessarily defined with respect to
1481 # ambient coordinates and not any particular basis.
1482 self
._rational
_algebra
= FiniteDimensionalEJA(
1487 orthonormalize
=False,
1491 super().__init
__(basis
,
1495 check_field
=check_field
,
1499 def _charpoly_coefficients(self
):
1503 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
1504 ....: JordanSpinEJA)
1508 The base ring of the resulting polynomial coefficients is what
1509 it should be, and not the rationals (unless the algebra was
1510 already over the rationals)::
1512 sage: J = JordanSpinEJA(3)
1513 sage: J._charpoly_coefficients()
1514 (X1^2 - X2^2 - X3^2, -2*X1)
1515 sage: a0 = J._charpoly_coefficients()[0]
1517 Algebraic Real Field
1518 sage: a0.base_ring()
1519 Algebraic Real Field
1522 if self
._rational
_algebra
is None:
1523 # There's no need to construct *another* algebra over the
1524 # rationals if this one is already over the
1525 # rationals. Likewise, if we never orthonormalized our
1526 # basis, we might as well just use the given one.
1527 return super()._charpoly
_coefficients
()
1529 # Do the computation over the rationals. The answer will be
1530 # the same, because all we've done is a change of basis.
1531 # Then, change back from QQ to our real base ring
1532 a
= ( a_i
.change_ring(self
.base_ring())
1533 for a_i
in self
._rational
_algebra
._charpoly
_coefficients
() )
1535 if self
._deortho
_matrix
is None:
1536 # This can happen if our base ring was, say, AA and we
1537 # chose not to (or didn't need to) orthonormalize. It's
1538 # still faster to do the computations over QQ even if
1539 # the numbers in the boxes stay the same.
1542 # Otherwise, convert the coordinate variables back to the
1543 # deorthonormalized ones.
1544 R
= self
.coordinate_polynomial_ring()
1545 from sage
.modules
.free_module_element
import vector
1546 X
= vector(R
, R
.gens())
1547 BX
= self
._deortho
_matrix
*X
1549 subs_dict
= { X[i]: BX[i] for i in range(len(X)) }
1550 return tuple( a_i
.subs(subs_dict
) for a_i
in a
)
1552 class ConcreteEJA(RationalBasisEJA
):
1554 A class for the Euclidean Jordan algebras that we know by name.
1556 These are the Jordan algebras whose basis, multiplication table,
1557 rank, and so on are known a priori. More to the point, they are
1558 the Euclidean Jordan algebras for which we are able to conjure up
1559 a "random instance."
1563 sage: from mjo.eja.eja_algebra import ConcreteEJA
1567 Our basis is normalized with respect to the algebra's inner
1568 product, unless we specify otherwise::
1570 sage: set_random_seed()
1571 sage: J = ConcreteEJA.random_instance()
1572 sage: all( b.norm() == 1 for b in J.gens() )
1575 Since our basis is orthonormal with respect to the algebra's inner
1576 product, and since we know that this algebra is an EJA, any
1577 left-multiplication operator's matrix will be symmetric because
1578 natural->EJA basis representation is an isometry and within the
1579 EJA the operator is self-adjoint by the Jordan axiom::
1581 sage: set_random_seed()
1582 sage: J = ConcreteEJA.random_instance()
1583 sage: x = J.random_element()
1584 sage: x.operator().is_self_adjoint()
1589 def _max_random_instance_size():
1591 Return an integer "size" that is an upper bound on the size of
1592 this algebra when it is used in a random test
1593 case. Unfortunately, the term "size" is ambiguous -- when
1594 dealing with `R^n` under either the Hadamard or Jordan spin
1595 product, the "size" refers to the dimension `n`. When dealing
1596 with a matrix algebra (real symmetric or complex/quaternion
1597 Hermitian), it refers to the size of the matrix, which is far
1598 less than the dimension of the underlying vector space.
1600 This method must be implemented in each subclass.
1602 raise NotImplementedError
1605 def random_instance(cls
, *args
, **kwargs
):
1607 Return a random instance of this type of algebra.
1609 This method should be implemented in each subclass.
1611 from sage
.misc
.prandom
import choice
1612 eja_class
= choice(cls
.__subclasses
__())
1614 # These all bubble up to the RationalBasisEJA superclass
1615 # constructor, so any (kw)args valid there are also valid
1617 return eja_class
.random_instance(*args
, **kwargs
)
1622 def dimension_over_reals():
1624 The dimension of this matrix's base ring over the reals.
1626 The reals are dimension one over themselves, obviously; that's
1627 just `\mathbb{R}^{1}`. Likewise, the complex numbers `a + bi`
1628 have dimension two. Finally, the quaternions have dimension
1629 four over the reals.
1631 This is used to determine the size of the matrix returned from
1632 :meth:`real_embed`, among other things.
1634 raise NotImplementedError
1637 def real_embed(cls
,M
):
1639 Embed the matrix ``M`` into a space of real matrices.
1641 The matrix ``M`` can have entries in any field at the moment:
1642 the real numbers, complex numbers, or quaternions. And although
1643 they are not a field, we can probably support octonions at some
1644 point, too. This function returns a real matrix that "acts like"
1645 the original with respect to matrix multiplication; i.e.
1647 real_embed(M*N) = real_embed(M)*real_embed(N)
1650 if M
.ncols() != M
.nrows():
1651 raise ValueError("the matrix 'M' must be square")
1656 def real_unembed(cls
,M
):
1658 The inverse of :meth:`real_embed`.
1660 if M
.ncols() != M
.nrows():
1661 raise ValueError("the matrix 'M' must be square")
1662 if not ZZ(M
.nrows()).mod(cls
.dimension_over_reals()).is_zero():
1663 raise ValueError("the matrix 'M' must be a real embedding")
1667 def jordan_product(X
,Y
):
1668 return (X
*Y
+ Y
*X
)/2
1671 def trace_inner_product(cls
,X
,Y
):
1673 Compute the trace inner-product of two real-embeddings.
1677 sage: from mjo.eja.eja_algebra import (RealSymmetricEJA,
1678 ....: ComplexHermitianEJA,
1679 ....: QuaternionHermitianEJA)
1683 This gives the same answer as it would if we computed the trace
1684 from the unembedded (original) matrices::
1686 sage: set_random_seed()
1687 sage: J = RealSymmetricEJA.random_instance()
1688 sage: x,y = J.random_elements(2)
1689 sage: Xe = x.to_matrix()
1690 sage: Ye = y.to_matrix()
1691 sage: X = J.real_unembed(Xe)
1692 sage: Y = J.real_unembed(Ye)
1693 sage: expected = (X*Y).trace()
1694 sage: actual = J.trace_inner_product(Xe,Ye)
1695 sage: actual == expected
1700 sage: set_random_seed()
1701 sage: J = ComplexHermitianEJA.random_instance()
1702 sage: x,y = J.random_elements(2)
1703 sage: Xe = x.to_matrix()
1704 sage: Ye = y.to_matrix()
1705 sage: X = J.real_unembed(Xe)
1706 sage: Y = J.real_unembed(Ye)
1707 sage: expected = (X*Y).trace().real()
1708 sage: actual = J.trace_inner_product(Xe,Ye)
1709 sage: actual == expected
1714 sage: set_random_seed()
1715 sage: J = QuaternionHermitianEJA.random_instance()
1716 sage: x,y = J.random_elements(2)
1717 sage: Xe = x.to_matrix()
1718 sage: Ye = y.to_matrix()
1719 sage: X = J.real_unembed(Xe)
1720 sage: Y = J.real_unembed(Ye)
1721 sage: expected = (X*Y).trace().coefficient_tuple()[0]
1722 sage: actual = J.trace_inner_product(Xe,Ye)
1723 sage: actual == expected
1727 Xu
= cls
.real_unembed(X
)
1728 Yu
= cls
.real_unembed(Y
)
1729 tr
= (Xu
*Yu
).trace()
1732 # Works in QQ, AA, RDF, et cetera.
1734 except AttributeError:
1735 # A quaternion doesn't have a real() method, but does
1736 # have coefficient_tuple() method that returns the
1737 # coefficients of 1, i, j, and k -- in that order.
1738 return tr
.coefficient_tuple()[0]
1741 class RealMatrixEJA(MatrixEJA
):
1743 def dimension_over_reals():
1747 class RealSymmetricEJA(ConcreteEJA
, RealMatrixEJA
):
1749 The rank-n simple EJA consisting of real symmetric n-by-n
1750 matrices, the usual symmetric Jordan product, and the trace inner
1751 product. It has dimension `(n^2 + n)/2` over the reals.
1755 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1759 sage: J = RealSymmetricEJA(2)
1760 sage: e0, e1, e2 = J.gens()
1768 In theory, our "field" can be any subfield of the reals::
1770 sage: RealSymmetricEJA(2, field=RDF)
1771 Euclidean Jordan algebra of dimension 3 over Real Double Field
1772 sage: RealSymmetricEJA(2, field=RR)
1773 Euclidean Jordan algebra of dimension 3 over Real Field with
1774 53 bits of precision
1778 The dimension of this algebra is `(n^2 + n) / 2`::
1780 sage: set_random_seed()
1781 sage: n_max = RealSymmetricEJA._max_random_instance_size()
1782 sage: n = ZZ.random_element(1, n_max)
1783 sage: J = RealSymmetricEJA(n)
1784 sage: J.dimension() == (n^2 + n)/2
1787 The Jordan multiplication is what we think it is::
1789 sage: set_random_seed()
1790 sage: J = RealSymmetricEJA.random_instance()
1791 sage: x,y = J.random_elements(2)
1792 sage: actual = (x*y).to_matrix()
1793 sage: X = x.to_matrix()
1794 sage: Y = y.to_matrix()
1795 sage: expected = (X*Y + Y*X)/2
1796 sage: actual == expected
1798 sage: J(expected) == x*y
1801 We can change the generator prefix::
1803 sage: RealSymmetricEJA(3, prefix='q').gens()
1804 (q0, q1, q2, q3, q4, q5)
1806 We can construct the (trivial) algebra of rank zero::
1808 sage: RealSymmetricEJA(0)
1809 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1813 def _denormalized_basis(cls
, n
):
1815 Return a basis for the space of real symmetric n-by-n matrices.
1819 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1823 sage: set_random_seed()
1824 sage: n = ZZ.random_element(1,5)
1825 sage: B = RealSymmetricEJA._denormalized_basis(n)
1826 sage: all( M.is_symmetric() for M in B)
1830 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1834 for j
in range(i
+1):
1835 Eij
= matrix(ZZ
, n
, lambda k
,l
: k
==i
and l
==j
)
1839 Sij
= Eij
+ Eij
.transpose()
1845 def _max_random_instance_size():
1846 return 4 # Dimension 10
1849 def random_instance(cls
, **kwargs
):
1851 Return a random instance of this type of algebra.
1853 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1854 return cls(n
, **kwargs
)
1856 def __init__(self
, n
, **kwargs
):
1857 # We know this is a valid EJA, but will double-check
1858 # if the user passes check_axioms=True.
1859 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
1861 super(RealSymmetricEJA
, self
).__init
__(self
._denormalized
_basis
(n
),
1862 self
.jordan_product
,
1863 self
.trace_inner_product
,
1866 # TODO: this could be factored out somehow, but is left here
1867 # because the MatrixEJA is not presently a subclass of the
1868 # FDEJA class that defines rank() and one().
1869 self
.rank
.set_cache(n
)
1870 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
1871 self
.one
.set_cache(self(idV
))
1875 class ComplexMatrixEJA(MatrixEJA
):
1876 # A manual dictionary-cache for the complex_extension() method,
1877 # since apparently @classmethods can't also be @cached_methods.
1878 _complex_extension
= {}
1881 def complex_extension(cls
,field
):
1883 The complex field that we embed/unembed, as an extension
1884 of the given ``field``.
1886 if field
in cls
._complex
_extension
:
1887 return cls
._complex
_extension
[field
]
1889 # Sage doesn't know how to adjoin the complex "i" (the root of
1890 # x^2 + 1) to a field in a general way. Here, we just enumerate
1891 # all of the cases that I have cared to support so far.
1893 # Sage doesn't know how to embed AA into QQbar, i.e. how
1894 # to adjoin sqrt(-1) to AA.
1896 elif not field
.is_exact():
1898 F
= field
.complex_field()
1900 # Works for QQ and... maybe some other fields.
1901 R
= PolynomialRing(field
, 'z')
1903 F
= field
.extension(z
**2 + 1, 'I', embedding
=CLF(-1).sqrt())
1905 cls
._complex
_extension
[field
] = F
1909 def dimension_over_reals():
1913 def real_embed(cls
,M
):
1915 Embed the n-by-n complex matrix ``M`` into the space of real
1916 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1917 bi` to the block matrix ``[[a,b],[-b,a]]``.
1921 sage: from mjo.eja.eja_algebra import ComplexMatrixEJA
1925 sage: F = QuadraticField(-1, 'I')
1926 sage: x1 = F(4 - 2*i)
1927 sage: x2 = F(1 + 2*i)
1930 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1931 sage: ComplexMatrixEJA.real_embed(M)
1940 Embedding is a homomorphism (isomorphism, in fact)::
1942 sage: set_random_seed()
1943 sage: n = ZZ.random_element(3)
1944 sage: F = QuadraticField(-1, 'I')
1945 sage: X = random_matrix(F, n)
1946 sage: Y = random_matrix(F, n)
1947 sage: Xe = ComplexMatrixEJA.real_embed(X)
1948 sage: Ye = ComplexMatrixEJA.real_embed(Y)
1949 sage: XYe = ComplexMatrixEJA.real_embed(X*Y)
1954 super(ComplexMatrixEJA
,cls
).real_embed(M
)
1957 # We don't need any adjoined elements...
1958 field
= M
.base_ring().base_ring()
1964 blocks
.append(matrix(field
, 2, [ [ a
, b
],
1967 return matrix
.block(field
, n
, blocks
)
1971 def real_unembed(cls
,M
):
1973 The inverse of _embed_complex_matrix().
1977 sage: from mjo.eja.eja_algebra import ComplexMatrixEJA
1981 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1982 ....: [-2, 1, -4, 3],
1983 ....: [ 9, 10, 11, 12],
1984 ....: [-10, 9, -12, 11] ])
1985 sage: ComplexMatrixEJA.real_unembed(A)
1987 [ 10*I + 9 12*I + 11]
1991 Unembedding is the inverse of embedding::
1993 sage: set_random_seed()
1994 sage: F = QuadraticField(-1, 'I')
1995 sage: M = random_matrix(F, 3)
1996 sage: Me = ComplexMatrixEJA.real_embed(M)
1997 sage: ComplexMatrixEJA.real_unembed(Me) == M
2001 super(ComplexMatrixEJA
,cls
).real_unembed(M
)
2003 d
= cls
.dimension_over_reals()
2004 F
= cls
.complex_extension(M
.base_ring())
2007 # Go top-left to bottom-right (reading order), converting every
2008 # 2-by-2 block we see to a single complex element.
2010 for k
in range(n
/d
):
2011 for j
in range(n
/d
):
2012 submat
= M
[d
*k
:d
*k
+d
,d
*j
:d
*j
+d
]
2013 if submat
[0,0] != submat
[1,1]:
2014 raise ValueError('bad on-diagonal submatrix')
2015 if submat
[0,1] != -submat
[1,0]:
2016 raise ValueError('bad off-diagonal submatrix')
2017 z
= submat
[0,0] + submat
[0,1]*i
2020 return matrix(F
, n
/d
, elements
)
2023 class ComplexHermitianEJA(ConcreteEJA
, ComplexMatrixEJA
):
2025 The rank-n simple EJA consisting of complex Hermitian n-by-n
2026 matrices over the real numbers, the usual symmetric Jordan product,
2027 and the real-part-of-trace inner product. It has dimension `n^2` over
2032 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
2036 In theory, our "field" can be any subfield of the reals::
2038 sage: ComplexHermitianEJA(2, field=RDF)
2039 Euclidean Jordan algebra of dimension 4 over Real Double Field
2040 sage: ComplexHermitianEJA(2, field=RR)
2041 Euclidean Jordan algebra of dimension 4 over Real Field with
2042 53 bits of precision
2046 The dimension of this algebra is `n^2`::
2048 sage: set_random_seed()
2049 sage: n_max = ComplexHermitianEJA._max_random_instance_size()
2050 sage: n = ZZ.random_element(1, n_max)
2051 sage: J = ComplexHermitianEJA(n)
2052 sage: J.dimension() == n^2
2055 The Jordan multiplication is what we think it is::
2057 sage: set_random_seed()
2058 sage: J = ComplexHermitianEJA.random_instance()
2059 sage: x,y = J.random_elements(2)
2060 sage: actual = (x*y).to_matrix()
2061 sage: X = x.to_matrix()
2062 sage: Y = y.to_matrix()
2063 sage: expected = (X*Y + Y*X)/2
2064 sage: actual == expected
2066 sage: J(expected) == x*y
2069 We can change the generator prefix::
2071 sage: ComplexHermitianEJA(2, prefix='z').gens()
2074 We can construct the (trivial) algebra of rank zero::
2076 sage: ComplexHermitianEJA(0)
2077 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2082 def _denormalized_basis(cls
, n
):
2084 Returns a basis for the space of complex Hermitian n-by-n matrices.
2086 Why do we embed these? Basically, because all of numerical linear
2087 algebra assumes that you're working with vectors consisting of `n`
2088 entries from a field and scalars from the same field. There's no way
2089 to tell SageMath that (for example) the vectors contain complex
2090 numbers, while the scalar field is real.
2094 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
2098 sage: set_random_seed()
2099 sage: n = ZZ.random_element(1,5)
2100 sage: B = ComplexHermitianEJA._denormalized_basis(n)
2101 sage: all( M.is_symmetric() for M in B)
2106 R
= PolynomialRing(field
, 'z')
2108 F
= field
.extension(z
**2 + 1, 'I')
2111 # This is like the symmetric case, but we need to be careful:
2113 # * We want conjugate-symmetry, not just symmetry.
2114 # * The diagonal will (as a result) be real.
2117 Eij
= matrix
.zero(F
,n
)
2119 for j
in range(i
+1):
2123 Sij
= cls
.real_embed(Eij
)
2126 # The second one has a minus because it's conjugated.
2127 Eij
[j
,i
] = 1 # Eij = Eij + Eij.transpose()
2128 Sij_real
= cls
.real_embed(Eij
)
2130 # Eij = I*Eij - I*Eij.transpose()
2133 Sij_imag
= cls
.real_embed(Eij
)
2139 # Since we embedded these, we can drop back to the "field" that we
2140 # started with instead of the complex extension "F".
2141 return tuple( s
.change_ring(field
) for s
in S
)
2144 def __init__(self
, n
, **kwargs
):
2145 # We know this is a valid EJA, but will double-check
2146 # if the user passes check_axioms=True.
2147 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2149 super(ComplexHermitianEJA
, self
).__init
__(self
._denormalized
_basis
(n
),
2150 self
.jordan_product
,
2151 self
.trace_inner_product
,
2153 # TODO: this could be factored out somehow, but is left here
2154 # because the MatrixEJA is not presently a subclass of the
2155 # FDEJA class that defines rank() and one().
2156 self
.rank
.set_cache(n
)
2157 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
2158 self
.one
.set_cache(self(idV
))
2161 def _max_random_instance_size():
2162 return 3 # Dimension 9
2165 def random_instance(cls
, **kwargs
):
2167 Return a random instance of this type of algebra.
2169 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2170 return cls(n
, **kwargs
)
2172 class QuaternionMatrixEJA(MatrixEJA
):
2174 # A manual dictionary-cache for the quaternion_extension() method,
2175 # since apparently @classmethods can't also be @cached_methods.
2176 _quaternion_extension
= {}
2179 def quaternion_extension(cls
,field
):
2181 The quaternion field that we embed/unembed, as an extension
2182 of the given ``field``.
2184 if field
in cls
._quaternion
_extension
:
2185 return cls
._quaternion
_extension
[field
]
2187 Q
= QuaternionAlgebra(field
,-1,-1)
2189 cls
._quaternion
_extension
[field
] = Q
2193 def dimension_over_reals():
2197 def real_embed(cls
,M
):
2199 Embed the n-by-n quaternion matrix ``M`` into the space of real
2200 matrices of size 4n-by-4n by first sending each quaternion entry `z
2201 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
2202 c+di],[-c + di, a-bi]]`, and then embedding those into a real
2207 sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA
2211 sage: Q = QuaternionAlgebra(QQ,-1,-1)
2212 sage: i,j,k = Q.gens()
2213 sage: x = 1 + 2*i + 3*j + 4*k
2214 sage: M = matrix(Q, 1, [[x]])
2215 sage: QuaternionMatrixEJA.real_embed(M)
2221 Embedding is a homomorphism (isomorphism, in fact)::
2223 sage: set_random_seed()
2224 sage: n = ZZ.random_element(2)
2225 sage: Q = QuaternionAlgebra(QQ,-1,-1)
2226 sage: X = random_matrix(Q, n)
2227 sage: Y = random_matrix(Q, n)
2228 sage: Xe = QuaternionMatrixEJA.real_embed(X)
2229 sage: Ye = QuaternionMatrixEJA.real_embed(Y)
2230 sage: XYe = QuaternionMatrixEJA.real_embed(X*Y)
2235 super(QuaternionMatrixEJA
,cls
).real_embed(M
)
2236 quaternions
= M
.base_ring()
2239 F
= QuadraticField(-1, 'I')
2244 t
= z
.coefficient_tuple()
2249 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
2250 [-c
+ d
*i
, a
- b
*i
]])
2251 realM
= ComplexMatrixEJA
.real_embed(cplxM
)
2252 blocks
.append(realM
)
2254 # We should have real entries by now, so use the realest field
2255 # we've got for the return value.
2256 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
2261 def real_unembed(cls
,M
):
2263 The inverse of _embed_quaternion_matrix().
2267 sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA
2271 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
2272 ....: [-2, 1, -4, 3],
2273 ....: [-3, 4, 1, -2],
2274 ....: [-4, -3, 2, 1]])
2275 sage: QuaternionMatrixEJA.real_unembed(M)
2276 [1 + 2*i + 3*j + 4*k]
2280 Unembedding is the inverse of embedding::
2282 sage: set_random_seed()
2283 sage: Q = QuaternionAlgebra(QQ, -1, -1)
2284 sage: M = random_matrix(Q, 3)
2285 sage: Me = QuaternionMatrixEJA.real_embed(M)
2286 sage: QuaternionMatrixEJA.real_unembed(Me) == M
2290 super(QuaternionMatrixEJA
,cls
).real_unembed(M
)
2292 d
= cls
.dimension_over_reals()
2294 # Use the base ring of the matrix to ensure that its entries can be
2295 # multiplied by elements of the quaternion algebra.
2296 Q
= cls
.quaternion_extension(M
.base_ring())
2299 # Go top-left to bottom-right (reading order), converting every
2300 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
2303 for l
in range(n
/d
):
2304 for m
in range(n
/d
):
2305 submat
= ComplexMatrixEJA
.real_unembed(
2306 M
[d
*l
:d
*l
+d
,d
*m
:d
*m
+d
] )
2307 if submat
[0,0] != submat
[1,1].conjugate():
2308 raise ValueError('bad on-diagonal submatrix')
2309 if submat
[0,1] != -submat
[1,0].conjugate():
2310 raise ValueError('bad off-diagonal submatrix')
2311 z
= submat
[0,0].real()
2312 z
+= submat
[0,0].imag()*i
2313 z
+= submat
[0,1].real()*j
2314 z
+= submat
[0,1].imag()*k
2317 return matrix(Q
, n
/d
, elements
)
2320 class QuaternionHermitianEJA(ConcreteEJA
, QuaternionMatrixEJA
):
2322 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
2323 matrices, the usual symmetric Jordan product, and the
2324 real-part-of-trace inner product. It has dimension `2n^2 - n` over
2329 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2333 In theory, our "field" can be any subfield of the reals::
2335 sage: QuaternionHermitianEJA(2, field=RDF)
2336 Euclidean Jordan algebra of dimension 6 over Real Double Field
2337 sage: QuaternionHermitianEJA(2, field=RR)
2338 Euclidean Jordan algebra of dimension 6 over Real Field with
2339 53 bits of precision
2343 The dimension of this algebra is `2*n^2 - n`::
2345 sage: set_random_seed()
2346 sage: n_max = QuaternionHermitianEJA._max_random_instance_size()
2347 sage: n = ZZ.random_element(1, n_max)
2348 sage: J = QuaternionHermitianEJA(n)
2349 sage: J.dimension() == 2*(n^2) - n
2352 The Jordan multiplication is what we think it is::
2354 sage: set_random_seed()
2355 sage: J = QuaternionHermitianEJA.random_instance()
2356 sage: x,y = J.random_elements(2)
2357 sage: actual = (x*y).to_matrix()
2358 sage: X = x.to_matrix()
2359 sage: Y = y.to_matrix()
2360 sage: expected = (X*Y + Y*X)/2
2361 sage: actual == expected
2363 sage: J(expected) == x*y
2366 We can change the generator prefix::
2368 sage: QuaternionHermitianEJA(2, prefix='a').gens()
2369 (a0, a1, a2, a3, a4, a5)
2371 We can construct the (trivial) algebra of rank zero::
2373 sage: QuaternionHermitianEJA(0)
2374 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2378 def _denormalized_basis(cls
, n
):
2380 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
2382 Why do we embed these? Basically, because all of numerical
2383 linear algebra assumes that you're working with vectors consisting
2384 of `n` entries from a field and scalars from the same field. There's
2385 no way to tell SageMath that (for example) the vectors contain
2386 complex numbers, while the scalar field is real.
2390 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2394 sage: set_random_seed()
2395 sage: n = ZZ.random_element(1,5)
2396 sage: B = QuaternionHermitianEJA._denormalized_basis(n)
2397 sage: all( M.is_symmetric() for M in B )
2402 Q
= QuaternionAlgebra(QQ
,-1,-1)
2405 # This is like the symmetric case, but we need to be careful:
2407 # * We want conjugate-symmetry, not just symmetry.
2408 # * The diagonal will (as a result) be real.
2411 Eij
= matrix
.zero(Q
,n
)
2413 for j
in range(i
+1):
2417 Sij
= cls
.real_embed(Eij
)
2420 # The second, third, and fourth ones have a minus
2421 # because they're conjugated.
2422 # Eij = Eij + Eij.transpose()
2424 Sij_real
= cls
.real_embed(Eij
)
2426 # Eij = I*(Eij - Eij.transpose())
2429 Sij_I
= cls
.real_embed(Eij
)
2431 # Eij = J*(Eij - Eij.transpose())
2434 Sij_J
= cls
.real_embed(Eij
)
2436 # Eij = K*(Eij - Eij.transpose())
2439 Sij_K
= cls
.real_embed(Eij
)
2445 # Since we embedded these, we can drop back to the "field" that we
2446 # started with instead of the quaternion algebra "Q".
2447 return tuple( s
.change_ring(field
) for s
in S
)
2450 def __init__(self
, n
, **kwargs
):
2451 # We know this is a valid EJA, but will double-check
2452 # if the user passes check_axioms=True.
2453 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2455 super(QuaternionHermitianEJA
, self
).__init
__(self
._denormalized
_basis
(n
),
2456 self
.jordan_product
,
2457 self
.trace_inner_product
,
2459 # TODO: this could be factored out somehow, but is left here
2460 # because the MatrixEJA is not presently a subclass of the
2461 # FDEJA class that defines rank() and one().
2462 self
.rank
.set_cache(n
)
2463 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
2464 self
.one
.set_cache(self(idV
))
2468 def _max_random_instance_size():
2470 The maximum rank of a random QuaternionHermitianEJA.
2472 return 2 # Dimension 6
2475 def random_instance(cls
, **kwargs
):
2477 Return a random instance of this type of algebra.
2479 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2480 return cls(n
, **kwargs
)
2483 class HadamardEJA(ConcreteEJA
):
2485 Return the Euclidean Jordan Algebra corresponding to the set
2486 `R^n` under the Hadamard product.
2488 Note: this is nothing more than the Cartesian product of ``n``
2489 copies of the spin algebra. Once Cartesian product algebras
2490 are implemented, this can go.
2494 sage: from mjo.eja.eja_algebra import HadamardEJA
2498 This multiplication table can be verified by hand::
2500 sage: J = HadamardEJA(3)
2501 sage: e0,e1,e2 = J.gens()
2517 We can change the generator prefix::
2519 sage: HadamardEJA(3, prefix='r').gens()
2523 def __init__(self
, n
, **kwargs
):
2525 jordan_product
= lambda x
,y
: x
2526 inner_product
= lambda x
,y
: x
2528 def jordan_product(x
,y
):
2530 return P( xi
*yi
for (xi
,yi
) in zip(x
,y
) )
2532 def inner_product(x
,y
):
2535 # New defaults for keyword arguments. Don't orthonormalize
2536 # because our basis is already orthonormal with respect to our
2537 # inner-product. Don't check the axioms, because we know this
2538 # is a valid EJA... but do double-check if the user passes
2539 # check_axioms=True. Note: we DON'T override the "check_field"
2540 # default here, because the user can pass in a field!
2541 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2542 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2544 column_basis
= tuple( b
.column() for b
in FreeModule(ZZ
, n
).basis() )
2545 super().__init
__(column_basis
,
2550 self
.rank
.set_cache(n
)
2553 self
.one
.set_cache( self
.zero() )
2555 self
.one
.set_cache( sum(self
.gens()) )
2558 def _max_random_instance_size():
2560 The maximum dimension of a random HadamardEJA.
2565 def random_instance(cls
, **kwargs
):
2567 Return a random instance of this type of algebra.
2569 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2570 return cls(n
, **kwargs
)
2573 class BilinearFormEJA(ConcreteEJA
):
2575 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2576 with the half-trace inner product and jordan product ``x*y =
2577 (<Bx,y>,y_bar>, x0*y_bar + y0*x_bar)`` where `B = 1 \times B22` is
2578 a symmetric positive-definite "bilinear form" matrix. Its
2579 dimension is the size of `B`, and it has rank two in dimensions
2580 larger than two. It reduces to the ``JordanSpinEJA`` when `B` is
2581 the identity matrix of order ``n``.
2583 We insist that the one-by-one upper-left identity block of `B` be
2584 passed in as well so that we can be passed a matrix of size zero
2585 to construct a trivial algebra.
2589 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
2590 ....: JordanSpinEJA)
2594 When no bilinear form is specified, the identity matrix is used,
2595 and the resulting algebra is the Jordan spin algebra::
2597 sage: B = matrix.identity(AA,3)
2598 sage: J0 = BilinearFormEJA(B)
2599 sage: J1 = JordanSpinEJA(3)
2600 sage: J0.multiplication_table() == J0.multiplication_table()
2603 An error is raised if the matrix `B` does not correspond to a
2604 positive-definite bilinear form::
2606 sage: B = matrix.random(QQ,2,3)
2607 sage: J = BilinearFormEJA(B)
2608 Traceback (most recent call last):
2610 ValueError: bilinear form is not positive-definite
2611 sage: B = matrix.zero(QQ,3)
2612 sage: J = BilinearFormEJA(B)
2613 Traceback (most recent call last):
2615 ValueError: bilinear form is not positive-definite
2619 We can create a zero-dimensional algebra::
2621 sage: B = matrix.identity(AA,0)
2622 sage: J = BilinearFormEJA(B)
2626 We can check the multiplication condition given in the Jordan, von
2627 Neumann, and Wigner paper (and also discussed on my "On the
2628 symmetry..." paper). Note that this relies heavily on the standard
2629 choice of basis, as does anything utilizing the bilinear form
2630 matrix. We opt not to orthonormalize the basis, because if we
2631 did, we would have to normalize the `s_{i}` in a similar manner::
2633 sage: set_random_seed()
2634 sage: n = ZZ.random_element(5)
2635 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2636 sage: B11 = matrix.identity(QQ,1)
2637 sage: B22 = M.transpose()*M
2638 sage: B = block_matrix(2,2,[ [B11,0 ],
2640 sage: J = BilinearFormEJA(B, orthonormalize=False)
2641 sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
2642 sage: V = J.vector_space()
2643 sage: sis = [ J( V([0] + (M.inverse()*ei).list()).column() )
2644 ....: for ei in eis ]
2645 sage: actual = [ sis[i]*sis[j]
2646 ....: for i in range(n-1)
2647 ....: for j in range(n-1) ]
2648 sage: expected = [ J.one() if i == j else J.zero()
2649 ....: for i in range(n-1)
2650 ....: for j in range(n-1) ]
2651 sage: actual == expected
2655 def __init__(self
, B
, **kwargs
):
2656 # The matrix "B" is supplied by the user in most cases,
2657 # so it makes sense to check whether or not its positive-
2658 # definite unless we are specifically asked not to...
2659 if ("check_axioms" not in kwargs
) or kwargs
["check_axioms"]:
2660 if not B
.is_positive_definite():
2661 raise ValueError("bilinear form is not positive-definite")
2663 # However, all of the other data for this EJA is computed
2664 # by us in manner that guarantees the axioms are
2665 # satisfied. So, again, unless we are specifically asked to
2666 # verify things, we'll skip the rest of the checks.
2667 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2669 def inner_product(x
,y
):
2670 return (y
.T
*B
*x
)[0,0]
2672 def jordan_product(x
,y
):
2678 z0
= inner_product(y
,x
)
2679 zbar
= y0
*xbar
+ x0
*ybar
2680 return P([z0
] + zbar
.list())
2683 column_basis
= tuple( b
.column() for b
in FreeModule(ZZ
, n
).basis() )
2684 super(BilinearFormEJA
, self
).__init
__(column_basis
,
2689 # The rank of this algebra is two, unless we're in a
2690 # one-dimensional ambient space (because the rank is bounded
2691 # by the ambient dimension).
2692 self
.rank
.set_cache(min(n
,2))
2695 self
.one
.set_cache( self
.zero() )
2697 self
.one
.set_cache( self
.monomial(0) )
2700 def _max_random_instance_size():
2702 The maximum dimension of a random BilinearFormEJA.
2707 def random_instance(cls
, **kwargs
):
2709 Return a random instance of this algebra.
2711 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2713 B
= matrix
.identity(ZZ
, n
)
2714 return cls(B
, **kwargs
)
2716 B11
= matrix
.identity(ZZ
, 1)
2717 M
= matrix
.random(ZZ
, n
-1)
2718 I
= matrix
.identity(ZZ
, n
-1)
2720 while alpha
.is_zero():
2721 alpha
= ZZ
.random_element().abs()
2722 B22
= M
.transpose()*M
+ alpha
*I
2724 from sage
.matrix
.special
import block_matrix
2725 B
= block_matrix(2,2, [ [B11
, ZZ(0) ],
2728 return cls(B
, **kwargs
)
2731 class JordanSpinEJA(BilinearFormEJA
):
2733 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2734 with the usual inner product and jordan product ``x*y =
2735 (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2740 sage: from mjo.eja.eja_algebra import JordanSpinEJA
2744 This multiplication table can be verified by hand::
2746 sage: J = JordanSpinEJA(4)
2747 sage: e0,e1,e2,e3 = J.gens()
2763 We can change the generator prefix::
2765 sage: JordanSpinEJA(2, prefix='B').gens()
2770 Ensure that we have the usual inner product on `R^n`::
2772 sage: set_random_seed()
2773 sage: J = JordanSpinEJA.random_instance()
2774 sage: x,y = J.random_elements(2)
2775 sage: actual = x.inner_product(y)
2776 sage: expected = x.to_vector().inner_product(y.to_vector())
2777 sage: actual == expected
2781 def __init__(self
, n
, **kwargs
):
2782 # This is a special case of the BilinearFormEJA with the
2783 # identity matrix as its bilinear form.
2784 B
= matrix
.identity(ZZ
, n
)
2786 # Don't orthonormalize because our basis is already
2787 # orthonormal with respect to our inner-product.
2788 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2790 # But also don't pass check_field=False here, because the user
2791 # can pass in a field!
2792 super(JordanSpinEJA
, self
).__init
__(B
, **kwargs
)
2795 def _max_random_instance_size():
2797 The maximum dimension of a random JordanSpinEJA.
2802 def random_instance(cls
, **kwargs
):
2804 Return a random instance of this type of algebra.
2806 Needed here to override the implementation for ``BilinearFormEJA``.
2808 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2809 return cls(n
, **kwargs
)
2812 class TrivialEJA(ConcreteEJA
):
2814 The trivial Euclidean Jordan algebra consisting of only a zero element.
2818 sage: from mjo.eja.eja_algebra import TrivialEJA
2822 sage: J = TrivialEJA()
2829 sage: 7*J.one()*12*J.one()
2831 sage: J.one().inner_product(J.one())
2833 sage: J.one().norm()
2835 sage: J.one().subalgebra_generated_by()
2836 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2841 def __init__(self
, **kwargs
):
2842 jordan_product
= lambda x
,y
: x
2843 inner_product
= lambda x
,y
: 0
2846 # New defaults for keyword arguments
2847 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2848 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2850 super(TrivialEJA
, self
).__init
__(basis
,
2854 # The rank is zero using my definition, namely the dimension of the
2855 # largest subalgebra generated by any element.
2856 self
.rank
.set_cache(0)
2857 self
.one
.set_cache( self
.zero() )
2860 def random_instance(cls
, **kwargs
):
2861 # We don't take a "size" argument so the superclass method is
2862 # inappropriate for us.
2863 return cls(**kwargs
)
2866 class CartesianProductEJA(CombinatorialFreeModule_CartesianProduct
,
2867 FiniteDimensionalEJA
):
2869 The external (orthogonal) direct sum of two or more Euclidean
2870 Jordan algebras. Every Euclidean Jordan algebra decomposes into an
2871 orthogonal direct sum of simple Euclidean Jordan algebras which is
2872 then isometric to a Cartesian product, so no generality is lost by
2873 providing only this construction.
2877 sage: from mjo.eja.eja_algebra import (random_eja,
2878 ....: CartesianProductEJA,
2880 ....: JordanSpinEJA,
2881 ....: RealSymmetricEJA)
2885 The Jordan product is inherited from our factors and implemented by
2886 our CombinatorialFreeModule Cartesian product superclass::
2888 sage: set_random_seed()
2889 sage: J1 = HadamardEJA(2)
2890 sage: J2 = RealSymmetricEJA(2)
2891 sage: J = cartesian_product([J1,J2])
2892 sage: x,y = J.random_elements(2)
2896 The ability to retrieve the original factors is implemented by our
2897 CombinatorialFreeModule Cartesian product superclass::
2899 sage: J1 = HadamardEJA(2, field=QQ)
2900 sage: J2 = JordanSpinEJA(3, field=QQ)
2901 sage: J = cartesian_product([J1,J2])
2902 sage: J.cartesian_factors()
2903 (Euclidean Jordan algebra of dimension 2 over Rational Field,
2904 Euclidean Jordan algebra of dimension 3 over Rational Field)
2906 You can provide more than two factors::
2908 sage: J1 = HadamardEJA(2)
2909 sage: J2 = JordanSpinEJA(3)
2910 sage: J3 = RealSymmetricEJA(3)
2911 sage: cartesian_product([J1,J2,J3])
2912 Euclidean Jordan algebra of dimension 2 over Algebraic Real
2913 Field (+) Euclidean Jordan algebra of dimension 3 over Algebraic
2914 Real Field (+) Euclidean Jordan algebra of dimension 6 over
2915 Algebraic Real Field
2917 Rank is additive on a Cartesian product::
2919 sage: J1 = HadamardEJA(1)
2920 sage: J2 = RealSymmetricEJA(2)
2921 sage: J = cartesian_product([J1,J2])
2922 sage: J1.rank.clear_cache()
2923 sage: J2.rank.clear_cache()
2924 sage: J.rank.clear_cache()
2927 sage: J.rank() == J1.rank() + J2.rank()
2930 The same rank computation works over the rationals, with whatever
2933 sage: J1 = HadamardEJA(1, field=QQ, orthonormalize=False)
2934 sage: J2 = RealSymmetricEJA(2, field=QQ, orthonormalize=False)
2935 sage: J = cartesian_product([J1,J2])
2936 sage: J1.rank.clear_cache()
2937 sage: J2.rank.clear_cache()
2938 sage: J.rank.clear_cache()
2941 sage: J.rank() == J1.rank() + J2.rank()
2944 The product algebra will be associative if and only if all of its
2945 components are associative::
2947 sage: J1 = HadamardEJA(2)
2948 sage: J1.is_associative()
2950 sage: J2 = HadamardEJA(3)
2951 sage: J2.is_associative()
2953 sage: J3 = RealSymmetricEJA(3)
2954 sage: J3.is_associative()
2956 sage: CP1 = cartesian_product([J1,J2])
2957 sage: CP1.is_associative()
2959 sage: CP2 = cartesian_product([J1,J3])
2960 sage: CP2.is_associative()
2965 All factors must share the same base field::
2967 sage: J1 = HadamardEJA(2, field=QQ)
2968 sage: J2 = RealSymmetricEJA(2)
2969 sage: CartesianProductEJA((J1,J2))
2970 Traceback (most recent call last):
2972 ValueError: all factors must share the same base field
2974 The cached unit element is the same one that would be computed::
2976 sage: set_random_seed() # long time
2977 sage: J1 = random_eja() # long time
2978 sage: J2 = random_eja() # long time
2979 sage: J = cartesian_product([J1,J2]) # long time
2980 sage: actual = J.one() # long time
2981 sage: J.one.clear_cache() # long time
2982 sage: expected = J.one() # long time
2983 sage: actual == expected # long time
2987 Element
= FiniteDimensionalEJAElement
2990 def __init__(self
, algebras
, **kwargs
):
2991 CombinatorialFreeModule_CartesianProduct
.__init
__(self
,
2994 field
= algebras
[0].base_ring()
2995 if not all( J
.base_ring() == field
for J
in algebras
):
2996 raise ValueError("all factors must share the same base field")
2998 associative
= all( m
.is_associative() for m
in algebras
)
3000 # The definition of matrix_space() and self.basis() relies
3001 # only on the stuff in the CFM_CartesianProduct class, which
3002 # we've already initialized.
3003 Js
= self
.cartesian_factors()
3005 MS
= self
.matrix_space()
3007 MS(tuple( self
.cartesian_projection(i
)(b
).to_matrix()
3008 for i
in range(m
) ))
3009 for b
in self
.basis()
3012 # Define jordan/inner products that operate on that matrix_basis.
3013 def jordan_product(x
,y
):
3015 (Js
[i
](x
[i
])*Js
[i
](y
[i
])).to_matrix() for i
in range(m
)
3018 def inner_product(x
, y
):
3020 Js
[i
](x
[i
]).inner_product(Js
[i
](y
[i
])) for i
in range(m
)
3023 # There's no need to check the field since it already came
3024 # from an EJA. Likewise the axioms are guaranteed to be
3025 # satisfied, unless the guy writing this class sucks.
3027 # If you want the basis to be orthonormalized, orthonormalize
3029 FiniteDimensionalEJA
.__init
__(self
,
3034 orthonormalize
=False,
3035 associative
=associative
,
3036 cartesian_product
=True,
3040 ones
= tuple(J
.one() for J
in algebras
)
3041 self
.one
.set_cache(self
._cartesian
_product
_of
_elements
(ones
))
3042 self
.rank
.set_cache(sum(J
.rank() for J
in algebras
))
3044 def matrix_space(self
):
3046 Return the space that our matrix basis lives in as a Cartesian
3051 sage: from mjo.eja.eja_algebra import (HadamardEJA,
3052 ....: RealSymmetricEJA)
3056 sage: J1 = HadamardEJA(1)
3057 sage: J2 = RealSymmetricEJA(2)
3058 sage: J = cartesian_product([J1,J2])
3059 sage: J.matrix_space()
3060 The Cartesian product of (Full MatrixSpace of 1 by 1 dense
3061 matrices over Algebraic Real Field, Full MatrixSpace of 2
3062 by 2 dense matrices over Algebraic Real Field)
3065 from sage
.categories
.cartesian_product
import cartesian_product
3066 return cartesian_product( [J
.matrix_space()
3067 for J
in self
.cartesian_factors()] )
3070 def cartesian_projection(self
, i
):
3074 sage: from mjo.eja.eja_algebra import (random_eja,
3075 ....: JordanSpinEJA,
3077 ....: RealSymmetricEJA,
3078 ....: ComplexHermitianEJA)
3082 The projection morphisms are Euclidean Jordan algebra
3085 sage: J1 = HadamardEJA(2)
3086 sage: J2 = RealSymmetricEJA(2)
3087 sage: J = cartesian_product([J1,J2])
3088 sage: J.cartesian_projection(0)
3089 Linear operator between finite-dimensional Euclidean Jordan
3090 algebras represented by the matrix:
3093 Domain: Euclidean Jordan algebra of dimension 2 over Algebraic
3094 Real Field (+) Euclidean Jordan algebra of dimension 3 over
3095 Algebraic Real Field
3096 Codomain: Euclidean Jordan algebra of dimension 2 over Algebraic
3098 sage: J.cartesian_projection(1)
3099 Linear operator between finite-dimensional Euclidean Jordan
3100 algebras represented by the matrix:
3104 Domain: Euclidean Jordan algebra of dimension 2 over Algebraic
3105 Real Field (+) Euclidean Jordan algebra of dimension 3 over
3106 Algebraic Real Field
3107 Codomain: Euclidean Jordan algebra of dimension 3 over Algebraic
3110 The projections work the way you'd expect on the vector
3111 representation of an element::
3113 sage: J1 = JordanSpinEJA(2)
3114 sage: J2 = ComplexHermitianEJA(2)
3115 sage: J = cartesian_product([J1,J2])
3116 sage: pi_left = J.cartesian_projection(0)
3117 sage: pi_right = J.cartesian_projection(1)
3118 sage: pi_left(J.one()).to_vector()
3120 sage: pi_right(J.one()).to_vector()
3122 sage: J.one().to_vector()
3127 The answer never changes::
3129 sage: set_random_seed()
3130 sage: J1 = random_eja()
3131 sage: J2 = random_eja()
3132 sage: J = cartesian_product([J1,J2])
3133 sage: P0 = J.cartesian_projection(0)
3134 sage: P1 = J.cartesian_projection(0)
3139 Ji
= self
.cartesian_factors()[i
]
3140 # Requires the fix on Trac 31421/31422 to work!
3141 Pi
= super().cartesian_projection(i
)
3142 return FiniteDimensionalEJAOperator(self
,Ji
,Pi
.matrix())
3145 def cartesian_embedding(self
, i
):
3149 sage: from mjo.eja.eja_algebra import (random_eja,
3150 ....: JordanSpinEJA,
3152 ....: RealSymmetricEJA)
3156 The embedding morphisms are Euclidean Jordan algebra
3159 sage: J1 = HadamardEJA(2)
3160 sage: J2 = RealSymmetricEJA(2)
3161 sage: J = cartesian_product([J1,J2])
3162 sage: J.cartesian_embedding(0)
3163 Linear operator between finite-dimensional Euclidean Jordan
3164 algebras represented by the matrix:
3170 Domain: Euclidean Jordan algebra of dimension 2 over
3171 Algebraic Real Field
3172 Codomain: Euclidean Jordan algebra of dimension 2 over
3173 Algebraic Real Field (+) Euclidean Jordan algebra of
3174 dimension 3 over Algebraic Real Field
3175 sage: J.cartesian_embedding(1)
3176 Linear operator between finite-dimensional Euclidean Jordan
3177 algebras represented by the matrix:
3183 Domain: Euclidean Jordan algebra of dimension 3 over
3184 Algebraic Real Field
3185 Codomain: Euclidean Jordan algebra of dimension 2 over
3186 Algebraic Real Field (+) Euclidean Jordan algebra of
3187 dimension 3 over Algebraic Real Field
3189 The embeddings work the way you'd expect on the vector
3190 representation of an element::
3192 sage: J1 = JordanSpinEJA(3)
3193 sage: J2 = RealSymmetricEJA(2)
3194 sage: J = cartesian_product([J1,J2])
3195 sage: iota_left = J.cartesian_embedding(0)
3196 sage: iota_right = J.cartesian_embedding(1)
3197 sage: iota_left(J1.zero()) == J.zero()
3199 sage: iota_right(J2.zero()) == J.zero()
3201 sage: J1.one().to_vector()
3203 sage: iota_left(J1.one()).to_vector()
3205 sage: J2.one().to_vector()
3207 sage: iota_right(J2.one()).to_vector()
3209 sage: J.one().to_vector()
3214 The answer never changes::
3216 sage: set_random_seed()
3217 sage: J1 = random_eja()
3218 sage: J2 = random_eja()
3219 sage: J = cartesian_product([J1,J2])
3220 sage: E0 = J.cartesian_embedding(0)
3221 sage: E1 = J.cartesian_embedding(0)
3225 Composing a projection with the corresponding inclusion should
3226 produce the identity map, and mismatching them should produce
3229 sage: set_random_seed()
3230 sage: J1 = random_eja()
3231 sage: J2 = random_eja()
3232 sage: J = cartesian_product([J1,J2])
3233 sage: iota_left = J.cartesian_embedding(0)
3234 sage: iota_right = J.cartesian_embedding(1)
3235 sage: pi_left = J.cartesian_projection(0)
3236 sage: pi_right = J.cartesian_projection(1)
3237 sage: pi_left*iota_left == J1.one().operator()
3239 sage: pi_right*iota_right == J2.one().operator()
3241 sage: (pi_left*iota_right).is_zero()
3243 sage: (pi_right*iota_left).is_zero()
3247 Ji
= self
.cartesian_factors()[i
]
3248 # Requires the fix on Trac 31421/31422 to work!
3249 Ei
= super().cartesian_embedding(i
)
3250 return FiniteDimensionalEJAOperator(Ji
,self
,Ei
.matrix())
3254 FiniteDimensionalEJA
.CartesianProduct
= CartesianProductEJA
3256 class RationalBasisCartesianProductEJA(CartesianProductEJA
,
3259 A separate class for products of algebras for which we know a
3264 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
3265 ....: RealSymmetricEJA)
3269 This gives us fast characteristic polynomial computations in
3270 product algebras, too::
3273 sage: J1 = JordanSpinEJA(2)
3274 sage: J2 = RealSymmetricEJA(3)
3275 sage: J = cartesian_product([J1,J2])
3276 sage: J.characteristic_polynomial_of().degree()
3282 def __init__(self
, algebras
, **kwargs
):
3283 CartesianProductEJA
.__init
__(self
, algebras
, **kwargs
)
3285 self
._rational
_algebra
= None
3286 if self
.vector_space().base_field() is not QQ
:
3287 self
._rational
_algebra
= cartesian_product([
3288 r
._rational
_algebra
for r
in algebras
3292 RationalBasisEJA
.CartesianProduct
= RationalBasisCartesianProductEJA
3294 random_eja
= ConcreteEJA
.random_instance