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 _mat2vec
38 class FiniteDimensionalEJA(CombinatorialFreeModule
):
40 A finite-dimensional Euclidean Jordan algebra.
44 - basis -- a tuple of basis elements in their matrix form.
46 - jordan_product -- function of two elements (in matrix form)
47 that returns their jordan product in this algebra; this will
48 be applied to ``basis`` to compute a multiplication table for
51 - inner_product -- function of two elements (in matrix form) that
52 returns their inner product. This will be applied to ``basis`` to
53 compute an inner-product table (basically a matrix) for this algebra.
56 Element
= FiniteDimensionalEJAElement
70 if not field
.is_subring(RR
):
71 # Note: this does return true for the real algebraic
72 # field, the rationals, and any quadratic field where
73 # we've specified a real embedding.
74 raise ValueError("scalar field is not real")
76 # If the basis given to us wasn't over the field that it's
77 # supposed to be over, fix that. Or, you know, crash.
78 basis
= tuple( b
.change_ring(field
) for b
in basis
)
81 # Check commutativity of the Jordan and inner-products.
82 # This has to be done before we build the multiplication
83 # and inner-product tables/matrices, because we take
84 # advantage of symmetry in the process.
85 if not all( jordan_product(bi
,bj
) == jordan_product(bj
,bi
)
88 raise ValueError("Jordan product is not commutative")
90 if not all( inner_product(bi
,bj
) == inner_product(bj
,bi
)
93 raise ValueError("inner-product is not commutative")
96 category
= MagmaticAlgebras(field
).FiniteDimensional()
97 category
= category
.WithBasis().Unital()
99 # Element subalgebras can take advantage of this.
100 category
= category
.Associative()
102 # Call the superclass constructor so that we can use its from_vector()
103 # method to build our multiplication table.
105 super().__init
__(field
,
111 # Now comes all of the hard work. We'll be constructing an
112 # ambient vector space V that our (vectorized) basis lives in,
113 # as well as a subspace W of V spanned by those (vectorized)
114 # basis elements. The W-coordinates are the coefficients that
115 # we see in things like x = 1*e1 + 2*e2.
120 # Works on both column and square matrices...
121 degree
= len(basis
[0].list())
123 # Build an ambient space that fits our matrix basis when
124 # written out as "long vectors."
125 V
= VectorSpace(field
, degree
)
127 # The matrix that will hole the orthonormal -> unorthonormal
128 # coordinate transformation.
129 self
._deortho
_matrix
= None
132 # Save a copy of the un-orthonormalized basis for later.
133 # Convert it to ambient V (vector) coordinates while we're
134 # at it, because we'd have to do it later anyway.
135 deortho_vector_basis
= tuple( V(b
.list()) for b
in basis
)
137 from mjo
.eja
.eja_utils
import gram_schmidt
138 basis
= tuple(gram_schmidt(basis
, inner_product
))
140 # Save the (possibly orthonormalized) matrix basis for
142 self
._matrix
_basis
= basis
144 # Now create the vector space for the algebra, which will have
145 # its own set of non-ambient coordinates (in terms of the
147 vector_basis
= tuple( V(b
.list()) for b
in basis
)
148 W
= V
.span_of_basis( vector_basis
, check
=check_axioms
)
151 # Now "W" is the vector space of our algebra coordinates. The
152 # variables "X1", "X2",... refer to the entries of vectors in
153 # W. Thus to convert back and forth between the orthonormal
154 # coordinates and the given ones, we need to stick the original
156 U
= V
.span_of_basis( deortho_vector_basis
, check
=check_axioms
)
157 self
._deortho
_matrix
= matrix( U
.coordinate_vector(q
)
158 for q
in vector_basis
)
161 # Now we actually compute the multiplication and inner-product
162 # tables/matrices using the possibly-orthonormalized basis.
163 self
._inner
_product
_matrix
= matrix
.identity(field
, n
)
164 self
._multiplication
_table
= [ [0 for j
in range(i
+1)]
167 # Note: the Jordan and inner-products are defined in terms
168 # of the ambient basis. It's important that their arguments
169 # are in ambient coordinates as well.
172 # ortho basis w.r.t. ambient coords
176 # The jordan product returns a matrixy answer, so we
177 # have to convert it to the algebra coordinates.
178 elt
= jordan_product(q_i
, q_j
)
179 elt
= W
.coordinate_vector(V(elt
.list()))
180 self
._multiplication
_table
[i
][j
] = self
.from_vector(elt
)
182 if not orthonormalize
:
183 # If we're orthonormalizing the basis with respect
184 # to an inner-product, then the inner-product
185 # matrix with respect to the resulting basis is
186 # just going to be the identity.
187 ip
= inner_product(q_i
, q_j
)
188 self
._inner
_product
_matrix
[i
,j
] = ip
189 self
._inner
_product
_matrix
[j
,i
] = ip
191 self
._inner
_product
_matrix
._cache
= {'hermitian': True}
192 self
._inner
_product
_matrix
.set_immutable()
195 if not self
._is
_jordanian
():
196 raise ValueError("Jordan identity does not hold")
197 if not self
._inner
_product
_is
_associative
():
198 raise ValueError("inner product is not associative")
201 def _coerce_map_from_base_ring(self
):
203 Disable the map from the base ring into the algebra.
205 Performing a nonsense conversion like this automatically
206 is counterpedagogical. The fallback is to try the usual
207 element constructor, which should also fail.
211 sage: from mjo.eja.eja_algebra import random_eja
215 sage: set_random_seed()
216 sage: J = random_eja()
218 Traceback (most recent call last):
220 ValueError: not an element of this algebra
226 def product_on_basis(self
, i
, j
):
227 # We only stored the lower-triangular portion of the
228 # multiplication table.
230 return self
._multiplication
_table
[i
][j
]
232 return self
._multiplication
_table
[j
][i
]
234 def inner_product(self
, x
, y
):
236 The inner product associated with this Euclidean Jordan algebra.
238 Defaults to the trace inner product, but can be overridden by
239 subclasses if they are sure that the necessary properties are
244 sage: from mjo.eja.eja_algebra import (random_eja,
246 ....: BilinearFormEJA)
250 Our inner product is "associative," which means the following for
251 a symmetric bilinear form::
253 sage: set_random_seed()
254 sage: J = random_eja()
255 sage: x,y,z = J.random_elements(3)
256 sage: (x*y).inner_product(z) == y.inner_product(x*z)
261 Ensure that this is the usual inner product for the algebras
264 sage: set_random_seed()
265 sage: J = HadamardEJA.random_instance()
266 sage: x,y = J.random_elements(2)
267 sage: actual = x.inner_product(y)
268 sage: expected = x.to_vector().inner_product(y.to_vector())
269 sage: actual == expected
272 Ensure that this is one-half of the trace inner-product in a
273 BilinearFormEJA that isn't just the reals (when ``n`` isn't
274 one). This is in Faraut and Koranyi, and also my "On the
277 sage: set_random_seed()
278 sage: J = BilinearFormEJA.random_instance()
279 sage: n = J.dimension()
280 sage: x = J.random_element()
281 sage: y = J.random_element()
282 sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
285 B
= self
._inner
_product
_matrix
286 return (B
*x
.to_vector()).inner_product(y
.to_vector())
289 def _is_commutative(self
):
291 Whether or not this algebra's multiplication table is commutative.
293 This method should of course always return ``True``, unless
294 this algebra was constructed with ``check_axioms=False`` and
295 passed an invalid multiplication table.
297 return all( self
.product_on_basis(i
,j
) == self
.product_on_basis(i
,j
)
298 for i
in range(self
.dimension())
299 for j
in range(self
.dimension()) )
301 def _is_jordanian(self
):
303 Whether or not this algebra's multiplication table respects the
304 Jordan identity `(x^{2})(xy) = x(x^{2}y)`.
306 We only check one arrangement of `x` and `y`, so for a
307 ``True`` result to be truly true, you should also check
308 :meth:`_is_commutative`. This method should of course always
309 return ``True``, unless this algebra was constructed with
310 ``check_axioms=False`` and passed an invalid multiplication table.
312 return all( (self
.monomial(i
)**2)*(self
.monomial(i
)*self
.monomial(j
))
314 (self
.monomial(i
))*((self
.monomial(i
)**2)*self
.monomial(j
))
315 for i
in range(self
.dimension())
316 for j
in range(self
.dimension()) )
318 def _inner_product_is_associative(self
):
320 Return whether or not this algebra's inner product `B` is
321 associative; that is, whether or not `B(xy,z) = B(x,yz)`.
323 This method should of course always return ``True``, unless
324 this algebra was constructed with ``check_axioms=False`` and
325 passed an invalid Jordan or inner-product.
328 # Used to check whether or not something is zero in an inexact
329 # ring. This number is sufficient to allow the construction of
330 # QuaternionHermitianEJA(2, field=RDF) with check_axioms=True.
333 for i
in range(self
.dimension()):
334 for j
in range(self
.dimension()):
335 for k
in range(self
.dimension()):
339 diff
= (x
*y
).inner_product(z
) - x
.inner_product(y
*z
)
341 if self
.base_ring().is_exact():
345 if diff
.abs() > epsilon
:
350 def _element_constructor_(self
, elt
):
352 Construct an element of this algebra from its vector or matrix
355 This gets called only after the parent element _call_ method
356 fails to find a coercion for the argument.
360 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
362 ....: RealSymmetricEJA)
366 The identity in `S^n` is converted to the identity in the EJA::
368 sage: J = RealSymmetricEJA(3)
369 sage: I = matrix.identity(QQ,3)
370 sage: J(I) == J.one()
373 This skew-symmetric matrix can't be represented in the EJA::
375 sage: J = RealSymmetricEJA(3)
376 sage: A = matrix(QQ,3, lambda i,j: i-j)
378 Traceback (most recent call last):
380 ValueError: not an element of this algebra
384 Ensure that we can convert any element of the two non-matrix
385 simple algebras (whose matrix representations are columns)
386 back and forth faithfully::
388 sage: set_random_seed()
389 sage: J = HadamardEJA.random_instance()
390 sage: x = J.random_element()
391 sage: J(x.to_vector().column()) == x
393 sage: J = JordanSpinEJA.random_instance()
394 sage: x = J.random_element()
395 sage: J(x.to_vector().column()) == x
399 msg
= "not an element of this algebra"
401 # The superclass implementation of random_element()
402 # needs to be able to coerce "0" into the algebra.
404 elif elt
in self
.base_ring():
405 # Ensure that no base ring -> algebra coercion is performed
406 # by this method. There's some stupidity in sage that would
407 # otherwise propagate to this method; for example, sage thinks
408 # that the integer 3 belongs to the space of 2-by-2 matrices.
409 raise ValueError(msg
)
413 except (AttributeError, TypeError):
414 # Try to convert a vector into a column-matrix
417 if elt
not in self
.matrix_space():
418 raise ValueError(msg
)
420 # Thanks for nothing! Matrix spaces aren't vector spaces in
421 # Sage, so we have to figure out its matrix-basis coordinates
422 # ourselves. We use the basis space's ring instead of the
423 # element's ring because the basis space might be an algebraic
424 # closure whereas the base ring of the 3-by-3 identity matrix
425 # could be QQ instead of QQbar.
427 # We pass check=False because the matrix basis is "guaranteed"
428 # to be linearly independent... right? Ha ha.
429 V
= VectorSpace(self
.base_ring(), elt
.nrows()*elt
.ncols())
430 W
= V
.span_of_basis( (_mat2vec(s
) for s
in self
.matrix_basis()),
434 coords
= W
.coordinate_vector(_mat2vec(elt
))
435 except ArithmeticError: # vector is not in free module
436 raise ValueError(msg
)
438 return self
.from_vector(coords
)
442 Return a string representation of ``self``.
446 sage: from mjo.eja.eja_algebra import JordanSpinEJA
450 Ensure that it says what we think it says::
452 sage: JordanSpinEJA(2, field=AA)
453 Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
454 sage: JordanSpinEJA(3, field=RDF)
455 Euclidean Jordan algebra of dimension 3 over Real Double Field
458 fmt
= "Euclidean Jordan algebra of dimension {} over {}"
459 return fmt
.format(self
.dimension(), self
.base_ring())
463 def characteristic_polynomial_of(self
):
465 Return the algebra's "characteristic polynomial of" function,
466 which is itself a multivariate polynomial that, when evaluated
467 at the coordinates of some algebra element, returns that
468 element's characteristic polynomial.
470 The resulting polynomial has `n+1` variables, where `n` is the
471 dimension of this algebra. The first `n` variables correspond to
472 the coordinates of an algebra element: when evaluated at the
473 coordinates of an algebra element with respect to a certain
474 basis, the result is a univariate polynomial (in the one
475 remaining variable ``t``), namely the characteristic polynomial
480 sage: from mjo.eja.eja_algebra import JordanSpinEJA, TrivialEJA
484 The characteristic polynomial in the spin algebra is given in
485 Alizadeh, Example 11.11::
487 sage: J = JordanSpinEJA(3)
488 sage: p = J.characteristic_polynomial_of(); p
489 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
490 sage: xvec = J.one().to_vector()
494 By definition, the characteristic polynomial is a monic
495 degree-zero polynomial in a rank-zero algebra. Note that
496 Cayley-Hamilton is indeed satisfied since the polynomial
497 ``1`` evaluates to the identity element of the algebra on
500 sage: J = TrivialEJA()
501 sage: J.characteristic_polynomial_of()
508 # The list of coefficient polynomials a_0, a_1, a_2, ..., a_(r-1).
509 a
= self
._charpoly
_coefficients
()
511 # We go to a bit of trouble here to reorder the
512 # indeterminates, so that it's easier to evaluate the
513 # characteristic polynomial at x's coordinates and get back
514 # something in terms of t, which is what we want.
515 S
= PolynomialRing(self
.base_ring(),'t')
519 S
= PolynomialRing(S
, R
.variable_names())
522 return (t
**r
+ sum( a
[k
]*(t
**k
) for k
in range(r
) ))
524 def coordinate_polynomial_ring(self
):
526 The multivariate polynomial ring in which this algebra's
527 :meth:`characteristic_polynomial_of` lives.
531 sage: from mjo.eja.eja_algebra import (HadamardEJA,
532 ....: RealSymmetricEJA)
536 sage: J = HadamardEJA(2)
537 sage: J.coordinate_polynomial_ring()
538 Multivariate Polynomial Ring in X1, X2...
539 sage: J = RealSymmetricEJA(3,field=QQ,orthonormalize=False)
540 sage: J.coordinate_polynomial_ring()
541 Multivariate Polynomial Ring in X1, X2, X3, X4, X5, X6...
544 var_names
= tuple( "X%d" % z
for z
in range(1, self
.dimension()+1) )
545 return PolynomialRing(self
.base_ring(), var_names
)
547 def inner_product(self
, x
, y
):
549 The inner product associated with this Euclidean Jordan algebra.
551 Defaults to the trace inner product, but can be overridden by
552 subclasses if they are sure that the necessary properties are
557 sage: from mjo.eja.eja_algebra import (random_eja,
559 ....: BilinearFormEJA)
563 Our inner product is "associative," which means the following for
564 a symmetric bilinear form::
566 sage: set_random_seed()
567 sage: J = random_eja()
568 sage: x,y,z = J.random_elements(3)
569 sage: (x*y).inner_product(z) == y.inner_product(x*z)
574 Ensure that this is the usual inner product for the algebras
577 sage: set_random_seed()
578 sage: J = HadamardEJA.random_instance()
579 sage: x,y = J.random_elements(2)
580 sage: actual = x.inner_product(y)
581 sage: expected = x.to_vector().inner_product(y.to_vector())
582 sage: actual == expected
585 Ensure that this is one-half of the trace inner-product in a
586 BilinearFormEJA that isn't just the reals (when ``n`` isn't
587 one). This is in Faraut and Koranyi, and also my "On the
590 sage: set_random_seed()
591 sage: J = BilinearFormEJA.random_instance()
592 sage: n = J.dimension()
593 sage: x = J.random_element()
594 sage: y = J.random_element()
595 sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
598 B
= self
._inner
_product
_matrix
599 return (B
*x
.to_vector()).inner_product(y
.to_vector())
602 def is_trivial(self
):
604 Return whether or not this algebra is trivial.
606 A trivial algebra contains only the zero element.
610 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
615 sage: J = ComplexHermitianEJA(3)
621 sage: J = TrivialEJA()
626 return self
.dimension() == 0
629 def multiplication_table(self
):
631 Return a visual representation of this algebra's multiplication
632 table (on basis elements).
636 sage: from mjo.eja.eja_algebra import JordanSpinEJA
640 sage: J = JordanSpinEJA(4)
641 sage: J.multiplication_table()
642 +----++----+----+----+----+
643 | * || e0 | e1 | e2 | e3 |
644 +====++====+====+====+====+
645 | e0 || e0 | e1 | e2 | e3 |
646 +----++----+----+----+----+
647 | e1 || e1 | e0 | 0 | 0 |
648 +----++----+----+----+----+
649 | e2 || e2 | 0 | e0 | 0 |
650 +----++----+----+----+----+
651 | e3 || e3 | 0 | 0 | e0 |
652 +----++----+----+----+----+
656 # Prepend the header row.
657 M
= [["*"] + list(self
.gens())]
659 # And to each subsequent row, prepend an entry that belongs to
660 # the left-side "header column."
661 M
+= [ [self
.monomial(i
)] + [ self
.product_on_basis(i
,j
)
665 return table(M
, header_row
=True, header_column
=True, frame
=True)
668 def matrix_basis(self
):
670 Return an (often more natural) representation of this algebras
671 basis as an ordered tuple of matrices.
673 Every finite-dimensional Euclidean Jordan Algebra is a, up to
674 Jordan isomorphism, a direct sum of five simple
675 algebras---four of which comprise Hermitian matrices. And the
676 last type of algebra can of course be thought of as `n`-by-`1`
677 column matrices (ambiguusly called column vectors) to avoid
678 special cases. As a result, matrices (and column vectors) are
679 a natural representation format for Euclidean Jordan algebra
682 But, when we construct an algebra from a basis of matrices,
683 those matrix representations are lost in favor of coordinate
684 vectors *with respect to* that basis. We could eventually
685 convert back if we tried hard enough, but having the original
686 representations handy is valuable enough that we simply store
687 them and return them from this method.
689 Why implement this for non-matrix algebras? Avoiding special
690 cases for the :class:`BilinearFormEJA` pays with simplicity in
691 its own right. But mainly, we would like to be able to assume
692 that elements of a :class:`CartesianProductEJA` can be displayed
693 nicely, without having to have special classes for direct sums
694 one of whose components was a matrix algebra.
698 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
699 ....: RealSymmetricEJA)
703 sage: J = RealSymmetricEJA(2)
705 Finite family {0: e0, 1: e1, 2: e2}
706 sage: J.matrix_basis()
708 [1 0] [ 0 0.7071067811865475?] [0 0]
709 [0 0], [0.7071067811865475? 0], [0 1]
714 sage: J = JordanSpinEJA(2)
716 Finite family {0: e0, 1: e1}
717 sage: J.matrix_basis()
723 return self
._matrix
_basis
726 def matrix_space(self
):
728 Return the matrix space in which this algebra's elements live, if
729 we think of them as matrices (including column vectors of the
732 Generally this will be an `n`-by-`1` column-vector space,
733 except when the algebra is trivial. There it's `n`-by-`n`
734 (where `n` is zero), to ensure that two elements of the matrix
735 space (empty matrices) can be multiplied.
737 Matrix algebras override this with something more useful.
739 if self
.is_trivial():
740 return MatrixSpace(self
.base_ring(), 0)
742 return self
.matrix_basis()[0].parent()
748 Return the unit element of this algebra.
752 sage: from mjo.eja.eja_algebra import (HadamardEJA,
757 We can compute unit element in the Hadamard EJA::
759 sage: J = HadamardEJA(5)
761 e0 + e1 + e2 + e3 + e4
763 The unit element in the Hadamard EJA is inherited in the
764 subalgebras generated by its elements::
766 sage: J = HadamardEJA(5)
768 e0 + e1 + e2 + e3 + e4
769 sage: x = sum(J.gens())
770 sage: A = x.subalgebra_generated_by(orthonormalize=False)
773 sage: A.one().superalgebra_element()
774 e0 + e1 + e2 + e3 + e4
778 The identity element acts like the identity, regardless of
779 whether or not we orthonormalize::
781 sage: set_random_seed()
782 sage: J = random_eja()
783 sage: x = J.random_element()
784 sage: J.one()*x == x and x*J.one() == x
786 sage: A = x.subalgebra_generated_by()
787 sage: y = A.random_element()
788 sage: A.one()*y == y and y*A.one() == y
793 sage: set_random_seed()
794 sage: J = random_eja(field=QQ, orthonormalize=False)
795 sage: x = J.random_element()
796 sage: J.one()*x == x and x*J.one() == x
798 sage: A = x.subalgebra_generated_by(orthonormalize=False)
799 sage: y = A.random_element()
800 sage: A.one()*y == y and y*A.one() == y
803 The matrix of the unit element's operator is the identity,
804 regardless of the base field and whether or not we
807 sage: set_random_seed()
808 sage: J = random_eja()
809 sage: actual = J.one().operator().matrix()
810 sage: expected = matrix.identity(J.base_ring(), J.dimension())
811 sage: actual == expected
813 sage: x = J.random_element()
814 sage: A = x.subalgebra_generated_by()
815 sage: actual = A.one().operator().matrix()
816 sage: expected = matrix.identity(A.base_ring(), A.dimension())
817 sage: actual == expected
822 sage: set_random_seed()
823 sage: J = random_eja(field=QQ, orthonormalize=False)
824 sage: actual = J.one().operator().matrix()
825 sage: expected = matrix.identity(J.base_ring(), J.dimension())
826 sage: actual == expected
828 sage: x = J.random_element()
829 sage: A = x.subalgebra_generated_by(orthonormalize=False)
830 sage: actual = A.one().operator().matrix()
831 sage: expected = matrix.identity(A.base_ring(), A.dimension())
832 sage: actual == expected
835 Ensure that the cached unit element (often precomputed by
836 hand) agrees with the computed one::
838 sage: set_random_seed()
839 sage: J = random_eja()
840 sage: cached = J.one()
841 sage: J.one.clear_cache()
842 sage: J.one() == cached
847 sage: set_random_seed()
848 sage: J = random_eja(field=QQ, orthonormalize=False)
849 sage: cached = J.one()
850 sage: J.one.clear_cache()
851 sage: J.one() == cached
855 # We can brute-force compute the matrices of the operators
856 # that correspond to the basis elements of this algebra.
857 # If some linear combination of those basis elements is the
858 # algebra identity, then the same linear combination of
859 # their matrices has to be the identity matrix.
861 # Of course, matrices aren't vectors in sage, so we have to
862 # appeal to the "long vectors" isometry.
863 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
865 # Now we use basic linear algebra to find the coefficients,
866 # of the matrices-as-vectors-linear-combination, which should
867 # work for the original algebra basis too.
868 A
= matrix(self
.base_ring(), oper_vecs
)
870 # We used the isometry on the left-hand side already, but we
871 # still need to do it for the right-hand side. Recall that we
872 # wanted something that summed to the identity matrix.
873 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
875 # Now if there's an identity element in the algebra, this
876 # should work. We solve on the left to avoid having to
877 # transpose the matrix "A".
878 return self
.from_vector(A
.solve_left(b
))
881 def peirce_decomposition(self
, c
):
883 The Peirce decomposition of this algebra relative to the
886 In the future, this can be extended to a complete system of
887 orthogonal idempotents.
891 - ``c`` -- an idempotent of this algebra.
895 A triple (J0, J5, J1) containing two subalgebras and one subspace
898 - ``J0`` -- the algebra on the eigenspace of ``c.operator()``
899 corresponding to the eigenvalue zero.
901 - ``J5`` -- the eigenspace (NOT a subalgebra) of ``c.operator()``
902 corresponding to the eigenvalue one-half.
904 - ``J1`` -- the algebra on the eigenspace of ``c.operator()``
905 corresponding to the eigenvalue one.
907 These are the only possible eigenspaces for that operator, and this
908 algebra is a direct sum of them. The spaces ``J0`` and ``J1`` are
909 orthogonal, and are subalgebras of this algebra with the appropriate
914 sage: from mjo.eja.eja_algebra import random_eja, RealSymmetricEJA
918 The canonical example comes from the symmetric matrices, which
919 decompose into diagonal and off-diagonal parts::
921 sage: J = RealSymmetricEJA(3)
922 sage: C = matrix(QQ, [ [1,0,0],
926 sage: J0,J5,J1 = J.peirce_decomposition(c)
928 Euclidean Jordan algebra of dimension 1...
930 Vector space of degree 6 and dimension 2...
932 Euclidean Jordan algebra of dimension 3...
933 sage: J0.one().to_matrix()
937 sage: orig_df = AA.options.display_format
938 sage: AA.options.display_format = 'radical'
939 sage: J.from_vector(J5.basis()[0]).to_matrix()
943 sage: J.from_vector(J5.basis()[1]).to_matrix()
947 sage: AA.options.display_format = orig_df
948 sage: J1.one().to_matrix()
955 Every algebra decomposes trivially with respect to its identity
958 sage: set_random_seed()
959 sage: J = random_eja()
960 sage: J0,J5,J1 = J.peirce_decomposition(J.one())
961 sage: J0.dimension() == 0 and J5.dimension() == 0
963 sage: J1.superalgebra() == J and J1.dimension() == J.dimension()
966 The decomposition is into eigenspaces, and its components are
967 therefore necessarily orthogonal. Moreover, the identity
968 elements in the two subalgebras are the projections onto their
969 respective subspaces of the superalgebra's identity element::
971 sage: set_random_seed()
972 sage: J = random_eja()
973 sage: x = J.random_element()
974 sage: if not J.is_trivial():
975 ....: while x.is_nilpotent():
976 ....: x = J.random_element()
977 sage: c = x.subalgebra_idempotent()
978 sage: J0,J5,J1 = J.peirce_decomposition(c)
980 sage: for (w,y,z) in zip(J0.basis(), J5.basis(), J1.basis()):
981 ....: w = w.superalgebra_element()
982 ....: y = J.from_vector(y)
983 ....: z = z.superalgebra_element()
984 ....: ipsum += w.inner_product(y).abs()
985 ....: ipsum += w.inner_product(z).abs()
986 ....: ipsum += y.inner_product(z).abs()
989 sage: J1(c) == J1.one()
991 sage: J0(J.one() - c) == J0.one()
995 if not c
.is_idempotent():
996 raise ValueError("element is not idempotent: %s" % c
)
998 from mjo
.eja
.eja_subalgebra
import FiniteDimensionalEJASubalgebra
1000 # Default these to what they should be if they turn out to be
1001 # trivial, because eigenspaces_left() won't return eigenvalues
1002 # corresponding to trivial spaces (e.g. it returns only the
1003 # eigenspace corresponding to lambda=1 if you take the
1004 # decomposition relative to the identity element).
1005 trivial
= FiniteDimensionalEJASubalgebra(self
, ())
1006 J0
= trivial
# eigenvalue zero
1007 J5
= VectorSpace(self
.base_ring(), 0) # eigenvalue one-half
1008 J1
= trivial
# eigenvalue one
1010 for (eigval
, eigspace
) in c
.operator().matrix().right_eigenspaces():
1011 if eigval
== ~
(self
.base_ring()(2)):
1014 gens
= tuple( self
.from_vector(b
) for b
in eigspace
.basis() )
1015 subalg
= FiniteDimensionalEJASubalgebra(self
,
1023 raise ValueError("unexpected eigenvalue: %s" % eigval
)
1028 def random_element(self
, thorough
=False):
1030 Return a random element of this algebra.
1032 Our algebra superclass method only returns a linear
1033 combination of at most two basis elements. We instead
1034 want the vector space "random element" method that
1035 returns a more diverse selection.
1039 - ``thorough`` -- (boolean; default False) whether or not we
1040 should generate irrational coefficients for the random
1041 element when our base ring is irrational; this slows the
1042 algebra operations to a crawl, but any truly random method
1046 # For a general base ring... maybe we can trust this to do the
1047 # right thing? Unlikely, but.
1048 V
= self
.vector_space()
1049 v
= V
.random_element()
1051 if self
.base_ring() is AA
:
1052 # The "random element" method of the algebraic reals is
1053 # stupid at the moment, and only returns integers between
1054 # -2 and 2, inclusive:
1056 # https://trac.sagemath.org/ticket/30875
1058 # Instead, we implement our own "random vector" method,
1059 # and then coerce that into the algebra. We use the vector
1060 # space degree here instead of the dimension because a
1061 # subalgebra could (for example) be spanned by only two
1062 # vectors, each with five coordinates. We need to
1063 # generate all five coordinates.
1065 v
*= QQbar
.random_element().real()
1067 v
*= QQ
.random_element()
1069 return self
.from_vector(V
.coordinate_vector(v
))
1071 def random_elements(self
, count
, thorough
=False):
1073 Return ``count`` random elements as a tuple.
1077 - ``thorough`` -- (boolean; default False) whether or not we
1078 should generate irrational coefficients for the random
1079 elements when our base ring is irrational; this slows the
1080 algebra operations to a crawl, but any truly random method
1085 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1089 sage: J = JordanSpinEJA(3)
1090 sage: x,y,z = J.random_elements(3)
1091 sage: all( [ x in J, y in J, z in J ])
1093 sage: len( J.random_elements(10) ) == 10
1097 return tuple( self
.random_element(thorough
)
1098 for idx
in range(count
) )
1102 def _charpoly_coefficients(self
):
1104 The `r` polynomial coefficients of the "characteristic polynomial
1109 sage: from mjo.eja.eja_algebra import random_eja
1113 The theory shows that these are all homogeneous polynomials of
1116 sage: set_random_seed()
1117 sage: J = random_eja()
1118 sage: all(p.is_homogeneous() for p in J._charpoly_coefficients())
1122 n
= self
.dimension()
1123 R
= self
.coordinate_polynomial_ring()
1125 F
= R
.fraction_field()
1128 # From a result in my book, these are the entries of the
1129 # basis representation of L_x.
1130 return sum( vars[k
]*self
.monomial(k
).operator().matrix()[i
,j
]
1133 L_x
= matrix(F
, n
, n
, L_x_i_j
)
1136 if self
.rank
.is_in_cache():
1138 # There's no need to pad the system with redundant
1139 # columns if we *know* they'll be redundant.
1142 # Compute an extra power in case the rank is equal to
1143 # the dimension (otherwise, we would stop at x^(r-1)).
1144 x_powers
= [ (L_x
**k
)*self
.one().to_vector()
1145 for k
in range(n
+1) ]
1146 A
= matrix
.column(F
, x_powers
[:n
])
1147 AE
= A
.extended_echelon_form()
1154 # The theory says that only the first "r" coefficients are
1155 # nonzero, and they actually live in the original polynomial
1156 # ring and not the fraction field. We negate them because in
1157 # the actual characteristic polynomial, they get moved to the
1158 # other side where x^r lives. We don't bother to trim A_rref
1159 # down to a square matrix and solve the resulting system,
1160 # because the upper-left r-by-r portion of A_rref is
1161 # guaranteed to be the identity matrix, so e.g.
1163 # A_rref.solve_right(Y)
1165 # would just be returning Y.
1166 return (-E
*b
)[:r
].change_ring(R
)
1171 Return the rank of this EJA.
1173 This is a cached method because we know the rank a priori for
1174 all of the algebras we can construct. Thus we can avoid the
1175 expensive ``_charpoly_coefficients()`` call unless we truly
1176 need to compute the whole characteristic polynomial.
1180 sage: from mjo.eja.eja_algebra import (HadamardEJA,
1181 ....: JordanSpinEJA,
1182 ....: RealSymmetricEJA,
1183 ....: ComplexHermitianEJA,
1184 ....: QuaternionHermitianEJA,
1189 The rank of the Jordan spin algebra is always two::
1191 sage: JordanSpinEJA(2).rank()
1193 sage: JordanSpinEJA(3).rank()
1195 sage: JordanSpinEJA(4).rank()
1198 The rank of the `n`-by-`n` Hermitian real, complex, or
1199 quaternion matrices is `n`::
1201 sage: RealSymmetricEJA(4).rank()
1203 sage: ComplexHermitianEJA(3).rank()
1205 sage: QuaternionHermitianEJA(2).rank()
1210 Ensure that every EJA that we know how to construct has a
1211 positive integer rank, unless the algebra is trivial in
1212 which case its rank will be zero::
1214 sage: set_random_seed()
1215 sage: J = random_eja()
1219 sage: r > 0 or (r == 0 and J.is_trivial())
1222 Ensure that computing the rank actually works, since the ranks
1223 of all simple algebras are known and will be cached by default::
1225 sage: set_random_seed() # long time
1226 sage: J = random_eja() # long time
1227 sage: cached = J.rank() # long time
1228 sage: J.rank.clear_cache() # long time
1229 sage: J.rank() == cached # long time
1233 return len(self
._charpoly
_coefficients
())
1236 def vector_space(self
):
1238 Return the vector space that underlies this algebra.
1242 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1246 sage: J = RealSymmetricEJA(2)
1247 sage: J.vector_space()
1248 Vector space of dimension 3 over...
1251 return self
.zero().to_vector().parent().ambient_vector_space()
1254 Element
= FiniteDimensionalEJAElement
1256 class RationalBasisEJA(FiniteDimensionalEJA
):
1258 New class for algebras whose supplied basis elements have all rational entries.
1262 sage: from mjo.eja.eja_algebra import BilinearFormEJA
1266 The supplied basis is orthonormalized by default::
1268 sage: B = matrix(QQ, [[1, 0, 0], [0, 25, -32], [0, -32, 41]])
1269 sage: J = BilinearFormEJA(B)
1270 sage: J.matrix_basis()
1287 # Abuse the check_field parameter to check that the entries of
1288 # out basis (in ambient coordinates) are in the field QQ.
1289 if not all( all(b_i
in QQ
for b_i
in b
.list()) for b
in basis
):
1290 raise TypeError("basis not rational")
1292 self
._rational
_algebra
= None
1294 # There's no point in constructing the extra algebra if this
1295 # one is already rational.
1297 # Note: the same Jordan and inner-products work here,
1298 # because they are necessarily defined with respect to
1299 # ambient coordinates and not any particular basis.
1300 self
._rational
_algebra
= FiniteDimensionalEJA(
1305 orthonormalize
=False,
1309 super().__init
__(basis
,
1313 check_field
=check_field
,
1317 def _charpoly_coefficients(self
):
1321 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
1322 ....: JordanSpinEJA)
1326 The base ring of the resulting polynomial coefficients is what
1327 it should be, and not the rationals (unless the algebra was
1328 already over the rationals)::
1330 sage: J = JordanSpinEJA(3)
1331 sage: J._charpoly_coefficients()
1332 (X1^2 - X2^2 - X3^2, -2*X1)
1333 sage: a0 = J._charpoly_coefficients()[0]
1335 Algebraic Real Field
1336 sage: a0.base_ring()
1337 Algebraic Real Field
1340 if self
._rational
_algebra
is None:
1341 # There's no need to construct *another* algebra over the
1342 # rationals if this one is already over the
1343 # rationals. Likewise, if we never orthonormalized our
1344 # basis, we might as well just use the given one.
1345 return super()._charpoly
_coefficients
()
1347 # Do the computation over the rationals. The answer will be
1348 # the same, because all we've done is a change of basis.
1349 # Then, change back from QQ to our real base ring
1350 a
= ( a_i
.change_ring(self
.base_ring())
1351 for a_i
in self
._rational
_algebra
._charpoly
_coefficients
() )
1353 if self
._deortho
_matrix
is None:
1354 # This can happen if our base ring was, say, AA and we
1355 # chose not to (or didn't need to) orthonormalize. It's
1356 # still faster to do the computations over QQ even if
1357 # the numbers in the boxes stay the same.
1360 # Otherwise, convert the coordinate variables back to the
1361 # deorthonormalized ones.
1362 R
= self
.coordinate_polynomial_ring()
1363 from sage
.modules
.free_module_element
import vector
1364 X
= vector(R
, R
.gens())
1365 BX
= self
._deortho
_matrix
*X
1367 subs_dict
= { X[i]: BX[i] for i in range(len(X)) }
1368 return tuple( a_i
.subs(subs_dict
) for a_i
in a
)
1370 class ConcreteEJA(RationalBasisEJA
):
1372 A class for the Euclidean Jordan algebras that we know by name.
1374 These are the Jordan algebras whose basis, multiplication table,
1375 rank, and so on are known a priori. More to the point, they are
1376 the Euclidean Jordan algebras for which we are able to conjure up
1377 a "random instance."
1381 sage: from mjo.eja.eja_algebra import ConcreteEJA
1385 Our basis is normalized with respect to the algebra's inner
1386 product, unless we specify otherwise::
1388 sage: set_random_seed()
1389 sage: J = ConcreteEJA.random_instance()
1390 sage: all( b.norm() == 1 for b in J.gens() )
1393 Since our basis is orthonormal with respect to the algebra's inner
1394 product, and since we know that this algebra is an EJA, any
1395 left-multiplication operator's matrix will be symmetric because
1396 natural->EJA basis representation is an isometry and within the
1397 EJA the operator is self-adjoint by the Jordan axiom::
1399 sage: set_random_seed()
1400 sage: J = ConcreteEJA.random_instance()
1401 sage: x = J.random_element()
1402 sage: x.operator().is_self_adjoint()
1407 def _max_random_instance_size():
1409 Return an integer "size" that is an upper bound on the size of
1410 this algebra when it is used in a random test
1411 case. Unfortunately, the term "size" is ambiguous -- when
1412 dealing with `R^n` under either the Hadamard or Jordan spin
1413 product, the "size" refers to the dimension `n`. When dealing
1414 with a matrix algebra (real symmetric or complex/quaternion
1415 Hermitian), it refers to the size of the matrix, which is far
1416 less than the dimension of the underlying vector space.
1418 This method must be implemented in each subclass.
1420 raise NotImplementedError
1423 def random_instance(cls
, *args
, **kwargs
):
1425 Return a random instance of this type of algebra.
1427 This method should be implemented in each subclass.
1429 from sage
.misc
.prandom
import choice
1430 eja_class
= choice(cls
.__subclasses
__())
1432 # These all bubble up to the RationalBasisEJA superclass
1433 # constructor, so any (kw)args valid there are also valid
1435 return eja_class
.random_instance(*args
, **kwargs
)
1440 def dimension_over_reals():
1442 The dimension of this matrix's base ring over the reals.
1444 The reals are dimension one over themselves, obviously; that's
1445 just `\mathbb{R}^{1}`. Likewise, the complex numbers `a + bi`
1446 have dimension two. Finally, the quaternions have dimension
1447 four over the reals.
1449 This is used to determine the size of the matrix returned from
1450 :meth:`real_embed`, among other things.
1452 raise NotImplementedError
1455 def real_embed(cls
,M
):
1457 Embed the matrix ``M`` into a space of real matrices.
1459 The matrix ``M`` can have entries in any field at the moment:
1460 the real numbers, complex numbers, or quaternions. And although
1461 they are not a field, we can probably support octonions at some
1462 point, too. This function returns a real matrix that "acts like"
1463 the original with respect to matrix multiplication; i.e.
1465 real_embed(M*N) = real_embed(M)*real_embed(N)
1468 if M
.ncols() != M
.nrows():
1469 raise ValueError("the matrix 'M' must be square")
1474 def real_unembed(cls
,M
):
1476 The inverse of :meth:`real_embed`.
1478 if M
.ncols() != M
.nrows():
1479 raise ValueError("the matrix 'M' must be square")
1480 if not ZZ(M
.nrows()).mod(cls
.dimension_over_reals()).is_zero():
1481 raise ValueError("the matrix 'M' must be a real embedding")
1485 def jordan_product(X
,Y
):
1486 return (X
*Y
+ Y
*X
)/2
1489 def trace_inner_product(cls
,X
,Y
):
1491 Compute the trace inner-product of two real-embeddings.
1495 sage: from mjo.eja.eja_algebra import (RealSymmetricEJA,
1496 ....: ComplexHermitianEJA,
1497 ....: QuaternionHermitianEJA)
1501 This gives the same answer as it would if we computed the trace
1502 from the unembedded (original) matrices::
1504 sage: set_random_seed()
1505 sage: J = RealSymmetricEJA.random_instance()
1506 sage: x,y = J.random_elements(2)
1507 sage: Xe = x.to_matrix()
1508 sage: Ye = y.to_matrix()
1509 sage: X = J.real_unembed(Xe)
1510 sage: Y = J.real_unembed(Ye)
1511 sage: expected = (X*Y).trace()
1512 sage: actual = J.trace_inner_product(Xe,Ye)
1513 sage: actual == expected
1518 sage: set_random_seed()
1519 sage: J = ComplexHermitianEJA.random_instance()
1520 sage: x,y = J.random_elements(2)
1521 sage: Xe = x.to_matrix()
1522 sage: Ye = y.to_matrix()
1523 sage: X = J.real_unembed(Xe)
1524 sage: Y = J.real_unembed(Ye)
1525 sage: expected = (X*Y).trace().real()
1526 sage: actual = J.trace_inner_product(Xe,Ye)
1527 sage: actual == expected
1532 sage: set_random_seed()
1533 sage: J = QuaternionHermitianEJA.random_instance()
1534 sage: x,y = J.random_elements(2)
1535 sage: Xe = x.to_matrix()
1536 sage: Ye = y.to_matrix()
1537 sage: X = J.real_unembed(Xe)
1538 sage: Y = J.real_unembed(Ye)
1539 sage: expected = (X*Y).trace().coefficient_tuple()[0]
1540 sage: actual = J.trace_inner_product(Xe,Ye)
1541 sage: actual == expected
1545 Xu
= cls
.real_unembed(X
)
1546 Yu
= cls
.real_unembed(Y
)
1547 tr
= (Xu
*Yu
).trace()
1550 # Works in QQ, AA, RDF, et cetera.
1552 except AttributeError:
1553 # A quaternion doesn't have a real() method, but does
1554 # have coefficient_tuple() method that returns the
1555 # coefficients of 1, i, j, and k -- in that order.
1556 return tr
.coefficient_tuple()[0]
1559 class RealMatrixEJA(MatrixEJA
):
1561 def dimension_over_reals():
1565 class RealSymmetricEJA(ConcreteEJA
, RealMatrixEJA
):
1567 The rank-n simple EJA consisting of real symmetric n-by-n
1568 matrices, the usual symmetric Jordan product, and the trace inner
1569 product. It has dimension `(n^2 + n)/2` over the reals.
1573 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1577 sage: J = RealSymmetricEJA(2)
1578 sage: e0, e1, e2 = J.gens()
1586 In theory, our "field" can be any subfield of the reals::
1588 sage: RealSymmetricEJA(2, field=RDF)
1589 Euclidean Jordan algebra of dimension 3 over Real Double Field
1590 sage: RealSymmetricEJA(2, field=RR)
1591 Euclidean Jordan algebra of dimension 3 over Real Field with
1592 53 bits of precision
1596 The dimension of this algebra is `(n^2 + n) / 2`::
1598 sage: set_random_seed()
1599 sage: n_max = RealSymmetricEJA._max_random_instance_size()
1600 sage: n = ZZ.random_element(1, n_max)
1601 sage: J = RealSymmetricEJA(n)
1602 sage: J.dimension() == (n^2 + n)/2
1605 The Jordan multiplication is what we think it is::
1607 sage: set_random_seed()
1608 sage: J = RealSymmetricEJA.random_instance()
1609 sage: x,y = J.random_elements(2)
1610 sage: actual = (x*y).to_matrix()
1611 sage: X = x.to_matrix()
1612 sage: Y = y.to_matrix()
1613 sage: expected = (X*Y + Y*X)/2
1614 sage: actual == expected
1616 sage: J(expected) == x*y
1619 We can change the generator prefix::
1621 sage: RealSymmetricEJA(3, prefix='q').gens()
1622 (q0, q1, q2, q3, q4, q5)
1624 We can construct the (trivial) algebra of rank zero::
1626 sage: RealSymmetricEJA(0)
1627 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1631 def _denormalized_basis(cls
, n
):
1633 Return a basis for the space of real symmetric n-by-n matrices.
1637 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1641 sage: set_random_seed()
1642 sage: n = ZZ.random_element(1,5)
1643 sage: B = RealSymmetricEJA._denormalized_basis(n)
1644 sage: all( M.is_symmetric() for M in B)
1648 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1652 for j
in range(i
+1):
1653 Eij
= matrix(ZZ
, n
, lambda k
,l
: k
==i
and l
==j
)
1657 Sij
= Eij
+ Eij
.transpose()
1663 def _max_random_instance_size():
1664 return 4 # Dimension 10
1667 def random_instance(cls
, **kwargs
):
1669 Return a random instance of this type of algebra.
1671 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1672 return cls(n
, **kwargs
)
1674 def __init__(self
, n
, **kwargs
):
1675 # We know this is a valid EJA, but will double-check
1676 # if the user passes check_axioms=True.
1677 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
1679 super(RealSymmetricEJA
, self
).__init
__(self
._denormalized
_basis
(n
),
1680 self
.jordan_product
,
1681 self
.trace_inner_product
,
1684 # TODO: this could be factored out somehow, but is left here
1685 # because the MatrixEJA is not presently a subclass of the
1686 # FDEJA class that defines rank() and one().
1687 self
.rank
.set_cache(n
)
1688 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
1689 self
.one
.set_cache(self(idV
))
1693 class ComplexMatrixEJA(MatrixEJA
):
1694 # A manual dictionary-cache for the complex_extension() method,
1695 # since apparently @classmethods can't also be @cached_methods.
1696 _complex_extension
= {}
1699 def complex_extension(cls
,field
):
1701 The complex field that we embed/unembed, as an extension
1702 of the given ``field``.
1704 if field
in cls
._complex
_extension
:
1705 return cls
._complex
_extension
[field
]
1707 # Sage doesn't know how to adjoin the complex "i" (the root of
1708 # x^2 + 1) to a field in a general way. Here, we just enumerate
1709 # all of the cases that I have cared to support so far.
1711 # Sage doesn't know how to embed AA into QQbar, i.e. how
1712 # to adjoin sqrt(-1) to AA.
1714 elif not field
.is_exact():
1716 F
= field
.complex_field()
1718 # Works for QQ and... maybe some other fields.
1719 R
= PolynomialRing(field
, 'z')
1721 F
= field
.extension(z
**2 + 1, 'I', embedding
=CLF(-1).sqrt())
1723 cls
._complex
_extension
[field
] = F
1727 def dimension_over_reals():
1731 def real_embed(cls
,M
):
1733 Embed the n-by-n complex matrix ``M`` into the space of real
1734 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1735 bi` to the block matrix ``[[a,b],[-b,a]]``.
1739 sage: from mjo.eja.eja_algebra import ComplexMatrixEJA
1743 sage: F = QuadraticField(-1, 'I')
1744 sage: x1 = F(4 - 2*i)
1745 sage: x2 = F(1 + 2*i)
1748 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1749 sage: ComplexMatrixEJA.real_embed(M)
1758 Embedding is a homomorphism (isomorphism, in fact)::
1760 sage: set_random_seed()
1761 sage: n = ZZ.random_element(3)
1762 sage: F = QuadraticField(-1, 'I')
1763 sage: X = random_matrix(F, n)
1764 sage: Y = random_matrix(F, n)
1765 sage: Xe = ComplexMatrixEJA.real_embed(X)
1766 sage: Ye = ComplexMatrixEJA.real_embed(Y)
1767 sage: XYe = ComplexMatrixEJA.real_embed(X*Y)
1772 super(ComplexMatrixEJA
,cls
).real_embed(M
)
1775 # We don't need any adjoined elements...
1776 field
= M
.base_ring().base_ring()
1782 blocks
.append(matrix(field
, 2, [ [ a
, b
],
1785 return matrix
.block(field
, n
, blocks
)
1789 def real_unembed(cls
,M
):
1791 The inverse of _embed_complex_matrix().
1795 sage: from mjo.eja.eja_algebra import ComplexMatrixEJA
1799 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1800 ....: [-2, 1, -4, 3],
1801 ....: [ 9, 10, 11, 12],
1802 ....: [-10, 9, -12, 11] ])
1803 sage: ComplexMatrixEJA.real_unembed(A)
1805 [ 10*I + 9 12*I + 11]
1809 Unembedding is the inverse of embedding::
1811 sage: set_random_seed()
1812 sage: F = QuadraticField(-1, 'I')
1813 sage: M = random_matrix(F, 3)
1814 sage: Me = ComplexMatrixEJA.real_embed(M)
1815 sage: ComplexMatrixEJA.real_unembed(Me) == M
1819 super(ComplexMatrixEJA
,cls
).real_unembed(M
)
1821 d
= cls
.dimension_over_reals()
1822 F
= cls
.complex_extension(M
.base_ring())
1825 # Go top-left to bottom-right (reading order), converting every
1826 # 2-by-2 block we see to a single complex element.
1828 for k
in range(n
/d
):
1829 for j
in range(n
/d
):
1830 submat
= M
[d
*k
:d
*k
+d
,d
*j
:d
*j
+d
]
1831 if submat
[0,0] != submat
[1,1]:
1832 raise ValueError('bad on-diagonal submatrix')
1833 if submat
[0,1] != -submat
[1,0]:
1834 raise ValueError('bad off-diagonal submatrix')
1835 z
= submat
[0,0] + submat
[0,1]*i
1838 return matrix(F
, n
/d
, elements
)
1841 class ComplexHermitianEJA(ConcreteEJA
, ComplexMatrixEJA
):
1843 The rank-n simple EJA consisting of complex Hermitian n-by-n
1844 matrices over the real numbers, the usual symmetric Jordan product,
1845 and the real-part-of-trace inner product. It has dimension `n^2` over
1850 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1854 In theory, our "field" can be any subfield of the reals::
1856 sage: ComplexHermitianEJA(2, field=RDF)
1857 Euclidean Jordan algebra of dimension 4 over Real Double Field
1858 sage: ComplexHermitianEJA(2, field=RR)
1859 Euclidean Jordan algebra of dimension 4 over Real Field with
1860 53 bits of precision
1864 The dimension of this algebra is `n^2`::
1866 sage: set_random_seed()
1867 sage: n_max = ComplexHermitianEJA._max_random_instance_size()
1868 sage: n = ZZ.random_element(1, n_max)
1869 sage: J = ComplexHermitianEJA(n)
1870 sage: J.dimension() == n^2
1873 The Jordan multiplication is what we think it is::
1875 sage: set_random_seed()
1876 sage: J = ComplexHermitianEJA.random_instance()
1877 sage: x,y = J.random_elements(2)
1878 sage: actual = (x*y).to_matrix()
1879 sage: X = x.to_matrix()
1880 sage: Y = y.to_matrix()
1881 sage: expected = (X*Y + Y*X)/2
1882 sage: actual == expected
1884 sage: J(expected) == x*y
1887 We can change the generator prefix::
1889 sage: ComplexHermitianEJA(2, prefix='z').gens()
1892 We can construct the (trivial) algebra of rank zero::
1894 sage: ComplexHermitianEJA(0)
1895 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1900 def _denormalized_basis(cls
, n
):
1902 Returns a basis for the space of complex Hermitian n-by-n matrices.
1904 Why do we embed these? Basically, because all of numerical linear
1905 algebra assumes that you're working with vectors consisting of `n`
1906 entries from a field and scalars from the same field. There's no way
1907 to tell SageMath that (for example) the vectors contain complex
1908 numbers, while the scalar field is real.
1912 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1916 sage: set_random_seed()
1917 sage: n = ZZ.random_element(1,5)
1918 sage: B = ComplexHermitianEJA._denormalized_basis(n)
1919 sage: all( M.is_symmetric() for M in B)
1924 R
= PolynomialRing(field
, 'z')
1926 F
= field
.extension(z
**2 + 1, 'I')
1929 # This is like the symmetric case, but we need to be careful:
1931 # * We want conjugate-symmetry, not just symmetry.
1932 # * The diagonal will (as a result) be real.
1935 Eij
= matrix
.zero(F
,n
)
1937 for j
in range(i
+1):
1941 Sij
= cls
.real_embed(Eij
)
1944 # The second one has a minus because it's conjugated.
1945 Eij
[j
,i
] = 1 # Eij = Eij + Eij.transpose()
1946 Sij_real
= cls
.real_embed(Eij
)
1948 # Eij = I*Eij - I*Eij.transpose()
1951 Sij_imag
= cls
.real_embed(Eij
)
1957 # Since we embedded these, we can drop back to the "field" that we
1958 # started with instead of the complex extension "F".
1959 return tuple( s
.change_ring(field
) for s
in S
)
1962 def __init__(self
, n
, **kwargs
):
1963 # We know this is a valid EJA, but will double-check
1964 # if the user passes check_axioms=True.
1965 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
1967 super(ComplexHermitianEJA
, self
).__init
__(self
._denormalized
_basis
(n
),
1968 self
.jordan_product
,
1969 self
.trace_inner_product
,
1971 # TODO: this could be factored out somehow, but is left here
1972 # because the MatrixEJA is not presently a subclass of the
1973 # FDEJA class that defines rank() and one().
1974 self
.rank
.set_cache(n
)
1975 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
1976 self
.one
.set_cache(self(idV
))
1979 def _max_random_instance_size():
1980 return 3 # Dimension 9
1983 def random_instance(cls
, **kwargs
):
1985 Return a random instance of this type of algebra.
1987 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1988 return cls(n
, **kwargs
)
1990 class QuaternionMatrixEJA(MatrixEJA
):
1992 # A manual dictionary-cache for the quaternion_extension() method,
1993 # since apparently @classmethods can't also be @cached_methods.
1994 _quaternion_extension
= {}
1997 def quaternion_extension(cls
,field
):
1999 The quaternion field that we embed/unembed, as an extension
2000 of the given ``field``.
2002 if field
in cls
._quaternion
_extension
:
2003 return cls
._quaternion
_extension
[field
]
2005 Q
= QuaternionAlgebra(field
,-1,-1)
2007 cls
._quaternion
_extension
[field
] = Q
2011 def dimension_over_reals():
2015 def real_embed(cls
,M
):
2017 Embed the n-by-n quaternion matrix ``M`` into the space of real
2018 matrices of size 4n-by-4n by first sending each quaternion entry `z
2019 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
2020 c+di],[-c + di, a-bi]]`, and then embedding those into a real
2025 sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA
2029 sage: Q = QuaternionAlgebra(QQ,-1,-1)
2030 sage: i,j,k = Q.gens()
2031 sage: x = 1 + 2*i + 3*j + 4*k
2032 sage: M = matrix(Q, 1, [[x]])
2033 sage: QuaternionMatrixEJA.real_embed(M)
2039 Embedding is a homomorphism (isomorphism, in fact)::
2041 sage: set_random_seed()
2042 sage: n = ZZ.random_element(2)
2043 sage: Q = QuaternionAlgebra(QQ,-1,-1)
2044 sage: X = random_matrix(Q, n)
2045 sage: Y = random_matrix(Q, n)
2046 sage: Xe = QuaternionMatrixEJA.real_embed(X)
2047 sage: Ye = QuaternionMatrixEJA.real_embed(Y)
2048 sage: XYe = QuaternionMatrixEJA.real_embed(X*Y)
2053 super(QuaternionMatrixEJA
,cls
).real_embed(M
)
2054 quaternions
= M
.base_ring()
2057 F
= QuadraticField(-1, 'I')
2062 t
= z
.coefficient_tuple()
2067 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
2068 [-c
+ d
*i
, a
- b
*i
]])
2069 realM
= ComplexMatrixEJA
.real_embed(cplxM
)
2070 blocks
.append(realM
)
2072 # We should have real entries by now, so use the realest field
2073 # we've got for the return value.
2074 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
2079 def real_unembed(cls
,M
):
2081 The inverse of _embed_quaternion_matrix().
2085 sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA
2089 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
2090 ....: [-2, 1, -4, 3],
2091 ....: [-3, 4, 1, -2],
2092 ....: [-4, -3, 2, 1]])
2093 sage: QuaternionMatrixEJA.real_unembed(M)
2094 [1 + 2*i + 3*j + 4*k]
2098 Unembedding is the inverse of embedding::
2100 sage: set_random_seed()
2101 sage: Q = QuaternionAlgebra(QQ, -1, -1)
2102 sage: M = random_matrix(Q, 3)
2103 sage: Me = QuaternionMatrixEJA.real_embed(M)
2104 sage: QuaternionMatrixEJA.real_unembed(Me) == M
2108 super(QuaternionMatrixEJA
,cls
).real_unembed(M
)
2110 d
= cls
.dimension_over_reals()
2112 # Use the base ring of the matrix to ensure that its entries can be
2113 # multiplied by elements of the quaternion algebra.
2114 Q
= cls
.quaternion_extension(M
.base_ring())
2117 # Go top-left to bottom-right (reading order), converting every
2118 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
2121 for l
in range(n
/d
):
2122 for m
in range(n
/d
):
2123 submat
= ComplexMatrixEJA
.real_unembed(
2124 M
[d
*l
:d
*l
+d
,d
*m
:d
*m
+d
] )
2125 if submat
[0,0] != submat
[1,1].conjugate():
2126 raise ValueError('bad on-diagonal submatrix')
2127 if submat
[0,1] != -submat
[1,0].conjugate():
2128 raise ValueError('bad off-diagonal submatrix')
2129 z
= submat
[0,0].real()
2130 z
+= submat
[0,0].imag()*i
2131 z
+= submat
[0,1].real()*j
2132 z
+= submat
[0,1].imag()*k
2135 return matrix(Q
, n
/d
, elements
)
2138 class QuaternionHermitianEJA(ConcreteEJA
, QuaternionMatrixEJA
):
2140 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
2141 matrices, the usual symmetric Jordan product, and the
2142 real-part-of-trace inner product. It has dimension `2n^2 - n` over
2147 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2151 In theory, our "field" can be any subfield of the reals::
2153 sage: QuaternionHermitianEJA(2, field=RDF)
2154 Euclidean Jordan algebra of dimension 6 over Real Double Field
2155 sage: QuaternionHermitianEJA(2, field=RR)
2156 Euclidean Jordan algebra of dimension 6 over Real Field with
2157 53 bits of precision
2161 The dimension of this algebra is `2*n^2 - n`::
2163 sage: set_random_seed()
2164 sage: n_max = QuaternionHermitianEJA._max_random_instance_size()
2165 sage: n = ZZ.random_element(1, n_max)
2166 sage: J = QuaternionHermitianEJA(n)
2167 sage: J.dimension() == 2*(n^2) - n
2170 The Jordan multiplication is what we think it is::
2172 sage: set_random_seed()
2173 sage: J = QuaternionHermitianEJA.random_instance()
2174 sage: x,y = J.random_elements(2)
2175 sage: actual = (x*y).to_matrix()
2176 sage: X = x.to_matrix()
2177 sage: Y = y.to_matrix()
2178 sage: expected = (X*Y + Y*X)/2
2179 sage: actual == expected
2181 sage: J(expected) == x*y
2184 We can change the generator prefix::
2186 sage: QuaternionHermitianEJA(2, prefix='a').gens()
2187 (a0, a1, a2, a3, a4, a5)
2189 We can construct the (trivial) algebra of rank zero::
2191 sage: QuaternionHermitianEJA(0)
2192 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2196 def _denormalized_basis(cls
, n
):
2198 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
2200 Why do we embed these? Basically, because all of numerical
2201 linear algebra assumes that you're working with vectors consisting
2202 of `n` entries from a field and scalars from the same field. There's
2203 no way to tell SageMath that (for example) the vectors contain
2204 complex numbers, while the scalar field is real.
2208 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2212 sage: set_random_seed()
2213 sage: n = ZZ.random_element(1,5)
2214 sage: B = QuaternionHermitianEJA._denormalized_basis(n)
2215 sage: all( M.is_symmetric() for M in B )
2220 Q
= QuaternionAlgebra(QQ
,-1,-1)
2223 # This is like the symmetric case, but we need to be careful:
2225 # * We want conjugate-symmetry, not just symmetry.
2226 # * The diagonal will (as a result) be real.
2229 Eij
= matrix
.zero(Q
,n
)
2231 for j
in range(i
+1):
2235 Sij
= cls
.real_embed(Eij
)
2238 # The second, third, and fourth ones have a minus
2239 # because they're conjugated.
2240 # Eij = Eij + Eij.transpose()
2242 Sij_real
= cls
.real_embed(Eij
)
2244 # Eij = I*(Eij - Eij.transpose())
2247 Sij_I
= cls
.real_embed(Eij
)
2249 # Eij = J*(Eij - Eij.transpose())
2252 Sij_J
= cls
.real_embed(Eij
)
2254 # Eij = K*(Eij - Eij.transpose())
2257 Sij_K
= cls
.real_embed(Eij
)
2263 # Since we embedded these, we can drop back to the "field" that we
2264 # started with instead of the quaternion algebra "Q".
2265 return tuple( s
.change_ring(field
) for s
in S
)
2268 def __init__(self
, n
, **kwargs
):
2269 # We know this is a valid EJA, but will double-check
2270 # if the user passes check_axioms=True.
2271 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2273 super(QuaternionHermitianEJA
, self
).__init
__(self
._denormalized
_basis
(n
),
2274 self
.jordan_product
,
2275 self
.trace_inner_product
,
2277 # TODO: this could be factored out somehow, but is left here
2278 # because the MatrixEJA is not presently a subclass of the
2279 # FDEJA class that defines rank() and one().
2280 self
.rank
.set_cache(n
)
2281 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
2282 self
.one
.set_cache(self(idV
))
2286 def _max_random_instance_size():
2288 The maximum rank of a random QuaternionHermitianEJA.
2290 return 2 # Dimension 6
2293 def random_instance(cls
, **kwargs
):
2295 Return a random instance of this type of algebra.
2297 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2298 return cls(n
, **kwargs
)
2301 class HadamardEJA(ConcreteEJA
):
2303 Return the Euclidean Jordan Algebra corresponding to the set
2304 `R^n` under the Hadamard product.
2306 Note: this is nothing more than the Cartesian product of ``n``
2307 copies of the spin algebra. Once Cartesian product algebras
2308 are implemented, this can go.
2312 sage: from mjo.eja.eja_algebra import HadamardEJA
2316 This multiplication table can be verified by hand::
2318 sage: J = HadamardEJA(3)
2319 sage: e0,e1,e2 = J.gens()
2335 We can change the generator prefix::
2337 sage: HadamardEJA(3, prefix='r').gens()
2341 def __init__(self
, n
, **kwargs
):
2343 jordan_product
= lambda x
,y
: x
2344 inner_product
= lambda x
,y
: x
2346 def jordan_product(x
,y
):
2348 return P( xi
*yi
for (xi
,yi
) in zip(x
,y
) )
2350 def inner_product(x
,y
):
2353 # New defaults for keyword arguments. Don't orthonormalize
2354 # because our basis is already orthonormal with respect to our
2355 # inner-product. Don't check the axioms, because we know this
2356 # is a valid EJA... but do double-check if the user passes
2357 # check_axioms=True. Note: we DON'T override the "check_field"
2358 # default here, because the user can pass in a field!
2359 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2360 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2362 column_basis
= tuple( b
.column() for b
in FreeModule(ZZ
, n
).basis() )
2363 super().__init
__(column_basis
, jordan_product
, inner_product
, **kwargs
)
2364 self
.rank
.set_cache(n
)
2367 self
.one
.set_cache( self
.zero() )
2369 self
.one
.set_cache( sum(self
.gens()) )
2372 def _max_random_instance_size():
2374 The maximum dimension of a random HadamardEJA.
2379 def random_instance(cls
, **kwargs
):
2381 Return a random instance of this type of algebra.
2383 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2384 return cls(n
, **kwargs
)
2387 class BilinearFormEJA(ConcreteEJA
):
2389 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2390 with the half-trace inner product and jordan product ``x*y =
2391 (<Bx,y>,y_bar>, x0*y_bar + y0*x_bar)`` where `B = 1 \times B22` is
2392 a symmetric positive-definite "bilinear form" matrix. Its
2393 dimension is the size of `B`, and it has rank two in dimensions
2394 larger than two. It reduces to the ``JordanSpinEJA`` when `B` is
2395 the identity matrix of order ``n``.
2397 We insist that the one-by-one upper-left identity block of `B` be
2398 passed in as well so that we can be passed a matrix of size zero
2399 to construct a trivial algebra.
2403 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
2404 ....: JordanSpinEJA)
2408 When no bilinear form is specified, the identity matrix is used,
2409 and the resulting algebra is the Jordan spin algebra::
2411 sage: B = matrix.identity(AA,3)
2412 sage: J0 = BilinearFormEJA(B)
2413 sage: J1 = JordanSpinEJA(3)
2414 sage: J0.multiplication_table() == J0.multiplication_table()
2417 An error is raised if the matrix `B` does not correspond to a
2418 positive-definite bilinear form::
2420 sage: B = matrix.random(QQ,2,3)
2421 sage: J = BilinearFormEJA(B)
2422 Traceback (most recent call last):
2424 ValueError: bilinear form is not positive-definite
2425 sage: B = matrix.zero(QQ,3)
2426 sage: J = BilinearFormEJA(B)
2427 Traceback (most recent call last):
2429 ValueError: bilinear form is not positive-definite
2433 We can create a zero-dimensional algebra::
2435 sage: B = matrix.identity(AA,0)
2436 sage: J = BilinearFormEJA(B)
2440 We can check the multiplication condition given in the Jordan, von
2441 Neumann, and Wigner paper (and also discussed on my "On the
2442 symmetry..." paper). Note that this relies heavily on the standard
2443 choice of basis, as does anything utilizing the bilinear form
2444 matrix. We opt not to orthonormalize the basis, because if we
2445 did, we would have to normalize the `s_{i}` in a similar manner::
2447 sage: set_random_seed()
2448 sage: n = ZZ.random_element(5)
2449 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2450 sage: B11 = matrix.identity(QQ,1)
2451 sage: B22 = M.transpose()*M
2452 sage: B = block_matrix(2,2,[ [B11,0 ],
2454 sage: J = BilinearFormEJA(B, orthonormalize=False)
2455 sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
2456 sage: V = J.vector_space()
2457 sage: sis = [ J( V([0] + (M.inverse()*ei).list()).column() )
2458 ....: for ei in eis ]
2459 sage: actual = [ sis[i]*sis[j]
2460 ....: for i in range(n-1)
2461 ....: for j in range(n-1) ]
2462 sage: expected = [ J.one() if i == j else J.zero()
2463 ....: for i in range(n-1)
2464 ....: for j in range(n-1) ]
2465 sage: actual == expected
2469 def __init__(self
, B
, **kwargs
):
2470 # The matrix "B" is supplied by the user in most cases,
2471 # so it makes sense to check whether or not its positive-
2472 # definite unless we are specifically asked not to...
2473 if ("check_axioms" not in kwargs
) or kwargs
["check_axioms"]:
2474 if not B
.is_positive_definite():
2475 raise ValueError("bilinear form is not positive-definite")
2477 # However, all of the other data for this EJA is computed
2478 # by us in manner that guarantees the axioms are
2479 # satisfied. So, again, unless we are specifically asked to
2480 # verify things, we'll skip the rest of the checks.
2481 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2483 def inner_product(x
,y
):
2484 return (y
.T
*B
*x
)[0,0]
2486 def jordan_product(x
,y
):
2492 z0
= inner_product(y
,x
)
2493 zbar
= y0
*xbar
+ x0
*ybar
2494 return P([z0
] + zbar
.list())
2497 column_basis
= tuple( b
.column() for b
in FreeModule(ZZ
, n
).basis() )
2498 super(BilinearFormEJA
, self
).__init
__(column_basis
,
2503 # The rank of this algebra is two, unless we're in a
2504 # one-dimensional ambient space (because the rank is bounded
2505 # by the ambient dimension).
2506 self
.rank
.set_cache(min(n
,2))
2509 self
.one
.set_cache( self
.zero() )
2511 self
.one
.set_cache( self
.monomial(0) )
2514 def _max_random_instance_size():
2516 The maximum dimension of a random BilinearFormEJA.
2521 def random_instance(cls
, **kwargs
):
2523 Return a random instance of this algebra.
2525 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2527 B
= matrix
.identity(ZZ
, n
)
2528 return cls(B
, **kwargs
)
2530 B11
= matrix
.identity(ZZ
, 1)
2531 M
= matrix
.random(ZZ
, n
-1)
2532 I
= matrix
.identity(ZZ
, n
-1)
2534 while alpha
.is_zero():
2535 alpha
= ZZ
.random_element().abs()
2536 B22
= M
.transpose()*M
+ alpha
*I
2538 from sage
.matrix
.special
import block_matrix
2539 B
= block_matrix(2,2, [ [B11
, ZZ(0) ],
2542 return cls(B
, **kwargs
)
2545 class JordanSpinEJA(BilinearFormEJA
):
2547 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2548 with the usual inner product and jordan product ``x*y =
2549 (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2554 sage: from mjo.eja.eja_algebra import JordanSpinEJA
2558 This multiplication table can be verified by hand::
2560 sage: J = JordanSpinEJA(4)
2561 sage: e0,e1,e2,e3 = J.gens()
2577 We can change the generator prefix::
2579 sage: JordanSpinEJA(2, prefix='B').gens()
2584 Ensure that we have the usual inner product on `R^n`::
2586 sage: set_random_seed()
2587 sage: J = JordanSpinEJA.random_instance()
2588 sage: x,y = J.random_elements(2)
2589 sage: actual = x.inner_product(y)
2590 sage: expected = x.to_vector().inner_product(y.to_vector())
2591 sage: actual == expected
2595 def __init__(self
, n
, **kwargs
):
2596 # This is a special case of the BilinearFormEJA with the
2597 # identity matrix as its bilinear form.
2598 B
= matrix
.identity(ZZ
, n
)
2600 # Don't orthonormalize because our basis is already
2601 # orthonormal with respect to our inner-product.
2602 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2604 # But also don't pass check_field=False here, because the user
2605 # can pass in a field!
2606 super(JordanSpinEJA
, self
).__init
__(B
, **kwargs
)
2609 def _max_random_instance_size():
2611 The maximum dimension of a random JordanSpinEJA.
2616 def random_instance(cls
, **kwargs
):
2618 Return a random instance of this type of algebra.
2620 Needed here to override the implementation for ``BilinearFormEJA``.
2622 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2623 return cls(n
, **kwargs
)
2626 class TrivialEJA(ConcreteEJA
):
2628 The trivial Euclidean Jordan algebra consisting of only a zero element.
2632 sage: from mjo.eja.eja_algebra import TrivialEJA
2636 sage: J = TrivialEJA()
2643 sage: 7*J.one()*12*J.one()
2645 sage: J.one().inner_product(J.one())
2647 sage: J.one().norm()
2649 sage: J.one().subalgebra_generated_by()
2650 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2655 def __init__(self
, **kwargs
):
2656 jordan_product
= lambda x
,y
: x
2657 inner_product
= lambda x
,y
: 0
2660 # New defaults for keyword arguments
2661 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2662 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2664 super(TrivialEJA
, self
).__init
__(basis
,
2668 # The rank is zero using my definition, namely the dimension of the
2669 # largest subalgebra generated by any element.
2670 self
.rank
.set_cache(0)
2671 self
.one
.set_cache( self
.zero() )
2674 def random_instance(cls
, **kwargs
):
2675 # We don't take a "size" argument so the superclass method is
2676 # inappropriate for us.
2677 return cls(**kwargs
)
2680 class CartesianProductEJA(CombinatorialFreeModule_CartesianProduct
,
2681 FiniteDimensionalEJA
):
2683 The external (orthogonal) direct sum of two or more Euclidean
2684 Jordan algebras. Every Euclidean Jordan algebra decomposes into an
2685 orthogonal direct sum of simple Euclidean Jordan algebras which is
2686 then isometric to a Cartesian product, so no generality is lost by
2687 providing only this construction.
2691 sage: from mjo.eja.eja_algebra import (CartesianProductEJA,
2693 ....: JordanSpinEJA,
2694 ....: RealSymmetricEJA)
2698 The Jordan product is inherited from our factors and implemented by
2699 our CombinatorialFreeModule Cartesian product superclass::
2701 sage: set_random_seed()
2702 sage: J1 = HadamardEJA(2)
2703 sage: J2 = RealSymmetricEJA(2)
2704 sage: J = cartesian_product([J1,J2])
2705 sage: x,y = J.random_elements(2)
2709 The ability to retrieve the original factors is implemented by our
2710 CombinatorialFreeModule Cartesian product superclass::
2712 sage: J1 = HadamardEJA(2, field=QQ)
2713 sage: J2 = JordanSpinEJA(3, field=QQ)
2714 sage: J = cartesian_product([J1,J2])
2715 sage: J.cartesian_factors()
2716 (Euclidean Jordan algebra of dimension 2 over Rational Field,
2717 Euclidean Jordan algebra of dimension 3 over Rational Field)
2721 All factors must share the same base field::
2723 sage: J1 = HadamardEJA(2, field=QQ)
2724 sage: J2 = RealSymmetricEJA(2)
2725 sage: CartesianProductEJA((J1,J2))
2726 Traceback (most recent call last):
2728 ValueError: all factors must share the same base field
2730 def __init__(self
, modules
, **kwargs
):
2731 CombinatorialFreeModule_CartesianProduct
.__init
__(self
, modules
, **kwargs
)
2732 field
= modules
[0].base_ring()
2733 if not all( J
.base_ring() == field
for J
in modules
):
2734 raise ValueError("all factors must share the same base field")
2736 M
= cartesian_product( [J
.matrix_space() for J
in modules
] )
2739 W
= VectorSpace(field
,m
)
2740 self
._matrix
_basis
= []
2742 for a
in modules
[k
].matrix_basis():
2745 self
._matrix
_basis
.append(M(v
))
2747 self
._matrix
_basis
= tuple(self
._matrix
_basis
)
2749 n
= len(self
._matrix
_basis
)
2752 # Initialize the FDEJA class, too. Does this override the
2753 # initialization that we did for the
2754 # CombinatorialFreeModule_CartesianProduct class? If not, we
2755 # will probably have to duplicate some of the work (i.e. one
2756 # of the constructors). Since the CartesianProduct one is
2757 # smaller, that makes the most sense to copy/paste if it comes
2761 self
.rank
.set_cache(sum(J
.rank() for J
in modules
))
2764 def cartesian_projection(self
, i
):
2768 sage: from mjo.eja.eja_algebra import (random_eja,
2770 ....: RealSymmetricEJA)
2774 sage: J1 = HadamardEJA(2)
2775 sage: J2 = RealSymmetricEJA(2)
2776 sage: J = cartesian_product([J1,J2])
2777 sage: J.cartesian_projection(0)
2778 Linear operator between finite-dimensional Euclidean Jordan
2779 algebras represented by the matrix:
2782 Domain: Euclidean Jordan algebra of dimension 2 over Algebraic
2783 Real Field (+) Euclidean Jordan algebra of dimension 3 over
2784 Algebraic Real Field
2785 Codomain: Euclidean Jordan algebra of dimension 2 over Algebraic
2787 sage: J.cartesian_projection(1)
2788 Linear operator between finite-dimensional Euclidean Jordan
2789 algebras represented by the matrix:
2793 Domain: Euclidean Jordan algebra of dimension 2 over Algebraic
2794 Real Field (+) Euclidean Jordan algebra of dimension 3 over
2795 Algebraic Real Field
2796 Codomain: Euclidean Jordan algebra of dimension 3 over Algebraic
2801 The answer never changes::
2803 sage: set_random_seed()
2804 sage: J1 = random_eja()
2805 sage: J2 = random_eja()
2806 sage: J = cartesian_product([J1,J2])
2807 sage: P0 = J.cartesian_projection(0)
2808 sage: P1 = J.cartesian_projection(0)
2813 Ji
= self
.cartesian_factors()[i
]
2814 # We reimplement the CombinatorialFreeModule superclass method
2815 # because if we don't, something gets messed up with the caching
2816 # and the answer changes the second time you run it. See the TESTS.
2817 Pi
= self
._module
_morphism
(lambda j_t
: Ji
.monomial(j_t
[1])
2818 if i
== j_t
[0] else Ji
.zero(),
2820 return FiniteDimensionalEJAOperator(self
,Ji
,Pi
.matrix())
2823 def cartesian_embedding(self
, i
):
2827 sage: from mjo.eja.eja_algebra import (random_eja,
2829 ....: RealSymmetricEJA)
2833 sage: J1 = HadamardEJA(2)
2834 sage: J2 = RealSymmetricEJA(2)
2835 sage: J = cartesian_product([J1,J2])
2838 sage: J.cartesian_embedding
2840 sage: J.cartesian_embedding(0)
2841 Linear operator between finite-dimensional Euclidean Jordan
2842 algebras represented by the matrix:
2848 Domain: Euclidean Jordan algebra of dimension 2 over
2849 Algebraic Real Field
2850 Codomain: Euclidean Jordan algebra of dimension 2 over
2851 Algebraic Real Field (+) Euclidean Jordan algebra of
2852 dimension 3 over Algebraic Real Field
2853 sage: J.cartesian_embedding(1)
2854 Linear operator between finite-dimensional Euclidean Jordan
2855 algebras represented by the matrix:
2861 Domain: Euclidean Jordan algebra of dimension 3 over
2862 Algebraic Real Field
2863 Codomain: Euclidean Jordan algebra of dimension 2 over
2864 Algebraic Real Field (+) Euclidean Jordan algebra of
2865 dimension 3 over Algebraic Real Field
2869 The answer never changes::
2871 sage: set_random_seed()
2872 sage: J1 = random_eja()
2873 sage: J2 = random_eja()
2874 sage: J = cartesian_product([J1,J2])
2875 sage: E0 = J.cartesian_embedding(0)
2876 sage: E1 = J.cartesian_embedding(0)
2881 Ji
= self
.cartesian_factors()[i
]
2882 # We reimplement the CombinatorialFreeModule superclass method
2883 # because if we don't, something gets messed up with the caching
2884 # and the answer changes the second time you run it. See the TESTS.
2885 Ei
= Ji
._module
_morphism
(lambda t
: self
.monomial((i
, t
)), codomain
=self
)
2886 return FiniteDimensionalEJAOperator(Ji
,self
,Ei
.matrix())
2889 FiniteDimensionalEJA
.CartesianProduct
= CartesianProductEJA
2892 # def projections(self):
2894 # Return a pair of projections onto this algebra's factors.
2898 # sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
2899 # ....: ComplexHermitianEJA,
2900 # ....: DirectSumEJA)
2904 # sage: J1 = JordanSpinEJA(2)
2905 # sage: J2 = ComplexHermitianEJA(2)
2906 # sage: J = DirectSumEJA(J1,J2)
2907 # sage: (pi_left, pi_right) = J.projections()
2908 # sage: J.one().to_vector()
2909 # (1, 0, 1, 0, 0, 1)
2910 # sage: pi_left(J.one()).to_vector()
2912 # sage: pi_right(J.one()).to_vector()
2916 # (J1,J2) = self.factors()
2917 # m = J1.dimension()
2918 # n = J2.dimension()
2919 # V_basis = self.vector_space().basis()
2920 # # Need to specify the dimensions explicitly so that we don't
2921 # # wind up with a zero-by-zero matrix when we want e.g. a
2922 # # zero-by-two matrix (important for composing things).
2923 # P1 = matrix(self.base_ring(), m, m+n, V_basis[:m])
2924 # P2 = matrix(self.base_ring(), n, m+n, V_basis[m:])
2925 # pi_left = FiniteDimensionalEJAOperator(self,J1,P1)
2926 # pi_right = FiniteDimensionalEJAOperator(self,J2,P2)
2927 # return (pi_left, pi_right)
2929 # def inclusions(self):
2931 # Return the pair of inclusion maps from our factors into us.
2935 # sage: from mjo.eja.eja_algebra import (random_eja,
2936 # ....: JordanSpinEJA,
2937 # ....: RealSymmetricEJA,
2938 # ....: DirectSumEJA)
2942 # sage: J1 = JordanSpinEJA(3)
2943 # sage: J2 = RealSymmetricEJA(2)
2944 # sage: J = DirectSumEJA(J1,J2)
2945 # sage: (iota_left, iota_right) = J.inclusions()
2946 # sage: iota_left(J1.zero()) == J.zero()
2948 # sage: iota_right(J2.zero()) == J.zero()
2950 # sage: J1.one().to_vector()
2952 # sage: iota_left(J1.one()).to_vector()
2953 # (1, 0, 0, 0, 0, 0)
2954 # sage: J2.one().to_vector()
2956 # sage: iota_right(J2.one()).to_vector()
2957 # (0, 0, 0, 1, 0, 1)
2958 # sage: J.one().to_vector()
2959 # (1, 0, 0, 1, 0, 1)
2963 # Composing a projection with the corresponding inclusion should
2964 # produce the identity map, and mismatching them should produce
2967 # sage: set_random_seed()
2968 # sage: J1 = random_eja()
2969 # sage: J2 = random_eja()
2970 # sage: J = DirectSumEJA(J1,J2)
2971 # sage: (iota_left, iota_right) = J.inclusions()
2972 # sage: (pi_left, pi_right) = J.projections()
2973 # sage: pi_left*iota_left == J1.one().operator()
2975 # sage: pi_right*iota_right == J2.one().operator()
2977 # sage: (pi_left*iota_right).is_zero()
2979 # sage: (pi_right*iota_left).is_zero()
2983 # (J1,J2) = self.factors()
2984 # m = J1.dimension()
2985 # n = J2.dimension()
2986 # V_basis = self.vector_space().basis()
2987 # # Need to specify the dimensions explicitly so that we don't
2988 # # wind up with a zero-by-zero matrix when we want e.g. a
2989 # # two-by-zero matrix (important for composing things).
2990 # I1 = matrix.column(self.base_ring(), m, m+n, V_basis[:m])
2991 # I2 = matrix.column(self.base_ring(), n, m+n, V_basis[m:])
2992 # iota_left = FiniteDimensionalEJAOperator(J1,self,I1)
2993 # iota_right = FiniteDimensionalEJAOperator(J2,self,I2)
2994 # return (iota_left, iota_right)
2996 # def inner_product(self, x, y):
2998 # The standard Cartesian inner-product.
3000 # We project ``x`` and ``y`` onto our factors, and add up the
3001 # inner-products from the subalgebras.
3006 # sage: from mjo.eja.eja_algebra import (HadamardEJA,
3007 # ....: QuaternionHermitianEJA,
3008 # ....: DirectSumEJA)
3012 # sage: J1 = HadamardEJA(3,field=QQ)
3013 # sage: J2 = QuaternionHermitianEJA(2,field=QQ,orthonormalize=False)
3014 # sage: J = DirectSumEJA(J1,J2)
3015 # sage: x1 = J1.one()
3017 # sage: y1 = J2.one()
3019 # sage: x1.inner_product(x2)
3021 # sage: y1.inner_product(y2)
3023 # sage: J.one().inner_product(J.one())
3027 # (pi_left, pi_right) = self.projections()
3033 # return (x1.inner_product(y1) + x2.inner_product(y2))
3037 random_eja
= ConcreteEJA
.random_instance