2 Euclidean Jordan Algebras. These are formally-real Jordan Algebras;
3 specifically those where u^2 + v^2 = 0 implies that u = v = 0. They
4 are used in optimization, and have some additional nice methods beyond
5 what can be supported in a general Jordan Algebra.
10 sage: from mjo.eja.eja_algebra import random_eja
15 Euclidean Jordan algebra of dimension...
19 from itertools
import repeat
21 from sage
.algebras
.quatalg
.quaternion_algebra
import QuaternionAlgebra
22 from sage
.categories
.magmatic_algebras
import MagmaticAlgebras
23 from sage
.categories
.sets_cat
import cartesian_product
24 from sage
.combinat
.free_module
import (CombinatorialFreeModule
,
25 CombinatorialFreeModule_CartesianProduct
)
26 from sage
.matrix
.constructor
import matrix
27 from sage
.matrix
.matrix_space
import MatrixSpace
28 from sage
.misc
.cachefunc
import cached_method
29 from sage
.misc
.table
import table
30 from sage
.modules
.free_module
import FreeModule
, VectorSpace
31 from sage
.rings
.all
import (ZZ
, QQ
, AA
, QQbar
, RR
, RLF
, CLF
,
34 from mjo
.eja
.eja_element
import FiniteDimensionalEJAElement
35 from mjo
.eja
.eja_operator
import FiniteDimensionalEJAOperator
36 from mjo
.eja
.eja_utils
import _all2list
, _mat2vec
38 class FiniteDimensionalEJA(CombinatorialFreeModule
):
40 A finite-dimensional Euclidean Jordan algebra.
44 - basis -- a tuple of basis elements in "matrix form," which
45 must be the same form as the arguments to ``jordan_product``
46 and ``inner_product``. In reality, "matrix form" can be either
47 vectors, matrices, or a Cartesian product (ordered tuple)
48 of vectors or matrices. All of these would ideally be vector
49 spaces in sage with no special-casing needed; but in reality
50 we turn vectors into column-matrices and Cartesian products
51 `(a,b)` into column matrices `(a,b)^{T}` after converting
52 `a` and `b` themselves.
54 - jordan_product -- function of two elements (in matrix form)
55 that returns their jordan product in this algebra; this will
56 be applied to ``basis`` to compute a multiplication table for
59 - inner_product -- function of two elements (in matrix form) that
60 returns their inner product. This will be applied to ``basis`` to
61 compute an inner-product table (basically a matrix) for this algebra.
64 Element
= FiniteDimensionalEJAElement
73 cartesian_product
=False,
78 # Keep track of whether or not the matrix basis consists of
79 # tuples, since we need special cases for them damned near
80 # everywhere. This is INDEPENDENT of whether or not the
81 # algebra is a cartesian product, since a subalgebra of a
82 # cartesian product will have a basis of tuples, but will not
83 # in general itself be a cartesian product algebra.
84 self
._matrix
_basis
_is
_cartesian
= False
87 if hasattr(basis
[0], 'cartesian_factors'):
88 self
._matrix
_basis
_is
_cartesian
= True
91 if not field
.is_subring(RR
):
92 # Note: this does return true for the real algebraic
93 # field, the rationals, and any quadratic field where
94 # we've specified a real embedding.
95 raise ValueError("scalar field is not real")
97 # If the basis given to us wasn't over the field that it's
98 # supposed to be over, fix that. Or, you know, crash.
99 if not cartesian_product
:
100 # The field for a cartesian product algebra comes from one
101 # of its factors and is the same for all factors, so
102 # there's no need to "reapply" it on product algebras.
103 if self
._matrix
_basis
_is
_cartesian
:
104 # OK since if n == 0, the basis does not consist of tuples.
105 P
= basis
[0].parent()
106 basis
= tuple( P(tuple(b_i
.change_ring(field
) for b_i
in b
))
109 basis
= tuple( b
.change_ring(field
) for b
in basis
)
113 # Check commutativity of the Jordan and inner-products.
114 # This has to be done before we build the multiplication
115 # and inner-product tables/matrices, because we take
116 # advantage of symmetry in the process.
117 if not all( jordan_product(bi
,bj
) == jordan_product(bj
,bi
)
120 raise ValueError("Jordan product is not commutative")
122 if not all( inner_product(bi
,bj
) == inner_product(bj
,bi
)
125 raise ValueError("inner-product is not commutative")
128 category
= MagmaticAlgebras(field
).FiniteDimensional()
129 category
= category
.WithBasis().Unital()
131 # Element subalgebras can take advantage of this.
132 category
= category
.Associative()
133 if cartesian_product
:
134 category
= category
.CartesianProducts()
136 # Call the superclass constructor so that we can use its from_vector()
137 # method to build our multiplication table.
138 CombinatorialFreeModule
.__init
__(self
,
145 # Now comes all of the hard work. We'll be constructing an
146 # ambient vector space V that our (vectorized) basis lives in,
147 # as well as a subspace W of V spanned by those (vectorized)
148 # basis elements. The W-coordinates are the coefficients that
149 # we see in things like x = 1*e1 + 2*e2.
154 degree
= len(_all2list(basis
[0]))
156 # Build an ambient space that fits our matrix basis when
157 # written out as "long vectors."
158 V
= VectorSpace(field
, degree
)
160 # The matrix that will hole the orthonormal -> unorthonormal
161 # coordinate transformation.
162 self
._deortho
_matrix
= None
165 # Save a copy of the un-orthonormalized basis for later.
166 # Convert it to ambient V (vector) coordinates while we're
167 # at it, because we'd have to do it later anyway.
168 deortho_vector_basis
= tuple( V(_all2list(b
)) for b
in basis
)
170 from mjo
.eja
.eja_utils
import gram_schmidt
171 basis
= tuple(gram_schmidt(basis
, inner_product
))
173 # Save the (possibly orthonormalized) matrix basis for
175 self
._matrix
_basis
= basis
177 # Now create the vector space for the algebra, which will have
178 # its own set of non-ambient coordinates (in terms of the
180 vector_basis
= tuple( V(_all2list(b
)) for b
in basis
)
181 W
= V
.span_of_basis( vector_basis
, check
=check_axioms
)
184 # Now "W" is the vector space of our algebra coordinates. The
185 # variables "X1", "X2",... refer to the entries of vectors in
186 # W. Thus to convert back and forth between the orthonormal
187 # coordinates and the given ones, we need to stick the original
189 U
= V
.span_of_basis( deortho_vector_basis
, check
=check_axioms
)
190 self
._deortho
_matrix
= matrix( U
.coordinate_vector(q
)
191 for q
in vector_basis
)
194 # Now we actually compute the multiplication and inner-product
195 # tables/matrices using the possibly-orthonormalized basis.
196 self
._inner
_product
_matrix
= matrix
.identity(field
, n
)
197 self
._multiplication
_table
= [ [0 for j
in range(i
+1)]
200 # Note: the Jordan and inner-products are defined in terms
201 # of the ambient basis. It's important that their arguments
202 # are in ambient coordinates as well.
205 # ortho basis w.r.t. ambient coords
209 # The jordan product returns a matrixy answer, so we
210 # have to convert it to the algebra coordinates.
211 elt
= jordan_product(q_i
, q_j
)
212 elt
= W
.coordinate_vector(V(_all2list(elt
)))
213 self
._multiplication
_table
[i
][j
] = self
.from_vector(elt
)
215 if not orthonormalize
:
216 # If we're orthonormalizing the basis with respect
217 # to an inner-product, then the inner-product
218 # matrix with respect to the resulting basis is
219 # just going to be the identity.
220 ip
= inner_product(q_i
, q_j
)
221 self
._inner
_product
_matrix
[i
,j
] = ip
222 self
._inner
_product
_matrix
[j
,i
] = ip
224 self
._inner
_product
_matrix
._cache
= {'hermitian': True}
225 self
._inner
_product
_matrix
.set_immutable()
228 if not self
._is
_jordanian
():
229 raise ValueError("Jordan identity does not hold")
230 if not self
._inner
_product
_is
_associative
():
231 raise ValueError("inner product is not associative")
234 def _coerce_map_from_base_ring(self
):
236 Disable the map from the base ring into the algebra.
238 Performing a nonsense conversion like this automatically
239 is counterpedagogical. The fallback is to try the usual
240 element constructor, which should also fail.
244 sage: from mjo.eja.eja_algebra import random_eja
248 sage: set_random_seed()
249 sage: J = random_eja()
251 Traceback (most recent call last):
253 ValueError: not an element of this algebra
259 def product_on_basis(self
, i
, j
):
260 # We only stored the lower-triangular portion of the
261 # multiplication table.
263 return self
._multiplication
_table
[i
][j
]
265 return self
._multiplication
_table
[j
][i
]
267 def inner_product(self
, x
, y
):
269 The inner product associated with this Euclidean Jordan algebra.
271 Defaults to the trace inner product, but can be overridden by
272 subclasses if they are sure that the necessary properties are
277 sage: from mjo.eja.eja_algebra import (random_eja,
279 ....: BilinearFormEJA)
283 Our inner product is "associative," which means the following for
284 a symmetric bilinear form::
286 sage: set_random_seed()
287 sage: J = random_eja()
288 sage: x,y,z = J.random_elements(3)
289 sage: (x*y).inner_product(z) == y.inner_product(x*z)
294 Ensure that this is the usual inner product for the algebras
297 sage: set_random_seed()
298 sage: J = HadamardEJA.random_instance()
299 sage: x,y = J.random_elements(2)
300 sage: actual = x.inner_product(y)
301 sage: expected = x.to_vector().inner_product(y.to_vector())
302 sage: actual == expected
305 Ensure that this is one-half of the trace inner-product in a
306 BilinearFormEJA that isn't just the reals (when ``n`` isn't
307 one). This is in Faraut and Koranyi, and also my "On the
310 sage: set_random_seed()
311 sage: J = BilinearFormEJA.random_instance()
312 sage: n = J.dimension()
313 sage: x = J.random_element()
314 sage: y = J.random_element()
315 sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
319 B
= self
._inner
_product
_matrix
320 return (B
*x
.to_vector()).inner_product(y
.to_vector())
323 def is_associative(self
):
325 Return whether or not this algebra's Jordan product is associative.
329 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
333 sage: J = ComplexHermitianEJA(3, field=QQ, orthonormalize=False)
334 sage: J.is_associative()
336 sage: x = sum(J.gens())
337 sage: A = x.subalgebra_generated_by(orthonormalize=False)
338 sage: A.is_associative()
342 return "Associative" in self
.category().axioms()
344 def _is_jordanian(self
):
346 Whether or not this algebra's multiplication table respects the
347 Jordan identity `(x^{2})(xy) = x(x^{2}y)`.
349 We only check one arrangement of `x` and `y`, so for a
350 ``True`` result to be truly true, you should also check
351 :meth:`is_commutative`. This method should of course always
352 return ``True``, unless this algebra was constructed with
353 ``check_axioms=False`` and passed an invalid multiplication table.
355 return all( (self
.gens()[i
]**2)*(self
.gens()[i
]*self
.gens()[j
])
357 (self
.gens()[i
])*((self
.gens()[i
]**2)*self
.gens()[j
])
358 for i
in range(self
.dimension())
359 for j
in range(self
.dimension()) )
361 def _inner_product_is_associative(self
):
363 Return whether or not this algebra's inner product `B` is
364 associative; that is, whether or not `B(xy,z) = B(x,yz)`.
366 This method should of course always return ``True``, unless
367 this algebra was constructed with ``check_axioms=False`` and
368 passed an invalid Jordan or inner-product.
371 # Used to check whether or not something is zero in an inexact
372 # ring. This number is sufficient to allow the construction of
373 # QuaternionHermitianEJA(2, field=RDF) with check_axioms=True.
376 for i
in range(self
.dimension()):
377 for j
in range(self
.dimension()):
378 for k
in range(self
.dimension()):
382 diff
= (x
*y
).inner_product(z
) - x
.inner_product(y
*z
)
384 if self
.base_ring().is_exact():
388 if diff
.abs() > epsilon
:
393 def _element_constructor_(self
, elt
):
395 Construct an element of this algebra from its vector or matrix
398 This gets called only after the parent element _call_ method
399 fails to find a coercion for the argument.
403 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
405 ....: RealSymmetricEJA)
409 The identity in `S^n` is converted to the identity in the EJA::
411 sage: J = RealSymmetricEJA(3)
412 sage: I = matrix.identity(QQ,3)
413 sage: J(I) == J.one()
416 This skew-symmetric matrix can't be represented in the EJA::
418 sage: J = RealSymmetricEJA(3)
419 sage: A = matrix(QQ,3, lambda i,j: i-j)
421 Traceback (most recent call last):
423 ValueError: not an element of this algebra
425 Tuples work as well, provided that the matrix basis for the
426 algebra consists of them::
428 sage: J1 = HadamardEJA(3)
429 sage: J2 = RealSymmetricEJA(2)
430 sage: J = cartesian_product([J1,J2])
431 sage: J( (J1.matrix_basis()[1], J2.matrix_basis()[2]) )
436 Ensure that we can convert any element of the two non-matrix
437 simple algebras (whose matrix representations are columns)
438 back and forth faithfully::
440 sage: set_random_seed()
441 sage: J = HadamardEJA.random_instance()
442 sage: x = J.random_element()
443 sage: J(x.to_vector().column()) == x
445 sage: J = JordanSpinEJA.random_instance()
446 sage: x = J.random_element()
447 sage: J(x.to_vector().column()) == x
451 msg
= "not an element of this algebra"
453 # The superclass implementation of random_element()
454 # needs to be able to coerce "0" into the algebra.
456 elif elt
in self
.base_ring():
457 # Ensure that no base ring -> algebra coercion is performed
458 # by this method. There's some stupidity in sage that would
459 # otherwise propagate to this method; for example, sage thinks
460 # that the integer 3 belongs to the space of 2-by-2 matrices.
461 raise ValueError(msg
)
465 except (AttributeError, TypeError):
466 # Try to convert a vector into a column-matrix
469 if elt
not in self
.matrix_space():
470 raise ValueError(msg
)
472 # Thanks for nothing! Matrix spaces aren't vector spaces in
473 # Sage, so we have to figure out its matrix-basis coordinates
474 # ourselves. We use the basis space's ring instead of the
475 # element's ring because the basis space might be an algebraic
476 # closure whereas the base ring of the 3-by-3 identity matrix
477 # could be QQ instead of QQbar.
479 # And, we also have to handle Cartesian product bases (when
480 # the matric basis consists of tuples) here. The "good news"
481 # is that we're already converting everything to long vectors,
482 # and that strategy works for tuples as well.
484 # We pass check=False because the matrix basis is "guaranteed"
485 # to be linearly independent... right? Ha ha.
487 V
= VectorSpace(self
.base_ring(), len(elt
))
488 W
= V
.span_of_basis( (V(_all2list(s
)) for s
in self
.matrix_basis()),
492 coords
= W
.coordinate_vector(V(elt
))
493 except ArithmeticError: # vector is not in free module
494 raise ValueError(msg
)
496 return self
.from_vector(coords
)
500 Return a string representation of ``self``.
504 sage: from mjo.eja.eja_algebra import JordanSpinEJA
508 Ensure that it says what we think it says::
510 sage: JordanSpinEJA(2, field=AA)
511 Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
512 sage: JordanSpinEJA(3, field=RDF)
513 Euclidean Jordan algebra of dimension 3 over Real Double Field
516 fmt
= "Euclidean Jordan algebra of dimension {} over {}"
517 return fmt
.format(self
.dimension(), self
.base_ring())
521 def characteristic_polynomial_of(self
):
523 Return the algebra's "characteristic polynomial of" function,
524 which is itself a multivariate polynomial that, when evaluated
525 at the coordinates of some algebra element, returns that
526 element's characteristic polynomial.
528 The resulting polynomial has `n+1` variables, where `n` is the
529 dimension of this algebra. The first `n` variables correspond to
530 the coordinates of an algebra element: when evaluated at the
531 coordinates of an algebra element with respect to a certain
532 basis, the result is a univariate polynomial (in the one
533 remaining variable ``t``), namely the characteristic polynomial
538 sage: from mjo.eja.eja_algebra import JordanSpinEJA, TrivialEJA
542 The characteristic polynomial in the spin algebra is given in
543 Alizadeh, Example 11.11::
545 sage: J = JordanSpinEJA(3)
546 sage: p = J.characteristic_polynomial_of(); p
547 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
548 sage: xvec = J.one().to_vector()
552 By definition, the characteristic polynomial is a monic
553 degree-zero polynomial in a rank-zero algebra. Note that
554 Cayley-Hamilton is indeed satisfied since the polynomial
555 ``1`` evaluates to the identity element of the algebra on
558 sage: J = TrivialEJA()
559 sage: J.characteristic_polynomial_of()
566 # The list of coefficient polynomials a_0, a_1, a_2, ..., a_(r-1).
567 a
= self
._charpoly
_coefficients
()
569 # We go to a bit of trouble here to reorder the
570 # indeterminates, so that it's easier to evaluate the
571 # characteristic polynomial at x's coordinates and get back
572 # something in terms of t, which is what we want.
573 S
= PolynomialRing(self
.base_ring(),'t')
577 S
= PolynomialRing(S
, R
.variable_names())
580 return (t
**r
+ sum( a
[k
]*(t
**k
) for k
in range(r
) ))
582 def coordinate_polynomial_ring(self
):
584 The multivariate polynomial ring in which this algebra's
585 :meth:`characteristic_polynomial_of` lives.
589 sage: from mjo.eja.eja_algebra import (HadamardEJA,
590 ....: RealSymmetricEJA)
594 sage: J = HadamardEJA(2)
595 sage: J.coordinate_polynomial_ring()
596 Multivariate Polynomial Ring in X1, X2...
597 sage: J = RealSymmetricEJA(3,field=QQ,orthonormalize=False)
598 sage: J.coordinate_polynomial_ring()
599 Multivariate Polynomial Ring in X1, X2, X3, X4, X5, X6...
602 var_names
= tuple( "X%d" % z
for z
in range(1, self
.dimension()+1) )
603 return PolynomialRing(self
.base_ring(), var_names
)
605 def inner_product(self
, x
, y
):
607 The inner product associated with this Euclidean Jordan algebra.
609 Defaults to the trace inner product, but can be overridden by
610 subclasses if they are sure that the necessary properties are
615 sage: from mjo.eja.eja_algebra import (random_eja,
617 ....: BilinearFormEJA)
621 Our inner product is "associative," which means the following for
622 a symmetric bilinear form::
624 sage: set_random_seed()
625 sage: J = random_eja()
626 sage: x,y,z = J.random_elements(3)
627 sage: (x*y).inner_product(z) == y.inner_product(x*z)
632 Ensure that this is the usual inner product for the algebras
635 sage: set_random_seed()
636 sage: J = HadamardEJA.random_instance()
637 sage: x,y = J.random_elements(2)
638 sage: actual = x.inner_product(y)
639 sage: expected = x.to_vector().inner_product(y.to_vector())
640 sage: actual == expected
643 Ensure that this is one-half of the trace inner-product in a
644 BilinearFormEJA that isn't just the reals (when ``n`` isn't
645 one). This is in Faraut and Koranyi, and also my "On the
648 sage: set_random_seed()
649 sage: J = BilinearFormEJA.random_instance()
650 sage: n = J.dimension()
651 sage: x = J.random_element()
652 sage: y = J.random_element()
653 sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
656 B
= self
._inner
_product
_matrix
657 return (B
*x
.to_vector()).inner_product(y
.to_vector())
660 def is_trivial(self
):
662 Return whether or not this algebra is trivial.
664 A trivial algebra contains only the zero element.
668 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
673 sage: J = ComplexHermitianEJA(3)
679 sage: J = TrivialEJA()
684 return self
.dimension() == 0
687 def multiplication_table(self
):
689 Return a visual representation of this algebra's multiplication
690 table (on basis elements).
694 sage: from mjo.eja.eja_algebra import JordanSpinEJA
698 sage: J = JordanSpinEJA(4)
699 sage: J.multiplication_table()
700 +----++----+----+----+----+
701 | * || e0 | e1 | e2 | e3 |
702 +====++====+====+====+====+
703 | e0 || e0 | e1 | e2 | e3 |
704 +----++----+----+----+----+
705 | e1 || e1 | e0 | 0 | 0 |
706 +----++----+----+----+----+
707 | e2 || e2 | 0 | e0 | 0 |
708 +----++----+----+----+----+
709 | e3 || e3 | 0 | 0 | e0 |
710 +----++----+----+----+----+
714 # Prepend the header row.
715 M
= [["*"] + list(self
.gens())]
717 # And to each subsequent row, prepend an entry that belongs to
718 # the left-side "header column."
719 M
+= [ [self
.gens()[i
]] + [ self
.product_on_basis(i
,j
)
723 return table(M
, header_row
=True, header_column
=True, frame
=True)
726 def matrix_basis(self
):
728 Return an (often more natural) representation of this algebras
729 basis as an ordered tuple of matrices.
731 Every finite-dimensional Euclidean Jordan Algebra is a, up to
732 Jordan isomorphism, a direct sum of five simple
733 algebras---four of which comprise Hermitian matrices. And the
734 last type of algebra can of course be thought of as `n`-by-`1`
735 column matrices (ambiguusly called column vectors) to avoid
736 special cases. As a result, matrices (and column vectors) are
737 a natural representation format for Euclidean Jordan algebra
740 But, when we construct an algebra from a basis of matrices,
741 those matrix representations are lost in favor of coordinate
742 vectors *with respect to* that basis. We could eventually
743 convert back if we tried hard enough, but having the original
744 representations handy is valuable enough that we simply store
745 them and return them from this method.
747 Why implement this for non-matrix algebras? Avoiding special
748 cases for the :class:`BilinearFormEJA` pays with simplicity in
749 its own right. But mainly, we would like to be able to assume
750 that elements of a :class:`CartesianProductEJA` can be displayed
751 nicely, without having to have special classes for direct sums
752 one of whose components was a matrix algebra.
756 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
757 ....: RealSymmetricEJA)
761 sage: J = RealSymmetricEJA(2)
763 Finite family {0: e0, 1: e1, 2: e2}
764 sage: J.matrix_basis()
766 [1 0] [ 0 0.7071067811865475?] [0 0]
767 [0 0], [0.7071067811865475? 0], [0 1]
772 sage: J = JordanSpinEJA(2)
774 Finite family {0: e0, 1: e1}
775 sage: J.matrix_basis()
781 return self
._matrix
_basis
784 def matrix_space(self
):
786 Return the matrix space in which this algebra's elements live, if
787 we think of them as matrices (including column vectors of the
790 Generally this will be an `n`-by-`1` column-vector space,
791 except when the algebra is trivial. There it's `n`-by-`n`
792 (where `n` is zero), to ensure that two elements of the matrix
793 space (empty matrices) can be multiplied.
795 Matrix algebras override this with something more useful.
797 if self
.is_trivial():
798 return MatrixSpace(self
.base_ring(), 0)
800 return self
.matrix_basis()[0].parent()
806 Return the unit element of this algebra.
810 sage: from mjo.eja.eja_algebra import (HadamardEJA,
815 We can compute unit element in the Hadamard EJA::
817 sage: J = HadamardEJA(5)
819 e0 + e1 + e2 + e3 + e4
821 The unit element in the Hadamard EJA is inherited in the
822 subalgebras generated by its elements::
824 sage: J = HadamardEJA(5)
826 e0 + e1 + e2 + e3 + e4
827 sage: x = sum(J.gens())
828 sage: A = x.subalgebra_generated_by(orthonormalize=False)
831 sage: A.one().superalgebra_element()
832 e0 + e1 + e2 + e3 + e4
836 The identity element acts like the identity, regardless of
837 whether or not we orthonormalize::
839 sage: set_random_seed()
840 sage: J = random_eja()
841 sage: x = J.random_element()
842 sage: J.one()*x == x and x*J.one() == x
844 sage: A = x.subalgebra_generated_by()
845 sage: y = A.random_element()
846 sage: A.one()*y == y and y*A.one() == y
851 sage: set_random_seed()
852 sage: J = random_eja(field=QQ, orthonormalize=False)
853 sage: x = J.random_element()
854 sage: J.one()*x == x and x*J.one() == x
856 sage: A = x.subalgebra_generated_by(orthonormalize=False)
857 sage: y = A.random_element()
858 sage: A.one()*y == y and y*A.one() == y
861 The matrix of the unit element's operator is the identity,
862 regardless of the base field and whether or not we
865 sage: set_random_seed()
866 sage: J = random_eja()
867 sage: actual = J.one().operator().matrix()
868 sage: expected = matrix.identity(J.base_ring(), J.dimension())
869 sage: actual == expected
871 sage: x = J.random_element()
872 sage: A = x.subalgebra_generated_by()
873 sage: actual = A.one().operator().matrix()
874 sage: expected = matrix.identity(A.base_ring(), A.dimension())
875 sage: actual == expected
880 sage: set_random_seed()
881 sage: J = random_eja(field=QQ, orthonormalize=False)
882 sage: actual = J.one().operator().matrix()
883 sage: expected = matrix.identity(J.base_ring(), J.dimension())
884 sage: actual == expected
886 sage: x = J.random_element()
887 sage: A = x.subalgebra_generated_by(orthonormalize=False)
888 sage: actual = A.one().operator().matrix()
889 sage: expected = matrix.identity(A.base_ring(), A.dimension())
890 sage: actual == expected
893 Ensure that the cached unit element (often precomputed by
894 hand) agrees with the computed one::
896 sage: set_random_seed()
897 sage: J = random_eja()
898 sage: cached = J.one()
899 sage: J.one.clear_cache()
900 sage: J.one() == cached
905 sage: set_random_seed()
906 sage: J = random_eja(field=QQ, orthonormalize=False)
907 sage: cached = J.one()
908 sage: J.one.clear_cache()
909 sage: J.one() == cached
913 # We can brute-force compute the matrices of the operators
914 # that correspond to the basis elements of this algebra.
915 # If some linear combination of those basis elements is the
916 # algebra identity, then the same linear combination of
917 # their matrices has to be the identity matrix.
919 # Of course, matrices aren't vectors in sage, so we have to
920 # appeal to the "long vectors" isometry.
921 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
923 # Now we use basic linear algebra to find the coefficients,
924 # of the matrices-as-vectors-linear-combination, which should
925 # work for the original algebra basis too.
926 A
= matrix(self
.base_ring(), oper_vecs
)
928 # We used the isometry on the left-hand side already, but we
929 # still need to do it for the right-hand side. Recall that we
930 # wanted something that summed to the identity matrix.
931 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
933 # Now if there's an identity element in the algebra, this
934 # should work. We solve on the left to avoid having to
935 # transpose the matrix "A".
936 return self
.from_vector(A
.solve_left(b
))
939 def peirce_decomposition(self
, c
):
941 The Peirce decomposition of this algebra relative to the
944 In the future, this can be extended to a complete system of
945 orthogonal idempotents.
949 - ``c`` -- an idempotent of this algebra.
953 A triple (J0, J5, J1) containing two subalgebras and one subspace
956 - ``J0`` -- the algebra on the eigenspace of ``c.operator()``
957 corresponding to the eigenvalue zero.
959 - ``J5`` -- the eigenspace (NOT a subalgebra) of ``c.operator()``
960 corresponding to the eigenvalue one-half.
962 - ``J1`` -- the algebra on the eigenspace of ``c.operator()``
963 corresponding to the eigenvalue one.
965 These are the only possible eigenspaces for that operator, and this
966 algebra is a direct sum of them. The spaces ``J0`` and ``J1`` are
967 orthogonal, and are subalgebras of this algebra with the appropriate
972 sage: from mjo.eja.eja_algebra import random_eja, RealSymmetricEJA
976 The canonical example comes from the symmetric matrices, which
977 decompose into diagonal and off-diagonal parts::
979 sage: J = RealSymmetricEJA(3)
980 sage: C = matrix(QQ, [ [1,0,0],
984 sage: J0,J5,J1 = J.peirce_decomposition(c)
986 Euclidean Jordan algebra of dimension 1...
988 Vector space of degree 6 and dimension 2...
990 Euclidean Jordan algebra of dimension 3...
991 sage: J0.one().to_matrix()
995 sage: orig_df = AA.options.display_format
996 sage: AA.options.display_format = 'radical'
997 sage: J.from_vector(J5.basis()[0]).to_matrix()
1001 sage: J.from_vector(J5.basis()[1]).to_matrix()
1005 sage: AA.options.display_format = orig_df
1006 sage: J1.one().to_matrix()
1013 Every algebra decomposes trivially with respect to its identity
1016 sage: set_random_seed()
1017 sage: J = random_eja()
1018 sage: J0,J5,J1 = J.peirce_decomposition(J.one())
1019 sage: J0.dimension() == 0 and J5.dimension() == 0
1021 sage: J1.superalgebra() == J and J1.dimension() == J.dimension()
1024 The decomposition is into eigenspaces, and its components are
1025 therefore necessarily orthogonal. Moreover, the identity
1026 elements in the two subalgebras are the projections onto their
1027 respective subspaces of the superalgebra's identity element::
1029 sage: set_random_seed()
1030 sage: J = random_eja()
1031 sage: x = J.random_element()
1032 sage: if not J.is_trivial():
1033 ....: while x.is_nilpotent():
1034 ....: x = J.random_element()
1035 sage: c = x.subalgebra_idempotent()
1036 sage: J0,J5,J1 = J.peirce_decomposition(c)
1038 sage: for (w,y,z) in zip(J0.basis(), J5.basis(), J1.basis()):
1039 ....: w = w.superalgebra_element()
1040 ....: y = J.from_vector(y)
1041 ....: z = z.superalgebra_element()
1042 ....: ipsum += w.inner_product(y).abs()
1043 ....: ipsum += w.inner_product(z).abs()
1044 ....: ipsum += y.inner_product(z).abs()
1047 sage: J1(c) == J1.one()
1049 sage: J0(J.one() - c) == J0.one()
1053 if not c
.is_idempotent():
1054 raise ValueError("element is not idempotent: %s" % c
)
1056 # Default these to what they should be if they turn out to be
1057 # trivial, because eigenspaces_left() won't return eigenvalues
1058 # corresponding to trivial spaces (e.g. it returns only the
1059 # eigenspace corresponding to lambda=1 if you take the
1060 # decomposition relative to the identity element).
1061 trivial
= self
.subalgebra(())
1062 J0
= trivial
# eigenvalue zero
1063 J5
= VectorSpace(self
.base_ring(), 0) # eigenvalue one-half
1064 J1
= trivial
# eigenvalue one
1066 for (eigval
, eigspace
) in c
.operator().matrix().right_eigenspaces():
1067 if eigval
== ~
(self
.base_ring()(2)):
1070 gens
= tuple( self
.from_vector(b
) for b
in eigspace
.basis() )
1071 subalg
= self
.subalgebra(gens
, check_axioms
=False)
1077 raise ValueError("unexpected eigenvalue: %s" % eigval
)
1082 def random_element(self
, thorough
=False):
1084 Return a random element of this algebra.
1086 Our algebra superclass method only returns a linear
1087 combination of at most two basis elements. We instead
1088 want the vector space "random element" method that
1089 returns a more diverse selection.
1093 - ``thorough`` -- (boolean; default False) whether or not we
1094 should generate irrational coefficients for the random
1095 element when our base ring is irrational; this slows the
1096 algebra operations to a crawl, but any truly random method
1100 # For a general base ring... maybe we can trust this to do the
1101 # right thing? Unlikely, but.
1102 V
= self
.vector_space()
1103 v
= V
.random_element()
1105 if self
.base_ring() is AA
:
1106 # The "random element" method of the algebraic reals is
1107 # stupid at the moment, and only returns integers between
1108 # -2 and 2, inclusive:
1110 # https://trac.sagemath.org/ticket/30875
1112 # Instead, we implement our own "random vector" method,
1113 # and then coerce that into the algebra. We use the vector
1114 # space degree here instead of the dimension because a
1115 # subalgebra could (for example) be spanned by only two
1116 # vectors, each with five coordinates. We need to
1117 # generate all five coordinates.
1119 v
*= QQbar
.random_element().real()
1121 v
*= QQ
.random_element()
1123 return self
.from_vector(V
.coordinate_vector(v
))
1125 def random_elements(self
, count
, thorough
=False):
1127 Return ``count`` random elements as a tuple.
1131 - ``thorough`` -- (boolean; default False) whether or not we
1132 should generate irrational coefficients for the random
1133 elements when our base ring is irrational; this slows the
1134 algebra operations to a crawl, but any truly random method
1139 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1143 sage: J = JordanSpinEJA(3)
1144 sage: x,y,z = J.random_elements(3)
1145 sage: all( [ x in J, y in J, z in J ])
1147 sage: len( J.random_elements(10) ) == 10
1151 return tuple( self
.random_element(thorough
)
1152 for idx
in range(count
) )
1156 def _charpoly_coefficients(self
):
1158 The `r` polynomial coefficients of the "characteristic polynomial
1163 sage: from mjo.eja.eja_algebra import random_eja
1167 The theory shows that these are all homogeneous polynomials of
1170 sage: set_random_seed()
1171 sage: J = random_eja()
1172 sage: all(p.is_homogeneous() for p in J._charpoly_coefficients())
1176 n
= self
.dimension()
1177 R
= self
.coordinate_polynomial_ring()
1179 F
= R
.fraction_field()
1182 # From a result in my book, these are the entries of the
1183 # basis representation of L_x.
1184 return sum( vars[k
]*self
.gens()[k
].operator().matrix()[i
,j
]
1187 L_x
= matrix(F
, n
, n
, L_x_i_j
)
1190 if self
.rank
.is_in_cache():
1192 # There's no need to pad the system with redundant
1193 # columns if we *know* they'll be redundant.
1196 # Compute an extra power in case the rank is equal to
1197 # the dimension (otherwise, we would stop at x^(r-1)).
1198 x_powers
= [ (L_x
**k
)*self
.one().to_vector()
1199 for k
in range(n
+1) ]
1200 A
= matrix
.column(F
, x_powers
[:n
])
1201 AE
= A
.extended_echelon_form()
1208 # The theory says that only the first "r" coefficients are
1209 # nonzero, and they actually live in the original polynomial
1210 # ring and not the fraction field. We negate them because in
1211 # the actual characteristic polynomial, they get moved to the
1212 # other side where x^r lives. We don't bother to trim A_rref
1213 # down to a square matrix and solve the resulting system,
1214 # because the upper-left r-by-r portion of A_rref is
1215 # guaranteed to be the identity matrix, so e.g.
1217 # A_rref.solve_right(Y)
1219 # would just be returning Y.
1220 return (-E
*b
)[:r
].change_ring(R
)
1225 Return the rank of this EJA.
1227 This is a cached method because we know the rank a priori for
1228 all of the algebras we can construct. Thus we can avoid the
1229 expensive ``_charpoly_coefficients()`` call unless we truly
1230 need to compute the whole characteristic polynomial.
1234 sage: from mjo.eja.eja_algebra import (HadamardEJA,
1235 ....: JordanSpinEJA,
1236 ....: RealSymmetricEJA,
1237 ....: ComplexHermitianEJA,
1238 ....: QuaternionHermitianEJA,
1243 The rank of the Jordan spin algebra is always two::
1245 sage: JordanSpinEJA(2).rank()
1247 sage: JordanSpinEJA(3).rank()
1249 sage: JordanSpinEJA(4).rank()
1252 The rank of the `n`-by-`n` Hermitian real, complex, or
1253 quaternion matrices is `n`::
1255 sage: RealSymmetricEJA(4).rank()
1257 sage: ComplexHermitianEJA(3).rank()
1259 sage: QuaternionHermitianEJA(2).rank()
1264 Ensure that every EJA that we know how to construct has a
1265 positive integer rank, unless the algebra is trivial in
1266 which case its rank will be zero::
1268 sage: set_random_seed()
1269 sage: J = random_eja()
1273 sage: r > 0 or (r == 0 and J.is_trivial())
1276 Ensure that computing the rank actually works, since the ranks
1277 of all simple algebras are known and will be cached by default::
1279 sage: set_random_seed() # long time
1280 sage: J = random_eja() # long time
1281 sage: cached = J.rank() # long time
1282 sage: J.rank.clear_cache() # long time
1283 sage: J.rank() == cached # long time
1287 return len(self
._charpoly
_coefficients
())
1290 def subalgebra(self
, basis
, **kwargs
):
1292 Create a subalgebra of this algebra from the given basis.
1294 from mjo
.eja
.eja_subalgebra
import FiniteDimensionalEJASubalgebra
1295 return FiniteDimensionalEJASubalgebra(self
, basis
, **kwargs
)
1298 def vector_space(self
):
1300 Return the vector space that underlies this algebra.
1304 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1308 sage: J = RealSymmetricEJA(2)
1309 sage: J.vector_space()
1310 Vector space of dimension 3 over...
1313 return self
.zero().to_vector().parent().ambient_vector_space()
1317 class RationalBasisEJA(FiniteDimensionalEJA
):
1319 New class for algebras whose supplied basis elements have all rational entries.
1323 sage: from mjo.eja.eja_algebra import BilinearFormEJA
1327 The supplied basis is orthonormalized by default::
1329 sage: B = matrix(QQ, [[1, 0, 0], [0, 25, -32], [0, -32, 41]])
1330 sage: J = BilinearFormEJA(B)
1331 sage: J.matrix_basis()
1348 # Abuse the check_field parameter to check that the entries of
1349 # out basis (in ambient coordinates) are in the field QQ.
1350 if not all( all(b_i
in QQ
for b_i
in b
.list()) for b
in basis
):
1351 raise TypeError("basis not rational")
1353 self
._rational
_algebra
= None
1355 # There's no point in constructing the extra algebra if this
1356 # one is already rational.
1358 # Note: the same Jordan and inner-products work here,
1359 # because they are necessarily defined with respect to
1360 # ambient coordinates and not any particular basis.
1361 self
._rational
_algebra
= FiniteDimensionalEJA(
1366 orthonormalize
=False,
1370 super().__init
__(basis
,
1374 check_field
=check_field
,
1378 def _charpoly_coefficients(self
):
1382 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
1383 ....: JordanSpinEJA)
1387 The base ring of the resulting polynomial coefficients is what
1388 it should be, and not the rationals (unless the algebra was
1389 already over the rationals)::
1391 sage: J = JordanSpinEJA(3)
1392 sage: J._charpoly_coefficients()
1393 (X1^2 - X2^2 - X3^2, -2*X1)
1394 sage: a0 = J._charpoly_coefficients()[0]
1396 Algebraic Real Field
1397 sage: a0.base_ring()
1398 Algebraic Real Field
1401 if self
._rational
_algebra
is None:
1402 # There's no need to construct *another* algebra over the
1403 # rationals if this one is already over the
1404 # rationals. Likewise, if we never orthonormalized our
1405 # basis, we might as well just use the given one.
1406 return super()._charpoly
_coefficients
()
1408 # Do the computation over the rationals. The answer will be
1409 # the same, because all we've done is a change of basis.
1410 # Then, change back from QQ to our real base ring
1411 a
= ( a_i
.change_ring(self
.base_ring())
1412 for a_i
in self
._rational
_algebra
._charpoly
_coefficients
() )
1414 if self
._deortho
_matrix
is None:
1415 # This can happen if our base ring was, say, AA and we
1416 # chose not to (or didn't need to) orthonormalize. It's
1417 # still faster to do the computations over QQ even if
1418 # the numbers in the boxes stay the same.
1421 # Otherwise, convert the coordinate variables back to the
1422 # deorthonormalized ones.
1423 R
= self
.coordinate_polynomial_ring()
1424 from sage
.modules
.free_module_element
import vector
1425 X
= vector(R
, R
.gens())
1426 BX
= self
._deortho
_matrix
*X
1428 subs_dict
= { X[i]: BX[i] for i in range(len(X)) }
1429 return tuple( a_i
.subs(subs_dict
) for a_i
in a
)
1431 class ConcreteEJA(RationalBasisEJA
):
1433 A class for the Euclidean Jordan algebras that we know by name.
1435 These are the Jordan algebras whose basis, multiplication table,
1436 rank, and so on are known a priori. More to the point, they are
1437 the Euclidean Jordan algebras for which we are able to conjure up
1438 a "random instance."
1442 sage: from mjo.eja.eja_algebra import ConcreteEJA
1446 Our basis is normalized with respect to the algebra's inner
1447 product, unless we specify otherwise::
1449 sage: set_random_seed()
1450 sage: J = ConcreteEJA.random_instance()
1451 sage: all( b.norm() == 1 for b in J.gens() )
1454 Since our basis is orthonormal with respect to the algebra's inner
1455 product, and since we know that this algebra is an EJA, any
1456 left-multiplication operator's matrix will be symmetric because
1457 natural->EJA basis representation is an isometry and within the
1458 EJA the operator is self-adjoint by the Jordan axiom::
1460 sage: set_random_seed()
1461 sage: J = ConcreteEJA.random_instance()
1462 sage: x = J.random_element()
1463 sage: x.operator().is_self_adjoint()
1468 def _max_random_instance_size():
1470 Return an integer "size" that is an upper bound on the size of
1471 this algebra when it is used in a random test
1472 case. Unfortunately, the term "size" is ambiguous -- when
1473 dealing with `R^n` under either the Hadamard or Jordan spin
1474 product, the "size" refers to the dimension `n`. When dealing
1475 with a matrix algebra (real symmetric or complex/quaternion
1476 Hermitian), it refers to the size of the matrix, which is far
1477 less than the dimension of the underlying vector space.
1479 This method must be implemented in each subclass.
1481 raise NotImplementedError
1484 def random_instance(cls
, *args
, **kwargs
):
1486 Return a random instance of this type of algebra.
1488 This method should be implemented in each subclass.
1490 from sage
.misc
.prandom
import choice
1491 eja_class
= choice(cls
.__subclasses
__())
1493 # These all bubble up to the RationalBasisEJA superclass
1494 # constructor, so any (kw)args valid there are also valid
1496 return eja_class
.random_instance(*args
, **kwargs
)
1501 def dimension_over_reals():
1503 The dimension of this matrix's base ring over the reals.
1505 The reals are dimension one over themselves, obviously; that's
1506 just `\mathbb{R}^{1}`. Likewise, the complex numbers `a + bi`
1507 have dimension two. Finally, the quaternions have dimension
1508 four over the reals.
1510 This is used to determine the size of the matrix returned from
1511 :meth:`real_embed`, among other things.
1513 raise NotImplementedError
1516 def real_embed(cls
,M
):
1518 Embed the matrix ``M`` into a space of real matrices.
1520 The matrix ``M`` can have entries in any field at the moment:
1521 the real numbers, complex numbers, or quaternions. And although
1522 they are not a field, we can probably support octonions at some
1523 point, too. This function returns a real matrix that "acts like"
1524 the original with respect to matrix multiplication; i.e.
1526 real_embed(M*N) = real_embed(M)*real_embed(N)
1529 if M
.ncols() != M
.nrows():
1530 raise ValueError("the matrix 'M' must be square")
1535 def real_unembed(cls
,M
):
1537 The inverse of :meth:`real_embed`.
1539 if M
.ncols() != M
.nrows():
1540 raise ValueError("the matrix 'M' must be square")
1541 if not ZZ(M
.nrows()).mod(cls
.dimension_over_reals()).is_zero():
1542 raise ValueError("the matrix 'M' must be a real embedding")
1546 def jordan_product(X
,Y
):
1547 return (X
*Y
+ Y
*X
)/2
1550 def trace_inner_product(cls
,X
,Y
):
1552 Compute the trace inner-product of two real-embeddings.
1556 sage: from mjo.eja.eja_algebra import (RealSymmetricEJA,
1557 ....: ComplexHermitianEJA,
1558 ....: QuaternionHermitianEJA)
1562 This gives the same answer as it would if we computed the trace
1563 from the unembedded (original) matrices::
1565 sage: set_random_seed()
1566 sage: J = RealSymmetricEJA.random_instance()
1567 sage: x,y = J.random_elements(2)
1568 sage: Xe = x.to_matrix()
1569 sage: Ye = y.to_matrix()
1570 sage: X = J.real_unembed(Xe)
1571 sage: Y = J.real_unembed(Ye)
1572 sage: expected = (X*Y).trace()
1573 sage: actual = J.trace_inner_product(Xe,Ye)
1574 sage: actual == expected
1579 sage: set_random_seed()
1580 sage: J = ComplexHermitianEJA.random_instance()
1581 sage: x,y = J.random_elements(2)
1582 sage: Xe = x.to_matrix()
1583 sage: Ye = y.to_matrix()
1584 sage: X = J.real_unembed(Xe)
1585 sage: Y = J.real_unembed(Ye)
1586 sage: expected = (X*Y).trace().real()
1587 sage: actual = J.trace_inner_product(Xe,Ye)
1588 sage: actual == expected
1593 sage: set_random_seed()
1594 sage: J = QuaternionHermitianEJA.random_instance()
1595 sage: x,y = J.random_elements(2)
1596 sage: Xe = x.to_matrix()
1597 sage: Ye = y.to_matrix()
1598 sage: X = J.real_unembed(Xe)
1599 sage: Y = J.real_unembed(Ye)
1600 sage: expected = (X*Y).trace().coefficient_tuple()[0]
1601 sage: actual = J.trace_inner_product(Xe,Ye)
1602 sage: actual == expected
1606 Xu
= cls
.real_unembed(X
)
1607 Yu
= cls
.real_unembed(Y
)
1608 tr
= (Xu
*Yu
).trace()
1611 # Works in QQ, AA, RDF, et cetera.
1613 except AttributeError:
1614 # A quaternion doesn't have a real() method, but does
1615 # have coefficient_tuple() method that returns the
1616 # coefficients of 1, i, j, and k -- in that order.
1617 return tr
.coefficient_tuple()[0]
1620 class RealMatrixEJA(MatrixEJA
):
1622 def dimension_over_reals():
1626 class RealSymmetricEJA(ConcreteEJA
, RealMatrixEJA
):
1628 The rank-n simple EJA consisting of real symmetric n-by-n
1629 matrices, the usual symmetric Jordan product, and the trace inner
1630 product. It has dimension `(n^2 + n)/2` over the reals.
1634 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1638 sage: J = RealSymmetricEJA(2)
1639 sage: e0, e1, e2 = J.gens()
1647 In theory, our "field" can be any subfield of the reals::
1649 sage: RealSymmetricEJA(2, field=RDF)
1650 Euclidean Jordan algebra of dimension 3 over Real Double Field
1651 sage: RealSymmetricEJA(2, field=RR)
1652 Euclidean Jordan algebra of dimension 3 over Real Field with
1653 53 bits of precision
1657 The dimension of this algebra is `(n^2 + n) / 2`::
1659 sage: set_random_seed()
1660 sage: n_max = RealSymmetricEJA._max_random_instance_size()
1661 sage: n = ZZ.random_element(1, n_max)
1662 sage: J = RealSymmetricEJA(n)
1663 sage: J.dimension() == (n^2 + n)/2
1666 The Jordan multiplication is what we think it is::
1668 sage: set_random_seed()
1669 sage: J = RealSymmetricEJA.random_instance()
1670 sage: x,y = J.random_elements(2)
1671 sage: actual = (x*y).to_matrix()
1672 sage: X = x.to_matrix()
1673 sage: Y = y.to_matrix()
1674 sage: expected = (X*Y + Y*X)/2
1675 sage: actual == expected
1677 sage: J(expected) == x*y
1680 We can change the generator prefix::
1682 sage: RealSymmetricEJA(3, prefix='q').gens()
1683 (q0, q1, q2, q3, q4, q5)
1685 We can construct the (trivial) algebra of rank zero::
1687 sage: RealSymmetricEJA(0)
1688 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1692 def _denormalized_basis(cls
, n
):
1694 Return a basis for the space of real symmetric n-by-n matrices.
1698 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1702 sage: set_random_seed()
1703 sage: n = ZZ.random_element(1,5)
1704 sage: B = RealSymmetricEJA._denormalized_basis(n)
1705 sage: all( M.is_symmetric() for M in B)
1709 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1713 for j
in range(i
+1):
1714 Eij
= matrix(ZZ
, n
, lambda k
,l
: k
==i
and l
==j
)
1718 Sij
= Eij
+ Eij
.transpose()
1724 def _max_random_instance_size():
1725 return 4 # Dimension 10
1728 def random_instance(cls
, **kwargs
):
1730 Return a random instance of this type of algebra.
1732 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1733 return cls(n
, **kwargs
)
1735 def __init__(self
, n
, **kwargs
):
1736 # We know this is a valid EJA, but will double-check
1737 # if the user passes check_axioms=True.
1738 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
1740 super(RealSymmetricEJA
, self
).__init
__(self
._denormalized
_basis
(n
),
1741 self
.jordan_product
,
1742 self
.trace_inner_product
,
1745 # TODO: this could be factored out somehow, but is left here
1746 # because the MatrixEJA is not presently a subclass of the
1747 # FDEJA class that defines rank() and one().
1748 self
.rank
.set_cache(n
)
1749 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
1750 self
.one
.set_cache(self(idV
))
1754 class ComplexMatrixEJA(MatrixEJA
):
1755 # A manual dictionary-cache for the complex_extension() method,
1756 # since apparently @classmethods can't also be @cached_methods.
1757 _complex_extension
= {}
1760 def complex_extension(cls
,field
):
1762 The complex field that we embed/unembed, as an extension
1763 of the given ``field``.
1765 if field
in cls
._complex
_extension
:
1766 return cls
._complex
_extension
[field
]
1768 # Sage doesn't know how to adjoin the complex "i" (the root of
1769 # x^2 + 1) to a field in a general way. Here, we just enumerate
1770 # all of the cases that I have cared to support so far.
1772 # Sage doesn't know how to embed AA into QQbar, i.e. how
1773 # to adjoin sqrt(-1) to AA.
1775 elif not field
.is_exact():
1777 F
= field
.complex_field()
1779 # Works for QQ and... maybe some other fields.
1780 R
= PolynomialRing(field
, 'z')
1782 F
= field
.extension(z
**2 + 1, 'I', embedding
=CLF(-1).sqrt())
1784 cls
._complex
_extension
[field
] = F
1788 def dimension_over_reals():
1792 def real_embed(cls
,M
):
1794 Embed the n-by-n complex matrix ``M`` into the space of real
1795 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1796 bi` to the block matrix ``[[a,b],[-b,a]]``.
1800 sage: from mjo.eja.eja_algebra import ComplexMatrixEJA
1804 sage: F = QuadraticField(-1, 'I')
1805 sage: x1 = F(4 - 2*i)
1806 sage: x2 = F(1 + 2*i)
1809 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1810 sage: ComplexMatrixEJA.real_embed(M)
1819 Embedding is a homomorphism (isomorphism, in fact)::
1821 sage: set_random_seed()
1822 sage: n = ZZ.random_element(3)
1823 sage: F = QuadraticField(-1, 'I')
1824 sage: X = random_matrix(F, n)
1825 sage: Y = random_matrix(F, n)
1826 sage: Xe = ComplexMatrixEJA.real_embed(X)
1827 sage: Ye = ComplexMatrixEJA.real_embed(Y)
1828 sage: XYe = ComplexMatrixEJA.real_embed(X*Y)
1833 super(ComplexMatrixEJA
,cls
).real_embed(M
)
1836 # We don't need any adjoined elements...
1837 field
= M
.base_ring().base_ring()
1843 blocks
.append(matrix(field
, 2, [ [ a
, b
],
1846 return matrix
.block(field
, n
, blocks
)
1850 def real_unembed(cls
,M
):
1852 The inverse of _embed_complex_matrix().
1856 sage: from mjo.eja.eja_algebra import ComplexMatrixEJA
1860 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1861 ....: [-2, 1, -4, 3],
1862 ....: [ 9, 10, 11, 12],
1863 ....: [-10, 9, -12, 11] ])
1864 sage: ComplexMatrixEJA.real_unembed(A)
1866 [ 10*I + 9 12*I + 11]
1870 Unembedding is the inverse of embedding::
1872 sage: set_random_seed()
1873 sage: F = QuadraticField(-1, 'I')
1874 sage: M = random_matrix(F, 3)
1875 sage: Me = ComplexMatrixEJA.real_embed(M)
1876 sage: ComplexMatrixEJA.real_unembed(Me) == M
1880 super(ComplexMatrixEJA
,cls
).real_unembed(M
)
1882 d
= cls
.dimension_over_reals()
1883 F
= cls
.complex_extension(M
.base_ring())
1886 # Go top-left to bottom-right (reading order), converting every
1887 # 2-by-2 block we see to a single complex element.
1889 for k
in range(n
/d
):
1890 for j
in range(n
/d
):
1891 submat
= M
[d
*k
:d
*k
+d
,d
*j
:d
*j
+d
]
1892 if submat
[0,0] != submat
[1,1]:
1893 raise ValueError('bad on-diagonal submatrix')
1894 if submat
[0,1] != -submat
[1,0]:
1895 raise ValueError('bad off-diagonal submatrix')
1896 z
= submat
[0,0] + submat
[0,1]*i
1899 return matrix(F
, n
/d
, elements
)
1902 class ComplexHermitianEJA(ConcreteEJA
, ComplexMatrixEJA
):
1904 The rank-n simple EJA consisting of complex Hermitian n-by-n
1905 matrices over the real numbers, the usual symmetric Jordan product,
1906 and the real-part-of-trace inner product. It has dimension `n^2` over
1911 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1915 In theory, our "field" can be any subfield of the reals::
1917 sage: ComplexHermitianEJA(2, field=RDF)
1918 Euclidean Jordan algebra of dimension 4 over Real Double Field
1919 sage: ComplexHermitianEJA(2, field=RR)
1920 Euclidean Jordan algebra of dimension 4 over Real Field with
1921 53 bits of precision
1925 The dimension of this algebra is `n^2`::
1927 sage: set_random_seed()
1928 sage: n_max = ComplexHermitianEJA._max_random_instance_size()
1929 sage: n = ZZ.random_element(1, n_max)
1930 sage: J = ComplexHermitianEJA(n)
1931 sage: J.dimension() == n^2
1934 The Jordan multiplication is what we think it is::
1936 sage: set_random_seed()
1937 sage: J = ComplexHermitianEJA.random_instance()
1938 sage: x,y = J.random_elements(2)
1939 sage: actual = (x*y).to_matrix()
1940 sage: X = x.to_matrix()
1941 sage: Y = y.to_matrix()
1942 sage: expected = (X*Y + Y*X)/2
1943 sage: actual == expected
1945 sage: J(expected) == x*y
1948 We can change the generator prefix::
1950 sage: ComplexHermitianEJA(2, prefix='z').gens()
1953 We can construct the (trivial) algebra of rank zero::
1955 sage: ComplexHermitianEJA(0)
1956 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1961 def _denormalized_basis(cls
, n
):
1963 Returns a basis for the space of complex Hermitian n-by-n matrices.
1965 Why do we embed these? Basically, because all of numerical linear
1966 algebra assumes that you're working with vectors consisting of `n`
1967 entries from a field and scalars from the same field. There's no way
1968 to tell SageMath that (for example) the vectors contain complex
1969 numbers, while the scalar field is real.
1973 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1977 sage: set_random_seed()
1978 sage: n = ZZ.random_element(1,5)
1979 sage: B = ComplexHermitianEJA._denormalized_basis(n)
1980 sage: all( M.is_symmetric() for M in B)
1985 R
= PolynomialRing(field
, 'z')
1987 F
= field
.extension(z
**2 + 1, 'I')
1990 # This is like the symmetric case, but we need to be careful:
1992 # * We want conjugate-symmetry, not just symmetry.
1993 # * The diagonal will (as a result) be real.
1996 Eij
= matrix
.zero(F
,n
)
1998 for j
in range(i
+1):
2002 Sij
= cls
.real_embed(Eij
)
2005 # The second one has a minus because it's conjugated.
2006 Eij
[j
,i
] = 1 # Eij = Eij + Eij.transpose()
2007 Sij_real
= cls
.real_embed(Eij
)
2009 # Eij = I*Eij - I*Eij.transpose()
2012 Sij_imag
= cls
.real_embed(Eij
)
2018 # Since we embedded these, we can drop back to the "field" that we
2019 # started with instead of the complex extension "F".
2020 return tuple( s
.change_ring(field
) for s
in S
)
2023 def __init__(self
, n
, **kwargs
):
2024 # We know this is a valid EJA, but will double-check
2025 # if the user passes check_axioms=True.
2026 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2028 super(ComplexHermitianEJA
, self
).__init
__(self
._denormalized
_basis
(n
),
2029 self
.jordan_product
,
2030 self
.trace_inner_product
,
2032 # TODO: this could be factored out somehow, but is left here
2033 # because the MatrixEJA is not presently a subclass of the
2034 # FDEJA class that defines rank() and one().
2035 self
.rank
.set_cache(n
)
2036 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
2037 self
.one
.set_cache(self(idV
))
2040 def _max_random_instance_size():
2041 return 3 # Dimension 9
2044 def random_instance(cls
, **kwargs
):
2046 Return a random instance of this type of algebra.
2048 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2049 return cls(n
, **kwargs
)
2051 class QuaternionMatrixEJA(MatrixEJA
):
2053 # A manual dictionary-cache for the quaternion_extension() method,
2054 # since apparently @classmethods can't also be @cached_methods.
2055 _quaternion_extension
= {}
2058 def quaternion_extension(cls
,field
):
2060 The quaternion field that we embed/unembed, as an extension
2061 of the given ``field``.
2063 if field
in cls
._quaternion
_extension
:
2064 return cls
._quaternion
_extension
[field
]
2066 Q
= QuaternionAlgebra(field
,-1,-1)
2068 cls
._quaternion
_extension
[field
] = Q
2072 def dimension_over_reals():
2076 def real_embed(cls
,M
):
2078 Embed the n-by-n quaternion matrix ``M`` into the space of real
2079 matrices of size 4n-by-4n by first sending each quaternion entry `z
2080 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
2081 c+di],[-c + di, a-bi]]`, and then embedding those into a real
2086 sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA
2090 sage: Q = QuaternionAlgebra(QQ,-1,-1)
2091 sage: i,j,k = Q.gens()
2092 sage: x = 1 + 2*i + 3*j + 4*k
2093 sage: M = matrix(Q, 1, [[x]])
2094 sage: QuaternionMatrixEJA.real_embed(M)
2100 Embedding is a homomorphism (isomorphism, in fact)::
2102 sage: set_random_seed()
2103 sage: n = ZZ.random_element(2)
2104 sage: Q = QuaternionAlgebra(QQ,-1,-1)
2105 sage: X = random_matrix(Q, n)
2106 sage: Y = random_matrix(Q, n)
2107 sage: Xe = QuaternionMatrixEJA.real_embed(X)
2108 sage: Ye = QuaternionMatrixEJA.real_embed(Y)
2109 sage: XYe = QuaternionMatrixEJA.real_embed(X*Y)
2114 super(QuaternionMatrixEJA
,cls
).real_embed(M
)
2115 quaternions
= M
.base_ring()
2118 F
= QuadraticField(-1, 'I')
2123 t
= z
.coefficient_tuple()
2128 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
2129 [-c
+ d
*i
, a
- b
*i
]])
2130 realM
= ComplexMatrixEJA
.real_embed(cplxM
)
2131 blocks
.append(realM
)
2133 # We should have real entries by now, so use the realest field
2134 # we've got for the return value.
2135 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
2140 def real_unembed(cls
,M
):
2142 The inverse of _embed_quaternion_matrix().
2146 sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA
2150 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
2151 ....: [-2, 1, -4, 3],
2152 ....: [-3, 4, 1, -2],
2153 ....: [-4, -3, 2, 1]])
2154 sage: QuaternionMatrixEJA.real_unembed(M)
2155 [1 + 2*i + 3*j + 4*k]
2159 Unembedding is the inverse of embedding::
2161 sage: set_random_seed()
2162 sage: Q = QuaternionAlgebra(QQ, -1, -1)
2163 sage: M = random_matrix(Q, 3)
2164 sage: Me = QuaternionMatrixEJA.real_embed(M)
2165 sage: QuaternionMatrixEJA.real_unembed(Me) == M
2169 super(QuaternionMatrixEJA
,cls
).real_unembed(M
)
2171 d
= cls
.dimension_over_reals()
2173 # Use the base ring of the matrix to ensure that its entries can be
2174 # multiplied by elements of the quaternion algebra.
2175 Q
= cls
.quaternion_extension(M
.base_ring())
2178 # Go top-left to bottom-right (reading order), converting every
2179 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
2182 for l
in range(n
/d
):
2183 for m
in range(n
/d
):
2184 submat
= ComplexMatrixEJA
.real_unembed(
2185 M
[d
*l
:d
*l
+d
,d
*m
:d
*m
+d
] )
2186 if submat
[0,0] != submat
[1,1].conjugate():
2187 raise ValueError('bad on-diagonal submatrix')
2188 if submat
[0,1] != -submat
[1,0].conjugate():
2189 raise ValueError('bad off-diagonal submatrix')
2190 z
= submat
[0,0].real()
2191 z
+= submat
[0,0].imag()*i
2192 z
+= submat
[0,1].real()*j
2193 z
+= submat
[0,1].imag()*k
2196 return matrix(Q
, n
/d
, elements
)
2199 class QuaternionHermitianEJA(ConcreteEJA
, QuaternionMatrixEJA
):
2201 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
2202 matrices, the usual symmetric Jordan product, and the
2203 real-part-of-trace inner product. It has dimension `2n^2 - n` over
2208 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2212 In theory, our "field" can be any subfield of the reals::
2214 sage: QuaternionHermitianEJA(2, field=RDF)
2215 Euclidean Jordan algebra of dimension 6 over Real Double Field
2216 sage: QuaternionHermitianEJA(2, field=RR)
2217 Euclidean Jordan algebra of dimension 6 over Real Field with
2218 53 bits of precision
2222 The dimension of this algebra is `2*n^2 - n`::
2224 sage: set_random_seed()
2225 sage: n_max = QuaternionHermitianEJA._max_random_instance_size()
2226 sage: n = ZZ.random_element(1, n_max)
2227 sage: J = QuaternionHermitianEJA(n)
2228 sage: J.dimension() == 2*(n^2) - n
2231 The Jordan multiplication is what we think it is::
2233 sage: set_random_seed()
2234 sage: J = QuaternionHermitianEJA.random_instance()
2235 sage: x,y = J.random_elements(2)
2236 sage: actual = (x*y).to_matrix()
2237 sage: X = x.to_matrix()
2238 sage: Y = y.to_matrix()
2239 sage: expected = (X*Y + Y*X)/2
2240 sage: actual == expected
2242 sage: J(expected) == x*y
2245 We can change the generator prefix::
2247 sage: QuaternionHermitianEJA(2, prefix='a').gens()
2248 (a0, a1, a2, a3, a4, a5)
2250 We can construct the (trivial) algebra of rank zero::
2252 sage: QuaternionHermitianEJA(0)
2253 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2257 def _denormalized_basis(cls
, n
):
2259 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
2261 Why do we embed these? Basically, because all of numerical
2262 linear algebra assumes that you're working with vectors consisting
2263 of `n` entries from a field and scalars from the same field. There's
2264 no way to tell SageMath that (for example) the vectors contain
2265 complex numbers, while the scalar field is real.
2269 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2273 sage: set_random_seed()
2274 sage: n = ZZ.random_element(1,5)
2275 sage: B = QuaternionHermitianEJA._denormalized_basis(n)
2276 sage: all( M.is_symmetric() for M in B )
2281 Q
= QuaternionAlgebra(QQ
,-1,-1)
2284 # This is like the symmetric case, but we need to be careful:
2286 # * We want conjugate-symmetry, not just symmetry.
2287 # * The diagonal will (as a result) be real.
2290 Eij
= matrix
.zero(Q
,n
)
2292 for j
in range(i
+1):
2296 Sij
= cls
.real_embed(Eij
)
2299 # The second, third, and fourth ones have a minus
2300 # because they're conjugated.
2301 # Eij = Eij + Eij.transpose()
2303 Sij_real
= cls
.real_embed(Eij
)
2305 # Eij = I*(Eij - Eij.transpose())
2308 Sij_I
= cls
.real_embed(Eij
)
2310 # Eij = J*(Eij - Eij.transpose())
2313 Sij_J
= cls
.real_embed(Eij
)
2315 # Eij = K*(Eij - Eij.transpose())
2318 Sij_K
= cls
.real_embed(Eij
)
2324 # Since we embedded these, we can drop back to the "field" that we
2325 # started with instead of the quaternion algebra "Q".
2326 return tuple( s
.change_ring(field
) for s
in S
)
2329 def __init__(self
, n
, **kwargs
):
2330 # We know this is a valid EJA, but will double-check
2331 # if the user passes check_axioms=True.
2332 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2334 super(QuaternionHermitianEJA
, self
).__init
__(self
._denormalized
_basis
(n
),
2335 self
.jordan_product
,
2336 self
.trace_inner_product
,
2338 # TODO: this could be factored out somehow, but is left here
2339 # because the MatrixEJA is not presently a subclass of the
2340 # FDEJA class that defines rank() and one().
2341 self
.rank
.set_cache(n
)
2342 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
2343 self
.one
.set_cache(self(idV
))
2347 def _max_random_instance_size():
2349 The maximum rank of a random QuaternionHermitianEJA.
2351 return 2 # Dimension 6
2354 def random_instance(cls
, **kwargs
):
2356 Return a random instance of this type of algebra.
2358 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2359 return cls(n
, **kwargs
)
2362 class HadamardEJA(ConcreteEJA
):
2364 Return the Euclidean Jordan Algebra corresponding to the set
2365 `R^n` under the Hadamard product.
2367 Note: this is nothing more than the Cartesian product of ``n``
2368 copies of the spin algebra. Once Cartesian product algebras
2369 are implemented, this can go.
2373 sage: from mjo.eja.eja_algebra import HadamardEJA
2377 This multiplication table can be verified by hand::
2379 sage: J = HadamardEJA(3)
2380 sage: e0,e1,e2 = J.gens()
2396 We can change the generator prefix::
2398 sage: HadamardEJA(3, prefix='r').gens()
2402 def __init__(self
, n
, **kwargs
):
2404 jordan_product
= lambda x
,y
: x
2405 inner_product
= lambda x
,y
: x
2407 def jordan_product(x
,y
):
2409 return P( xi
*yi
for (xi
,yi
) in zip(x
,y
) )
2411 def inner_product(x
,y
):
2414 # New defaults for keyword arguments. Don't orthonormalize
2415 # because our basis is already orthonormal with respect to our
2416 # inner-product. Don't check the axioms, because we know this
2417 # is a valid EJA... but do double-check if the user passes
2418 # check_axioms=True. Note: we DON'T override the "check_field"
2419 # default here, because the user can pass in a field!
2420 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2421 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2423 column_basis
= tuple( b
.column() for b
in FreeModule(ZZ
, n
).basis() )
2424 super().__init
__(column_basis
,
2429 self
.rank
.set_cache(n
)
2432 self
.one
.set_cache( self
.zero() )
2434 self
.one
.set_cache( sum(self
.gens()) )
2437 def _max_random_instance_size():
2439 The maximum dimension of a random HadamardEJA.
2444 def random_instance(cls
, **kwargs
):
2446 Return a random instance of this type of algebra.
2448 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2449 return cls(n
, **kwargs
)
2452 class BilinearFormEJA(ConcreteEJA
):
2454 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2455 with the half-trace inner product and jordan product ``x*y =
2456 (<Bx,y>,y_bar>, x0*y_bar + y0*x_bar)`` where `B = 1 \times B22` is
2457 a symmetric positive-definite "bilinear form" matrix. Its
2458 dimension is the size of `B`, and it has rank two in dimensions
2459 larger than two. It reduces to the ``JordanSpinEJA`` when `B` is
2460 the identity matrix of order ``n``.
2462 We insist that the one-by-one upper-left identity block of `B` be
2463 passed in as well so that we can be passed a matrix of size zero
2464 to construct a trivial algebra.
2468 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
2469 ....: JordanSpinEJA)
2473 When no bilinear form is specified, the identity matrix is used,
2474 and the resulting algebra is the Jordan spin algebra::
2476 sage: B = matrix.identity(AA,3)
2477 sage: J0 = BilinearFormEJA(B)
2478 sage: J1 = JordanSpinEJA(3)
2479 sage: J0.multiplication_table() == J0.multiplication_table()
2482 An error is raised if the matrix `B` does not correspond to a
2483 positive-definite bilinear form::
2485 sage: B = matrix.random(QQ,2,3)
2486 sage: J = BilinearFormEJA(B)
2487 Traceback (most recent call last):
2489 ValueError: bilinear form is not positive-definite
2490 sage: B = matrix.zero(QQ,3)
2491 sage: J = BilinearFormEJA(B)
2492 Traceback (most recent call last):
2494 ValueError: bilinear form is not positive-definite
2498 We can create a zero-dimensional algebra::
2500 sage: B = matrix.identity(AA,0)
2501 sage: J = BilinearFormEJA(B)
2505 We can check the multiplication condition given in the Jordan, von
2506 Neumann, and Wigner paper (and also discussed on my "On the
2507 symmetry..." paper). Note that this relies heavily on the standard
2508 choice of basis, as does anything utilizing the bilinear form
2509 matrix. We opt not to orthonormalize the basis, because if we
2510 did, we would have to normalize the `s_{i}` in a similar manner::
2512 sage: set_random_seed()
2513 sage: n = ZZ.random_element(5)
2514 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2515 sage: B11 = matrix.identity(QQ,1)
2516 sage: B22 = M.transpose()*M
2517 sage: B = block_matrix(2,2,[ [B11,0 ],
2519 sage: J = BilinearFormEJA(B, orthonormalize=False)
2520 sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
2521 sage: V = J.vector_space()
2522 sage: sis = [ J( V([0] + (M.inverse()*ei).list()).column() )
2523 ....: for ei in eis ]
2524 sage: actual = [ sis[i]*sis[j]
2525 ....: for i in range(n-1)
2526 ....: for j in range(n-1) ]
2527 sage: expected = [ J.one() if i == j else J.zero()
2528 ....: for i in range(n-1)
2529 ....: for j in range(n-1) ]
2530 sage: actual == expected
2534 def __init__(self
, B
, **kwargs
):
2535 # The matrix "B" is supplied by the user in most cases,
2536 # so it makes sense to check whether or not its positive-
2537 # definite unless we are specifically asked not to...
2538 if ("check_axioms" not in kwargs
) or kwargs
["check_axioms"]:
2539 if not B
.is_positive_definite():
2540 raise ValueError("bilinear form is not positive-definite")
2542 # However, all of the other data for this EJA is computed
2543 # by us in manner that guarantees the axioms are
2544 # satisfied. So, again, unless we are specifically asked to
2545 # verify things, we'll skip the rest of the checks.
2546 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2548 def inner_product(x
,y
):
2549 return (y
.T
*B
*x
)[0,0]
2551 def jordan_product(x
,y
):
2557 z0
= inner_product(y
,x
)
2558 zbar
= y0
*xbar
+ x0
*ybar
2559 return P([z0
] + zbar
.list())
2562 column_basis
= tuple( b
.column() for b
in FreeModule(ZZ
, n
).basis() )
2563 super(BilinearFormEJA
, self
).__init
__(column_basis
,
2568 # The rank of this algebra is two, unless we're in a
2569 # one-dimensional ambient space (because the rank is bounded
2570 # by the ambient dimension).
2571 self
.rank
.set_cache(min(n
,2))
2574 self
.one
.set_cache( self
.zero() )
2576 self
.one
.set_cache( self
.monomial(0) )
2579 def _max_random_instance_size():
2581 The maximum dimension of a random BilinearFormEJA.
2586 def random_instance(cls
, **kwargs
):
2588 Return a random instance of this algebra.
2590 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2592 B
= matrix
.identity(ZZ
, n
)
2593 return cls(B
, **kwargs
)
2595 B11
= matrix
.identity(ZZ
, 1)
2596 M
= matrix
.random(ZZ
, n
-1)
2597 I
= matrix
.identity(ZZ
, n
-1)
2599 while alpha
.is_zero():
2600 alpha
= ZZ
.random_element().abs()
2601 B22
= M
.transpose()*M
+ alpha
*I
2603 from sage
.matrix
.special
import block_matrix
2604 B
= block_matrix(2,2, [ [B11
, ZZ(0) ],
2607 return cls(B
, **kwargs
)
2610 class JordanSpinEJA(BilinearFormEJA
):
2612 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2613 with the usual inner product and jordan product ``x*y =
2614 (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2619 sage: from mjo.eja.eja_algebra import JordanSpinEJA
2623 This multiplication table can be verified by hand::
2625 sage: J = JordanSpinEJA(4)
2626 sage: e0,e1,e2,e3 = J.gens()
2642 We can change the generator prefix::
2644 sage: JordanSpinEJA(2, prefix='B').gens()
2649 Ensure that we have the usual inner product on `R^n`::
2651 sage: set_random_seed()
2652 sage: J = JordanSpinEJA.random_instance()
2653 sage: x,y = J.random_elements(2)
2654 sage: actual = x.inner_product(y)
2655 sage: expected = x.to_vector().inner_product(y.to_vector())
2656 sage: actual == expected
2660 def __init__(self
, n
, **kwargs
):
2661 # This is a special case of the BilinearFormEJA with the
2662 # identity matrix as its bilinear form.
2663 B
= matrix
.identity(ZZ
, n
)
2665 # Don't orthonormalize because our basis is already
2666 # orthonormal with respect to our inner-product.
2667 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2669 # But also don't pass check_field=False here, because the user
2670 # can pass in a field!
2671 super(JordanSpinEJA
, self
).__init
__(B
, **kwargs
)
2674 def _max_random_instance_size():
2676 The maximum dimension of a random JordanSpinEJA.
2681 def random_instance(cls
, **kwargs
):
2683 Return a random instance of this type of algebra.
2685 Needed here to override the implementation for ``BilinearFormEJA``.
2687 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2688 return cls(n
, **kwargs
)
2691 class TrivialEJA(ConcreteEJA
):
2693 The trivial Euclidean Jordan algebra consisting of only a zero element.
2697 sage: from mjo.eja.eja_algebra import TrivialEJA
2701 sage: J = TrivialEJA()
2708 sage: 7*J.one()*12*J.one()
2710 sage: J.one().inner_product(J.one())
2712 sage: J.one().norm()
2714 sage: J.one().subalgebra_generated_by()
2715 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2720 def __init__(self
, **kwargs
):
2721 jordan_product
= lambda x
,y
: x
2722 inner_product
= lambda x
,y
: 0
2725 # New defaults for keyword arguments
2726 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2727 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2729 super(TrivialEJA
, self
).__init
__(basis
,
2733 # The rank is zero using my definition, namely the dimension of the
2734 # largest subalgebra generated by any element.
2735 self
.rank
.set_cache(0)
2736 self
.one
.set_cache( self
.zero() )
2739 def random_instance(cls
, **kwargs
):
2740 # We don't take a "size" argument so the superclass method is
2741 # inappropriate for us.
2742 return cls(**kwargs
)
2745 class CartesianProductEJA(CombinatorialFreeModule_CartesianProduct
,
2746 FiniteDimensionalEJA
):
2748 The external (orthogonal) direct sum of two or more Euclidean
2749 Jordan algebras. Every Euclidean Jordan algebra decomposes into an
2750 orthogonal direct sum of simple Euclidean Jordan algebras which is
2751 then isometric to a Cartesian product, so no generality is lost by
2752 providing only this construction.
2756 sage: from mjo.eja.eja_algebra import (random_eja,
2757 ....: CartesianProductEJA,
2759 ....: JordanSpinEJA,
2760 ....: RealSymmetricEJA)
2764 The Jordan product is inherited from our factors and implemented by
2765 our CombinatorialFreeModule Cartesian product superclass::
2767 sage: set_random_seed()
2768 sage: J1 = HadamardEJA(2)
2769 sage: J2 = RealSymmetricEJA(2)
2770 sage: J = cartesian_product([J1,J2])
2771 sage: x,y = J.random_elements(2)
2775 The ability to retrieve the original factors is implemented by our
2776 CombinatorialFreeModule Cartesian product superclass::
2778 sage: J1 = HadamardEJA(2, field=QQ)
2779 sage: J2 = JordanSpinEJA(3, field=QQ)
2780 sage: J = cartesian_product([J1,J2])
2781 sage: J.cartesian_factors()
2782 (Euclidean Jordan algebra of dimension 2 over Rational Field,
2783 Euclidean Jordan algebra of dimension 3 over Rational Field)
2785 You can provide more than two factors::
2787 sage: J1 = HadamardEJA(2)
2788 sage: J2 = JordanSpinEJA(3)
2789 sage: J3 = RealSymmetricEJA(3)
2790 sage: cartesian_product([J1,J2,J3])
2791 Euclidean Jordan algebra of dimension 2 over Algebraic Real
2792 Field (+) Euclidean Jordan algebra of dimension 3 over Algebraic
2793 Real Field (+) Euclidean Jordan algebra of dimension 6 over
2794 Algebraic Real Field
2796 Rank is additive on a Cartesian product::
2798 sage: J1 = HadamardEJA(1)
2799 sage: J2 = RealSymmetricEJA(2)
2800 sage: J = cartesian_product([J1,J2])
2801 sage: J1.rank.clear_cache()
2802 sage: J2.rank.clear_cache()
2803 sage: J.rank.clear_cache()
2806 sage: J.rank() == J1.rank() + J2.rank()
2809 The same rank computation works over the rationals, with whatever
2812 sage: J1 = HadamardEJA(1, field=QQ, orthonormalize=False)
2813 sage: J2 = RealSymmetricEJA(2, field=QQ, orthonormalize=False)
2814 sage: J = cartesian_product([J1,J2])
2815 sage: J1.rank.clear_cache()
2816 sage: J2.rank.clear_cache()
2817 sage: J.rank.clear_cache()
2820 sage: J.rank() == J1.rank() + J2.rank()
2823 The product algebra will be associative if and only if all of its
2824 components are associative::
2826 sage: J1 = HadamardEJA(2)
2827 sage: J1.is_associative()
2829 sage: J2 = HadamardEJA(3)
2830 sage: J2.is_associative()
2832 sage: J3 = RealSymmetricEJA(3)
2833 sage: J3.is_associative()
2835 sage: CP1 = cartesian_product([J1,J2])
2836 sage: CP1.is_associative()
2838 sage: CP2 = cartesian_product([J1,J3])
2839 sage: CP2.is_associative()
2844 All factors must share the same base field::
2846 sage: J1 = HadamardEJA(2, field=QQ)
2847 sage: J2 = RealSymmetricEJA(2)
2848 sage: CartesianProductEJA((J1,J2))
2849 Traceback (most recent call last):
2851 ValueError: all factors must share the same base field
2853 The cached unit element is the same one that would be computed::
2855 sage: set_random_seed() # long time
2856 sage: J1 = random_eja() # long time
2857 sage: J2 = random_eja() # long time
2858 sage: J = cartesian_product([J1,J2]) # long time
2859 sage: actual = J.one() # long time
2860 sage: J.one.clear_cache() # long time
2861 sage: expected = J.one() # long time
2862 sage: actual == expected # long time
2866 Element
= FiniteDimensionalEJAElement
2869 def __init__(self
, algebras
, **kwargs
):
2870 CombinatorialFreeModule_CartesianProduct
.__init
__(self
,
2873 field
= algebras
[0].base_ring()
2874 if not all( J
.base_ring() == field
for J
in algebras
):
2875 raise ValueError("all factors must share the same base field")
2877 associative
= all( m
.is_associative() for m
in algebras
)
2879 # The definition of matrix_space() and self.basis() relies
2880 # only on the stuff in the CFM_CartesianProduct class, which
2881 # we've already initialized.
2882 Js
= self
.cartesian_factors()
2884 MS
= self
.matrix_space()
2886 MS(tuple( self
.cartesian_projection(i
)(b
).to_matrix()
2887 for i
in range(m
) ))
2888 for b
in self
.basis()
2891 # Define jordan/inner products that operate on that matrix_basis.
2892 def jordan_product(x
,y
):
2894 (Js
[i
](x
[i
])*Js
[i
](y
[i
])).to_matrix() for i
in range(m
)
2897 def inner_product(x
, y
):
2899 Js
[i
](x
[i
]).inner_product(Js
[i
](y
[i
])) for i
in range(m
)
2902 # There's no need to check the field since it already came
2903 # from an EJA. Likewise the axioms are guaranteed to be
2904 # satisfied, unless the guy writing this class sucks.
2906 # If you want the basis to be orthonormalized, orthonormalize
2908 FiniteDimensionalEJA
.__init
__(self
,
2913 orthonormalize
=False,
2914 associative
=associative
,
2915 cartesian_product
=True,
2919 ones
= tuple(J
.one() for J
in algebras
)
2920 self
.one
.set_cache(self
._cartesian
_product
_of
_elements
(ones
))
2921 self
.rank
.set_cache(sum(J
.rank() for J
in algebras
))
2923 def matrix_space(self
):
2925 Return the space that our matrix basis lives in as a Cartesian
2930 sage: from mjo.eja.eja_algebra import (HadamardEJA,
2931 ....: RealSymmetricEJA)
2935 sage: J1 = HadamardEJA(1)
2936 sage: J2 = RealSymmetricEJA(2)
2937 sage: J = cartesian_product([J1,J2])
2938 sage: J.matrix_space()
2939 The Cartesian product of (Full MatrixSpace of 1 by 1 dense
2940 matrices over Algebraic Real Field, Full MatrixSpace of 2
2941 by 2 dense matrices over Algebraic Real Field)
2944 from sage
.categories
.cartesian_product
import cartesian_product
2945 return cartesian_product( [J
.matrix_space()
2946 for J
in self
.cartesian_factors()] )
2949 def cartesian_projection(self
, i
):
2953 sage: from mjo.eja.eja_algebra import (random_eja,
2954 ....: JordanSpinEJA,
2956 ....: RealSymmetricEJA,
2957 ....: ComplexHermitianEJA)
2961 The projection morphisms are Euclidean Jordan algebra
2964 sage: J1 = HadamardEJA(2)
2965 sage: J2 = RealSymmetricEJA(2)
2966 sage: J = cartesian_product([J1,J2])
2967 sage: J.cartesian_projection(0)
2968 Linear operator between finite-dimensional Euclidean Jordan
2969 algebras represented by the matrix:
2972 Domain: Euclidean Jordan algebra of dimension 2 over Algebraic
2973 Real Field (+) Euclidean Jordan algebra of dimension 3 over
2974 Algebraic Real Field
2975 Codomain: Euclidean Jordan algebra of dimension 2 over Algebraic
2977 sage: J.cartesian_projection(1)
2978 Linear operator between finite-dimensional Euclidean Jordan
2979 algebras represented by the matrix:
2983 Domain: Euclidean Jordan algebra of dimension 2 over Algebraic
2984 Real Field (+) Euclidean Jordan algebra of dimension 3 over
2985 Algebraic Real Field
2986 Codomain: Euclidean Jordan algebra of dimension 3 over Algebraic
2989 The projections work the way you'd expect on the vector
2990 representation of an element::
2992 sage: J1 = JordanSpinEJA(2)
2993 sage: J2 = ComplexHermitianEJA(2)
2994 sage: J = cartesian_product([J1,J2])
2995 sage: pi_left = J.cartesian_projection(0)
2996 sage: pi_right = J.cartesian_projection(1)
2997 sage: pi_left(J.one()).to_vector()
2999 sage: pi_right(J.one()).to_vector()
3001 sage: J.one().to_vector()
3006 The answer never changes::
3008 sage: set_random_seed()
3009 sage: J1 = random_eja()
3010 sage: J2 = random_eja()
3011 sage: J = cartesian_product([J1,J2])
3012 sage: P0 = J.cartesian_projection(0)
3013 sage: P1 = J.cartesian_projection(0)
3018 Ji
= self
.cartesian_factors()[i
]
3019 # Requires the fix on Trac 31421/31422 to work!
3020 Pi
= super().cartesian_projection(i
)
3021 return FiniteDimensionalEJAOperator(self
,Ji
,Pi
.matrix())
3024 def cartesian_embedding(self
, i
):
3028 sage: from mjo.eja.eja_algebra import (random_eja,
3029 ....: JordanSpinEJA,
3031 ....: RealSymmetricEJA)
3035 The embedding morphisms are Euclidean Jordan algebra
3038 sage: J1 = HadamardEJA(2)
3039 sage: J2 = RealSymmetricEJA(2)
3040 sage: J = cartesian_product([J1,J2])
3041 sage: J.cartesian_embedding(0)
3042 Linear operator between finite-dimensional Euclidean Jordan
3043 algebras represented by the matrix:
3049 Domain: Euclidean Jordan algebra of dimension 2 over
3050 Algebraic Real Field
3051 Codomain: Euclidean Jordan algebra of dimension 2 over
3052 Algebraic Real Field (+) Euclidean Jordan algebra of
3053 dimension 3 over Algebraic Real Field
3054 sage: J.cartesian_embedding(1)
3055 Linear operator between finite-dimensional Euclidean Jordan
3056 algebras represented by the matrix:
3062 Domain: Euclidean Jordan algebra of dimension 3 over
3063 Algebraic Real Field
3064 Codomain: Euclidean Jordan algebra of dimension 2 over
3065 Algebraic Real Field (+) Euclidean Jordan algebra of
3066 dimension 3 over Algebraic Real Field
3068 The embeddings work the way you'd expect on the vector
3069 representation of an element::
3071 sage: J1 = JordanSpinEJA(3)
3072 sage: J2 = RealSymmetricEJA(2)
3073 sage: J = cartesian_product([J1,J2])
3074 sage: iota_left = J.cartesian_embedding(0)
3075 sage: iota_right = J.cartesian_embedding(1)
3076 sage: iota_left(J1.zero()) == J.zero()
3078 sage: iota_right(J2.zero()) == J.zero()
3080 sage: J1.one().to_vector()
3082 sage: iota_left(J1.one()).to_vector()
3084 sage: J2.one().to_vector()
3086 sage: iota_right(J2.one()).to_vector()
3088 sage: J.one().to_vector()
3093 The answer never changes::
3095 sage: set_random_seed()
3096 sage: J1 = random_eja()
3097 sage: J2 = random_eja()
3098 sage: J = cartesian_product([J1,J2])
3099 sage: E0 = J.cartesian_embedding(0)
3100 sage: E1 = J.cartesian_embedding(0)
3104 Composing a projection with the corresponding inclusion should
3105 produce the identity map, and mismatching them should produce
3108 sage: set_random_seed()
3109 sage: J1 = random_eja()
3110 sage: J2 = random_eja()
3111 sage: J = cartesian_product([J1,J2])
3112 sage: iota_left = J.cartesian_embedding(0)
3113 sage: iota_right = J.cartesian_embedding(1)
3114 sage: pi_left = J.cartesian_projection(0)
3115 sage: pi_right = J.cartesian_projection(1)
3116 sage: pi_left*iota_left == J1.one().operator()
3118 sage: pi_right*iota_right == J2.one().operator()
3120 sage: (pi_left*iota_right).is_zero()
3122 sage: (pi_right*iota_left).is_zero()
3126 Ji
= self
.cartesian_factors()[i
]
3127 # Requires the fix on Trac 31421/31422 to work!
3128 Ei
= super().cartesian_embedding(i
)
3129 return FiniteDimensionalEJAOperator(Ji
,self
,Ei
.matrix())
3133 FiniteDimensionalEJA
.CartesianProduct
= CartesianProductEJA
3135 random_eja
= ConcreteEJA
.random_instance
3136 #def random_eja(*args, **kwargs):
3137 # from sage.categories.cartesian_product import cartesian_product
3138 # J1 = HadamardEJA(1, **kwargs)
3139 # J2 = RealSymmetricEJA(2, **kwargs)
3140 # J = cartesian_product([J1,J2])