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
.combinat
.free_module
import CombinatorialFreeModule
24 from sage
.matrix
.constructor
import matrix
25 from sage
.matrix
.matrix_space
import MatrixSpace
26 from sage
.misc
.cachefunc
import cached_method
27 from sage
.misc
.table
import table
28 from sage
.modules
.free_module
import FreeModule
, VectorSpace
29 from sage
.rings
.all
import (ZZ
, QQ
, AA
, QQbar
, RR
, RLF
, CLF
,
32 from mjo
.eja
.eja_element
import FiniteDimensionalEJAElement
33 from mjo
.eja
.eja_operator
import FiniteDimensionalEJAOperator
34 from mjo
.eja
.eja_utils
import _mat2vec
36 class FiniteDimensionalEJA(CombinatorialFreeModule
):
38 A finite-dimensional Euclidean Jordan algebra.
42 - basis -- a tuple of basis elements in their matrix form.
44 - jordan_product -- function of two elements (in matrix form)
45 that returns their jordan product in this algebra; this will
46 be applied to ``basis`` to compute a multiplication table for
49 - inner_product -- function of two elements (in matrix form) that
50 returns their inner product. This will be applied to ``basis`` to
51 compute an inner-product table (basically a matrix) for this algebra.
54 Element
= FiniteDimensionalEJAElement
68 if not field
.is_subring(RR
):
69 # Note: this does return true for the real algebraic
70 # field, the rationals, and any quadratic field where
71 # we've specified a real embedding.
72 raise ValueError("scalar field is not real")
74 # If the basis given to us wasn't over the field that it's
75 # supposed to be over, fix that. Or, you know, crash.
76 basis
= tuple( b
.change_ring(field
) for b
in basis
)
79 # Check commutativity of the Jordan and inner-products.
80 # This has to be done before we build the multiplication
81 # and inner-product tables/matrices, because we take
82 # advantage of symmetry in the process.
83 if not all( jordan_product(bi
,bj
) == jordan_product(bj
,bi
)
86 raise ValueError("Jordan product is not commutative")
88 if not all( inner_product(bi
,bj
) == inner_product(bj
,bi
)
91 raise ValueError("inner-product is not commutative")
94 category
= MagmaticAlgebras(field
).FiniteDimensional()
95 category
= category
.WithBasis().Unital()
97 # Element subalgebras can take advantage of this.
98 category
= category
.Associative()
100 # Call the superclass constructor so that we can use its from_vector()
101 # method to build our multiplication table.
103 super().__init
__(field
,
109 # Now comes all of the hard work. We'll be constructing an
110 # ambient vector space V that our (vectorized) basis lives in,
111 # as well as a subspace W of V spanned by those (vectorized)
112 # basis elements. The W-coordinates are the coefficients that
113 # we see in things like x = 1*e1 + 2*e2.
118 # Works on both column and square matrices...
119 degree
= len(basis
[0].list())
121 # Build an ambient space that fits our matrix basis when
122 # written out as "long vectors."
123 V
= VectorSpace(field
, degree
)
125 # The matrix that will hole the orthonormal -> unorthonormal
126 # coordinate transformation.
127 self
._deortho
_matrix
= None
130 # Save a copy of the un-orthonormalized basis for later.
131 # Convert it to ambient V (vector) coordinates while we're
132 # at it, because we'd have to do it later anyway.
133 deortho_vector_basis
= tuple( V(b
.list()) for b
in basis
)
135 from mjo
.eja
.eja_utils
import gram_schmidt
136 basis
= tuple(gram_schmidt(basis
, inner_product
))
138 # Save the (possibly orthonormalized) matrix basis for
140 self
._matrix
_basis
= basis
142 # Now create the vector space for the algebra, which will have
143 # its own set of non-ambient coordinates (in terms of the
145 vector_basis
= tuple( V(b
.list()) for b
in basis
)
146 W
= V
.span_of_basis( vector_basis
, check
=check_axioms
)
149 # Now "W" is the vector space of our algebra coordinates. The
150 # variables "X1", "X2",... refer to the entries of vectors in
151 # W. Thus to convert back and forth between the orthonormal
152 # coordinates and the given ones, we need to stick the original
154 U
= V
.span_of_basis( deortho_vector_basis
, check
=check_axioms
)
155 self
._deortho
_matrix
= matrix( U
.coordinate_vector(q
)
156 for q
in vector_basis
)
159 # Now we actually compute the multiplication and inner-product
160 # tables/matrices using the possibly-orthonormalized basis.
161 self
._inner
_product
_matrix
= matrix
.identity(field
, n
)
162 self
._multiplication
_table
= [ [0 for j
in range(i
+1)]
165 # Note: the Jordan and inner-products are defined in terms
166 # of the ambient basis. It's important that their arguments
167 # are in ambient coordinates as well.
170 # ortho basis w.r.t. ambient coords
174 # The jordan product returns a matrixy answer, so we
175 # have to convert it to the algebra coordinates.
176 elt
= jordan_product(q_i
, q_j
)
177 elt
= W
.coordinate_vector(V(elt
.list()))
178 self
._multiplication
_table
[i
][j
] = self
.from_vector(elt
)
180 if not orthonormalize
:
181 # If we're orthonormalizing the basis with respect
182 # to an inner-product, then the inner-product
183 # matrix with respect to the resulting basis is
184 # just going to be the identity.
185 ip
= inner_product(q_i
, q_j
)
186 self
._inner
_product
_matrix
[i
,j
] = ip
187 self
._inner
_product
_matrix
[j
,i
] = ip
189 self
._inner
_product
_matrix
._cache
= {'hermitian': True}
190 self
._inner
_product
_matrix
.set_immutable()
193 if not self
._is
_jordanian
():
194 raise ValueError("Jordan identity does not hold")
195 if not self
._inner
_product
_is
_associative
():
196 raise ValueError("inner product is not associative")
199 def _coerce_map_from_base_ring(self
):
201 Disable the map from the base ring into the algebra.
203 Performing a nonsense conversion like this automatically
204 is counterpedagogical. The fallback is to try the usual
205 element constructor, which should also fail.
209 sage: from mjo.eja.eja_algebra import random_eja
213 sage: set_random_seed()
214 sage: J = random_eja()
216 Traceback (most recent call last):
218 ValueError: not an element of this algebra
224 def product_on_basis(self
, i
, j
):
225 # We only stored the lower-triangular portion of the
226 # multiplication table.
228 return self
._multiplication
_table
[i
][j
]
230 return self
._multiplication
_table
[j
][i
]
232 def inner_product(self
, x
, y
):
234 The inner product associated with this Euclidean Jordan algebra.
236 Defaults to the trace inner product, but can be overridden by
237 subclasses if they are sure that the necessary properties are
242 sage: from mjo.eja.eja_algebra import (random_eja,
244 ....: BilinearFormEJA)
248 Our inner product is "associative," which means the following for
249 a symmetric bilinear form::
251 sage: set_random_seed()
252 sage: J = random_eja()
253 sage: x,y,z = J.random_elements(3)
254 sage: (x*y).inner_product(z) == y.inner_product(x*z)
259 Ensure that this is the usual inner product for the algebras
262 sage: set_random_seed()
263 sage: J = HadamardEJA.random_instance()
264 sage: x,y = J.random_elements(2)
265 sage: actual = x.inner_product(y)
266 sage: expected = x.to_vector().inner_product(y.to_vector())
267 sage: actual == expected
270 Ensure that this is one-half of the trace inner-product in a
271 BilinearFormEJA that isn't just the reals (when ``n`` isn't
272 one). This is in Faraut and Koranyi, and also my "On the
275 sage: set_random_seed()
276 sage: J = BilinearFormEJA.random_instance()
277 sage: n = J.dimension()
278 sage: x = J.random_element()
279 sage: y = J.random_element()
280 sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
283 B
= self
._inner
_product
_matrix
284 return (B
*x
.to_vector()).inner_product(y
.to_vector())
287 def _is_commutative(self
):
289 Whether or not this algebra's multiplication table is commutative.
291 This method should of course always return ``True``, unless
292 this algebra was constructed with ``check_axioms=False`` and
293 passed an invalid multiplication table.
295 return all( self
.product_on_basis(i
,j
) == self
.product_on_basis(i
,j
)
296 for i
in range(self
.dimension())
297 for j
in range(self
.dimension()) )
299 def _is_jordanian(self
):
301 Whether or not this algebra's multiplication table respects the
302 Jordan identity `(x^{2})(xy) = x(x^{2}y)`.
304 We only check one arrangement of `x` and `y`, so for a
305 ``True`` result to be truly true, you should also check
306 :meth:`_is_commutative`. This method should of course always
307 return ``True``, unless this algebra was constructed with
308 ``check_axioms=False`` and passed an invalid multiplication table.
310 return all( (self
.monomial(i
)**2)*(self
.monomial(i
)*self
.monomial(j
))
312 (self
.monomial(i
))*((self
.monomial(i
)**2)*self
.monomial(j
))
313 for i
in range(self
.dimension())
314 for j
in range(self
.dimension()) )
316 def _inner_product_is_associative(self
):
318 Return whether or not this algebra's inner product `B` is
319 associative; that is, whether or not `B(xy,z) = B(x,yz)`.
321 This method should of course always return ``True``, unless
322 this algebra was constructed with ``check_axioms=False`` and
323 passed an invalid Jordan or inner-product.
326 # Used to check whether or not something is zero in an inexact
327 # ring. This number is sufficient to allow the construction of
328 # QuaternionHermitianEJA(2, field=RDF) with check_axioms=True.
331 for i
in range(self
.dimension()):
332 for j
in range(self
.dimension()):
333 for k
in range(self
.dimension()):
337 diff
= (x
*y
).inner_product(z
) - x
.inner_product(y
*z
)
339 if self
.base_ring().is_exact():
343 if diff
.abs() > epsilon
:
348 def _element_constructor_(self
, elt
):
350 Construct an element of this algebra from its vector or matrix
353 This gets called only after the parent element _call_ method
354 fails to find a coercion for the argument.
358 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
360 ....: RealSymmetricEJA)
364 The identity in `S^n` is converted to the identity in the EJA::
366 sage: J = RealSymmetricEJA(3)
367 sage: I = matrix.identity(QQ,3)
368 sage: J(I) == J.one()
371 This skew-symmetric matrix can't be represented in the EJA::
373 sage: J = RealSymmetricEJA(3)
374 sage: A = matrix(QQ,3, lambda i,j: i-j)
376 Traceback (most recent call last):
378 ValueError: not an element of this algebra
382 Ensure that we can convert any element of the two non-matrix
383 simple algebras (whose matrix representations are columns)
384 back and forth faithfully::
386 sage: set_random_seed()
387 sage: J = HadamardEJA.random_instance()
388 sage: x = J.random_element()
389 sage: J(x.to_vector().column()) == x
391 sage: J = JordanSpinEJA.random_instance()
392 sage: x = J.random_element()
393 sage: J(x.to_vector().column()) == x
397 msg
= "not an element of this algebra"
399 # The superclass implementation of random_element()
400 # needs to be able to coerce "0" into the algebra.
402 elif elt
in self
.base_ring():
403 # Ensure that no base ring -> algebra coercion is performed
404 # by this method. There's some stupidity in sage that would
405 # otherwise propagate to this method; for example, sage thinks
406 # that the integer 3 belongs to the space of 2-by-2 matrices.
407 raise ValueError(msg
)
411 except (AttributeError, TypeError):
412 # Try to convert a vector into a column-matrix
415 if elt
not in self
.matrix_space():
416 raise ValueError(msg
)
418 # Thanks for nothing! Matrix spaces aren't vector spaces in
419 # Sage, so we have to figure out its matrix-basis coordinates
420 # ourselves. We use the basis space's ring instead of the
421 # element's ring because the basis space might be an algebraic
422 # closure whereas the base ring of the 3-by-3 identity matrix
423 # could be QQ instead of QQbar.
425 # We pass check=False because the matrix basis is "guaranteed"
426 # to be linearly independent... right? Ha ha.
427 V
= VectorSpace(self
.base_ring(), elt
.nrows()*elt
.ncols())
428 W
= V
.span_of_basis( (_mat2vec(s
) for s
in self
.matrix_basis()),
432 coords
= W
.coordinate_vector(_mat2vec(elt
))
433 except ArithmeticError: # vector is not in free module
434 raise ValueError(msg
)
436 return self
.from_vector(coords
)
440 Return a string representation of ``self``.
444 sage: from mjo.eja.eja_algebra import JordanSpinEJA
448 Ensure that it says what we think it says::
450 sage: JordanSpinEJA(2, field=AA)
451 Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
452 sage: JordanSpinEJA(3, field=RDF)
453 Euclidean Jordan algebra of dimension 3 over Real Double Field
456 fmt
= "Euclidean Jordan algebra of dimension {} over {}"
457 return fmt
.format(self
.dimension(), self
.base_ring())
461 def characteristic_polynomial_of(self
):
463 Return the algebra's "characteristic polynomial of" function,
464 which is itself a multivariate polynomial that, when evaluated
465 at the coordinates of some algebra element, returns that
466 element's characteristic polynomial.
468 The resulting polynomial has `n+1` variables, where `n` is the
469 dimension of this algebra. The first `n` variables correspond to
470 the coordinates of an algebra element: when evaluated at the
471 coordinates of an algebra element with respect to a certain
472 basis, the result is a univariate polynomial (in the one
473 remaining variable ``t``), namely the characteristic polynomial
478 sage: from mjo.eja.eja_algebra import JordanSpinEJA, TrivialEJA
482 The characteristic polynomial in the spin algebra is given in
483 Alizadeh, Example 11.11::
485 sage: J = JordanSpinEJA(3)
486 sage: p = J.characteristic_polynomial_of(); p
487 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
488 sage: xvec = J.one().to_vector()
492 By definition, the characteristic polynomial is a monic
493 degree-zero polynomial in a rank-zero algebra. Note that
494 Cayley-Hamilton is indeed satisfied since the polynomial
495 ``1`` evaluates to the identity element of the algebra on
498 sage: J = TrivialEJA()
499 sage: J.characteristic_polynomial_of()
506 # The list of coefficient polynomials a_0, a_1, a_2, ..., a_(r-1).
507 a
= self
._charpoly
_coefficients
()
509 # We go to a bit of trouble here to reorder the
510 # indeterminates, so that it's easier to evaluate the
511 # characteristic polynomial at x's coordinates and get back
512 # something in terms of t, which is what we want.
513 S
= PolynomialRing(self
.base_ring(),'t')
517 S
= PolynomialRing(S
, R
.variable_names())
520 return (t
**r
+ sum( a
[k
]*(t
**k
) for k
in range(r
) ))
522 def coordinate_polynomial_ring(self
):
524 The multivariate polynomial ring in which this algebra's
525 :meth:`characteristic_polynomial_of` lives.
529 sage: from mjo.eja.eja_algebra import (HadamardEJA,
530 ....: RealSymmetricEJA)
534 sage: J = HadamardEJA(2)
535 sage: J.coordinate_polynomial_ring()
536 Multivariate Polynomial Ring in X1, X2...
537 sage: J = RealSymmetricEJA(3,field=QQ,orthonormalize=False)
538 sage: J.coordinate_polynomial_ring()
539 Multivariate Polynomial Ring in X1, X2, X3, X4, X5, X6...
542 var_names
= tuple( "X%d" % z
for z
in range(1, self
.dimension()+1) )
543 return PolynomialRing(self
.base_ring(), var_names
)
545 def inner_product(self
, x
, y
):
547 The inner product associated with this Euclidean Jordan algebra.
549 Defaults to the trace inner product, but can be overridden by
550 subclasses if they are sure that the necessary properties are
555 sage: from mjo.eja.eja_algebra import (random_eja,
557 ....: BilinearFormEJA)
561 Our inner product is "associative," which means the following for
562 a symmetric bilinear form::
564 sage: set_random_seed()
565 sage: J = random_eja()
566 sage: x,y,z = J.random_elements(3)
567 sage: (x*y).inner_product(z) == y.inner_product(x*z)
572 Ensure that this is the usual inner product for the algebras
575 sage: set_random_seed()
576 sage: J = HadamardEJA.random_instance()
577 sage: x,y = J.random_elements(2)
578 sage: actual = x.inner_product(y)
579 sage: expected = x.to_vector().inner_product(y.to_vector())
580 sage: actual == expected
583 Ensure that this is one-half of the trace inner-product in a
584 BilinearFormEJA that isn't just the reals (when ``n`` isn't
585 one). This is in Faraut and Koranyi, and also my "On the
588 sage: set_random_seed()
589 sage: J = BilinearFormEJA.random_instance()
590 sage: n = J.dimension()
591 sage: x = J.random_element()
592 sage: y = J.random_element()
593 sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
596 B
= self
._inner
_product
_matrix
597 return (B
*x
.to_vector()).inner_product(y
.to_vector())
600 def is_trivial(self
):
602 Return whether or not this algebra is trivial.
604 A trivial algebra contains only the zero element.
608 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
613 sage: J = ComplexHermitianEJA(3)
619 sage: J = TrivialEJA()
624 return self
.dimension() == 0
627 def multiplication_table(self
):
629 Return a visual representation of this algebra's multiplication
630 table (on basis elements).
634 sage: from mjo.eja.eja_algebra import JordanSpinEJA
638 sage: J = JordanSpinEJA(4)
639 sage: J.multiplication_table()
640 +----++----+----+----+----+
641 | * || e0 | e1 | e2 | e3 |
642 +====++====+====+====+====+
643 | e0 || e0 | e1 | e2 | e3 |
644 +----++----+----+----+----+
645 | e1 || e1 | e0 | 0 | 0 |
646 +----++----+----+----+----+
647 | e2 || e2 | 0 | e0 | 0 |
648 +----++----+----+----+----+
649 | e3 || e3 | 0 | 0 | e0 |
650 +----++----+----+----+----+
654 # Prepend the header row.
655 M
= [["*"] + list(self
.gens())]
657 # And to each subsequent row, prepend an entry that belongs to
658 # the left-side "header column."
659 M
+= [ [self
.monomial(i
)] + [ self
.product_on_basis(i
,j
)
663 return table(M
, header_row
=True, header_column
=True, frame
=True)
666 def matrix_basis(self
):
668 Return an (often more natural) representation of this algebras
669 basis as an ordered tuple of matrices.
671 Every finite-dimensional Euclidean Jordan Algebra is a, up to
672 Jordan isomorphism, a direct sum of five simple
673 algebras---four of which comprise Hermitian matrices. And the
674 last type of algebra can of course be thought of as `n`-by-`1`
675 column matrices (ambiguusly called column vectors) to avoid
676 special cases. As a result, matrices (and column vectors) are
677 a natural representation format for Euclidean Jordan algebra
680 But, when we construct an algebra from a basis of matrices,
681 those matrix representations are lost in favor of coordinate
682 vectors *with respect to* that basis. We could eventually
683 convert back if we tried hard enough, but having the original
684 representations handy is valuable enough that we simply store
685 them and return them from this method.
687 Why implement this for non-matrix algebras? Avoiding special
688 cases for the :class:`BilinearFormEJA` pays with simplicity in
689 its own right. But mainly, we would like to be able to assume
690 that elements of a :class:`DirectSumEJA` can be displayed
691 nicely, without having to have special classes for direct sums
692 one of whose components was a matrix algebra.
696 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
697 ....: RealSymmetricEJA)
701 sage: J = RealSymmetricEJA(2)
703 Finite family {0: e0, 1: e1, 2: e2}
704 sage: J.matrix_basis()
706 [1 0] [ 0 0.7071067811865475?] [0 0]
707 [0 0], [0.7071067811865475? 0], [0 1]
712 sage: J = JordanSpinEJA(2)
714 Finite family {0: e0, 1: e1}
715 sage: J.matrix_basis()
721 return self
._matrix
_basis
724 def matrix_space(self
):
726 Return the matrix space in which this algebra's elements live, if
727 we think of them as matrices (including column vectors of the
730 Generally this will be an `n`-by-`1` column-vector space,
731 except when the algebra is trivial. There it's `n`-by-`n`
732 (where `n` is zero), to ensure that two elements of the matrix
733 space (empty matrices) can be multiplied.
735 Matrix algebras override this with something more useful.
737 if self
.is_trivial():
738 return MatrixSpace(self
.base_ring(), 0)
740 return self
.matrix_basis()[0].parent()
746 Return the unit element of this algebra.
750 sage: from mjo.eja.eja_algebra import (HadamardEJA,
755 We can compute unit element in the Hadamard EJA::
757 sage: J = HadamardEJA(5)
759 e0 + e1 + e2 + e3 + e4
761 The unit element in the Hadamard EJA is inherited in the
762 subalgebras generated by its elements::
764 sage: J = HadamardEJA(5)
766 e0 + e1 + e2 + e3 + e4
767 sage: x = sum(J.gens())
768 sage: A = x.subalgebra_generated_by(orthonormalize=False)
771 sage: A.one().superalgebra_element()
772 e0 + e1 + e2 + e3 + e4
776 The identity element acts like the identity, regardless of
777 whether or not we orthonormalize::
779 sage: set_random_seed()
780 sage: J = random_eja()
781 sage: x = J.random_element()
782 sage: J.one()*x == x and x*J.one() == x
784 sage: A = x.subalgebra_generated_by()
785 sage: y = A.random_element()
786 sage: A.one()*y == y and y*A.one() == y
791 sage: set_random_seed()
792 sage: J = random_eja(field=QQ, orthonormalize=False)
793 sage: x = J.random_element()
794 sage: J.one()*x == x and x*J.one() == x
796 sage: A = x.subalgebra_generated_by(orthonormalize=False)
797 sage: y = A.random_element()
798 sage: A.one()*y == y and y*A.one() == y
801 The matrix of the unit element's operator is the identity,
802 regardless of the base field and whether or not we
805 sage: set_random_seed()
806 sage: J = random_eja()
807 sage: actual = J.one().operator().matrix()
808 sage: expected = matrix.identity(J.base_ring(), J.dimension())
809 sage: actual == expected
811 sage: x = J.random_element()
812 sage: A = x.subalgebra_generated_by()
813 sage: actual = A.one().operator().matrix()
814 sage: expected = matrix.identity(A.base_ring(), A.dimension())
815 sage: actual == expected
820 sage: set_random_seed()
821 sage: J = random_eja(field=QQ, orthonormalize=False)
822 sage: actual = J.one().operator().matrix()
823 sage: expected = matrix.identity(J.base_ring(), J.dimension())
824 sage: actual == expected
826 sage: x = J.random_element()
827 sage: A = x.subalgebra_generated_by(orthonormalize=False)
828 sage: actual = A.one().operator().matrix()
829 sage: expected = matrix.identity(A.base_ring(), A.dimension())
830 sage: actual == expected
833 Ensure that the cached unit element (often precomputed by
834 hand) agrees with the computed one::
836 sage: set_random_seed()
837 sage: J = random_eja()
838 sage: cached = J.one()
839 sage: J.one.clear_cache()
840 sage: J.one() == cached
845 sage: set_random_seed()
846 sage: J = random_eja(field=QQ, orthonormalize=False)
847 sage: cached = J.one()
848 sage: J.one.clear_cache()
849 sage: J.one() == cached
853 # We can brute-force compute the matrices of the operators
854 # that correspond to the basis elements of this algebra.
855 # If some linear combination of those basis elements is the
856 # algebra identity, then the same linear combination of
857 # their matrices has to be the identity matrix.
859 # Of course, matrices aren't vectors in sage, so we have to
860 # appeal to the "long vectors" isometry.
861 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
863 # Now we use basic linear algebra to find the coefficients,
864 # of the matrices-as-vectors-linear-combination, which should
865 # work for the original algebra basis too.
866 A
= matrix(self
.base_ring(), oper_vecs
)
868 # We used the isometry on the left-hand side already, but we
869 # still need to do it for the right-hand side. Recall that we
870 # wanted something that summed to the identity matrix.
871 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
873 # Now if there's an identity element in the algebra, this
874 # should work. We solve on the left to avoid having to
875 # transpose the matrix "A".
876 return self
.from_vector(A
.solve_left(b
))
879 def peirce_decomposition(self
, c
):
881 The Peirce decomposition of this algebra relative to the
884 In the future, this can be extended to a complete system of
885 orthogonal idempotents.
889 - ``c`` -- an idempotent of this algebra.
893 A triple (J0, J5, J1) containing two subalgebras and one subspace
896 - ``J0`` -- the algebra on the eigenspace of ``c.operator()``
897 corresponding to the eigenvalue zero.
899 - ``J5`` -- the eigenspace (NOT a subalgebra) of ``c.operator()``
900 corresponding to the eigenvalue one-half.
902 - ``J1`` -- the algebra on the eigenspace of ``c.operator()``
903 corresponding to the eigenvalue one.
905 These are the only possible eigenspaces for that operator, and this
906 algebra is a direct sum of them. The spaces ``J0`` and ``J1`` are
907 orthogonal, and are subalgebras of this algebra with the appropriate
912 sage: from mjo.eja.eja_algebra import random_eja, RealSymmetricEJA
916 The canonical example comes from the symmetric matrices, which
917 decompose into diagonal and off-diagonal parts::
919 sage: J = RealSymmetricEJA(3)
920 sage: C = matrix(QQ, [ [1,0,0],
924 sage: J0,J5,J1 = J.peirce_decomposition(c)
926 Euclidean Jordan algebra of dimension 1...
928 Vector space of degree 6 and dimension 2...
930 Euclidean Jordan algebra of dimension 3...
931 sage: J0.one().to_matrix()
935 sage: orig_df = AA.options.display_format
936 sage: AA.options.display_format = 'radical'
937 sage: J.from_vector(J5.basis()[0]).to_matrix()
941 sage: J.from_vector(J5.basis()[1]).to_matrix()
945 sage: AA.options.display_format = orig_df
946 sage: J1.one().to_matrix()
953 Every algebra decomposes trivially with respect to its identity
956 sage: set_random_seed()
957 sage: J = random_eja()
958 sage: J0,J5,J1 = J.peirce_decomposition(J.one())
959 sage: J0.dimension() == 0 and J5.dimension() == 0
961 sage: J1.superalgebra() == J and J1.dimension() == J.dimension()
964 The decomposition is into eigenspaces, and its components are
965 therefore necessarily orthogonal. Moreover, the identity
966 elements in the two subalgebras are the projections onto their
967 respective subspaces of the superalgebra's identity element::
969 sage: set_random_seed()
970 sage: J = random_eja()
971 sage: x = J.random_element()
972 sage: if not J.is_trivial():
973 ....: while x.is_nilpotent():
974 ....: x = J.random_element()
975 sage: c = x.subalgebra_idempotent()
976 sage: J0,J5,J1 = J.peirce_decomposition(c)
978 sage: for (w,y,z) in zip(J0.basis(), J5.basis(), J1.basis()):
979 ....: w = w.superalgebra_element()
980 ....: y = J.from_vector(y)
981 ....: z = z.superalgebra_element()
982 ....: ipsum += w.inner_product(y).abs()
983 ....: ipsum += w.inner_product(z).abs()
984 ....: ipsum += y.inner_product(z).abs()
987 sage: J1(c) == J1.one()
989 sage: J0(J.one() - c) == J0.one()
993 if not c
.is_idempotent():
994 raise ValueError("element is not idempotent: %s" % c
)
996 from mjo
.eja
.eja_subalgebra
import FiniteDimensionalEJASubalgebra
998 # Default these to what they should be if they turn out to be
999 # trivial, because eigenspaces_left() won't return eigenvalues
1000 # corresponding to trivial spaces (e.g. it returns only the
1001 # eigenspace corresponding to lambda=1 if you take the
1002 # decomposition relative to the identity element).
1003 trivial
= FiniteDimensionalEJASubalgebra(self
, ())
1004 J0
= trivial
# eigenvalue zero
1005 J5
= VectorSpace(self
.base_ring(), 0) # eigenvalue one-half
1006 J1
= trivial
# eigenvalue one
1008 for (eigval
, eigspace
) in c
.operator().matrix().right_eigenspaces():
1009 if eigval
== ~
(self
.base_ring()(2)):
1012 gens
= tuple( self
.from_vector(b
) for b
in eigspace
.basis() )
1013 subalg
= FiniteDimensionalEJASubalgebra(self
,
1021 raise ValueError("unexpected eigenvalue: %s" % eigval
)
1026 def random_element(self
, thorough
=False):
1028 Return a random element of this algebra.
1030 Our algebra superclass method only returns a linear
1031 combination of at most two basis elements. We instead
1032 want the vector space "random element" method that
1033 returns a more diverse selection.
1037 - ``thorough`` -- (boolean; default False) whether or not we
1038 should generate irrational coefficients for the random
1039 element when our base ring is irrational; this slows the
1040 algebra operations to a crawl, but any truly random method
1044 # For a general base ring... maybe we can trust this to do the
1045 # right thing? Unlikely, but.
1046 V
= self
.vector_space()
1047 v
= V
.random_element()
1049 if self
.base_ring() is AA
:
1050 # The "random element" method of the algebraic reals is
1051 # stupid at the moment, and only returns integers between
1052 # -2 and 2, inclusive:
1054 # https://trac.sagemath.org/ticket/30875
1056 # Instead, we implement our own "random vector" method,
1057 # and then coerce that into the algebra. We use the vector
1058 # space degree here instead of the dimension because a
1059 # subalgebra could (for example) be spanned by only two
1060 # vectors, each with five coordinates. We need to
1061 # generate all five coordinates.
1063 v
*= QQbar
.random_element().real()
1065 v
*= QQ
.random_element()
1067 return self
.from_vector(V
.coordinate_vector(v
))
1069 def random_elements(self
, count
, thorough
=False):
1071 Return ``count`` random elements as a tuple.
1075 - ``thorough`` -- (boolean; default False) whether or not we
1076 should generate irrational coefficients for the random
1077 elements when our base ring is irrational; this slows the
1078 algebra operations to a crawl, but any truly random method
1083 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1087 sage: J = JordanSpinEJA(3)
1088 sage: x,y,z = J.random_elements(3)
1089 sage: all( [ x in J, y in J, z in J ])
1091 sage: len( J.random_elements(10) ) == 10
1095 return tuple( self
.random_element(thorough
)
1096 for idx
in range(count
) )
1100 def _charpoly_coefficients(self
):
1102 The `r` polynomial coefficients of the "characteristic polynomial
1107 sage: from mjo.eja.eja_algebra import random_eja
1111 The theory shows that these are all homogeneous polynomials of
1114 sage: set_random_seed()
1115 sage: J = random_eja()
1116 sage: all(p.is_homogeneous() for p in J._charpoly_coefficients())
1120 n
= self
.dimension()
1121 R
= self
.coordinate_polynomial_ring()
1123 F
= R
.fraction_field()
1126 # From a result in my book, these are the entries of the
1127 # basis representation of L_x.
1128 return sum( vars[k
]*self
.monomial(k
).operator().matrix()[i
,j
]
1131 L_x
= matrix(F
, n
, n
, L_x_i_j
)
1134 if self
.rank
.is_in_cache():
1136 # There's no need to pad the system with redundant
1137 # columns if we *know* they'll be redundant.
1140 # Compute an extra power in case the rank is equal to
1141 # the dimension (otherwise, we would stop at x^(r-1)).
1142 x_powers
= [ (L_x
**k
)*self
.one().to_vector()
1143 for k
in range(n
+1) ]
1144 A
= matrix
.column(F
, x_powers
[:n
])
1145 AE
= A
.extended_echelon_form()
1152 # The theory says that only the first "r" coefficients are
1153 # nonzero, and they actually live in the original polynomial
1154 # ring and not the fraction field. We negate them because in
1155 # the actual characteristic polynomial, they get moved to the
1156 # other side where x^r lives. We don't bother to trim A_rref
1157 # down to a square matrix and solve the resulting system,
1158 # because the upper-left r-by-r portion of A_rref is
1159 # guaranteed to be the identity matrix, so e.g.
1161 # A_rref.solve_right(Y)
1163 # would just be returning Y.
1164 return (-E
*b
)[:r
].change_ring(R
)
1169 Return the rank of this EJA.
1171 This is a cached method because we know the rank a priori for
1172 all of the algebras we can construct. Thus we can avoid the
1173 expensive ``_charpoly_coefficients()`` call unless we truly
1174 need to compute the whole characteristic polynomial.
1178 sage: from mjo.eja.eja_algebra import (HadamardEJA,
1179 ....: JordanSpinEJA,
1180 ....: RealSymmetricEJA,
1181 ....: ComplexHermitianEJA,
1182 ....: QuaternionHermitianEJA,
1187 The rank of the Jordan spin algebra is always two::
1189 sage: JordanSpinEJA(2).rank()
1191 sage: JordanSpinEJA(3).rank()
1193 sage: JordanSpinEJA(4).rank()
1196 The rank of the `n`-by-`n` Hermitian real, complex, or
1197 quaternion matrices is `n`::
1199 sage: RealSymmetricEJA(4).rank()
1201 sage: ComplexHermitianEJA(3).rank()
1203 sage: QuaternionHermitianEJA(2).rank()
1208 Ensure that every EJA that we know how to construct has a
1209 positive integer rank, unless the algebra is trivial in
1210 which case its rank will be zero::
1212 sage: set_random_seed()
1213 sage: J = random_eja()
1217 sage: r > 0 or (r == 0 and J.is_trivial())
1220 Ensure that computing the rank actually works, since the ranks
1221 of all simple algebras are known and will be cached by default::
1223 sage: set_random_seed() # long time
1224 sage: J = random_eja() # long time
1225 sage: cached = J.rank() # long time
1226 sage: J.rank.clear_cache() # long time
1227 sage: J.rank() == cached # long time
1231 return len(self
._charpoly
_coefficients
())
1234 def vector_space(self
):
1236 Return the vector space that underlies this algebra.
1240 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1244 sage: J = RealSymmetricEJA(2)
1245 sage: J.vector_space()
1246 Vector space of dimension 3 over...
1249 return self
.zero().to_vector().parent().ambient_vector_space()
1252 Element
= FiniteDimensionalEJAElement
1254 class RationalBasisEJA(FiniteDimensionalEJA
):
1256 New class for algebras whose supplied basis elements have all rational entries.
1260 sage: from mjo.eja.eja_algebra import BilinearFormEJA
1264 The supplied basis is orthonormalized by default::
1266 sage: B = matrix(QQ, [[1, 0, 0], [0, 25, -32], [0, -32, 41]])
1267 sage: J = BilinearFormEJA(B)
1268 sage: J.matrix_basis()
1285 # Abuse the check_field parameter to check that the entries of
1286 # out basis (in ambient coordinates) are in the field QQ.
1287 if not all( all(b_i
in QQ
for b_i
in b
.list()) for b
in basis
):
1288 raise TypeError("basis not rational")
1290 self
._rational
_algebra
= None
1292 # There's no point in constructing the extra algebra if this
1293 # one is already rational.
1295 # Note: the same Jordan and inner-products work here,
1296 # because they are necessarily defined with respect to
1297 # ambient coordinates and not any particular basis.
1298 self
._rational
_algebra
= FiniteDimensionalEJA(
1303 orthonormalize
=False,
1307 super().__init
__(basis
,
1311 check_field
=check_field
,
1315 def _charpoly_coefficients(self
):
1319 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
1320 ....: JordanSpinEJA)
1324 The base ring of the resulting polynomial coefficients is what
1325 it should be, and not the rationals (unless the algebra was
1326 already over the rationals)::
1328 sage: J = JordanSpinEJA(3)
1329 sage: J._charpoly_coefficients()
1330 (X1^2 - X2^2 - X3^2, -2*X1)
1331 sage: a0 = J._charpoly_coefficients()[0]
1333 Algebraic Real Field
1334 sage: a0.base_ring()
1335 Algebraic Real Field
1338 if self
._rational
_algebra
is None:
1339 # There's no need to construct *another* algebra over the
1340 # rationals if this one is already over the
1341 # rationals. Likewise, if we never orthonormalized our
1342 # basis, we might as well just use the given one.
1343 return super()._charpoly
_coefficients
()
1345 # Do the computation over the rationals. The answer will be
1346 # the same, because all we've done is a change of basis.
1347 # Then, change back from QQ to our real base ring
1348 a
= ( a_i
.change_ring(self
.base_ring())
1349 for a_i
in self
._rational
_algebra
._charpoly
_coefficients
() )
1351 if self
._deortho
_matrix
is None:
1352 # This can happen if our base ring was, say, AA and we
1353 # chose not to (or didn't need to) orthonormalize. It's
1354 # still faster to do the computations over QQ even if
1355 # the numbers in the boxes stay the same.
1358 # Otherwise, convert the coordinate variables back to the
1359 # deorthonormalized ones.
1360 R
= self
.coordinate_polynomial_ring()
1361 from sage
.modules
.free_module_element
import vector
1362 X
= vector(R
, R
.gens())
1363 BX
= self
._deortho
_matrix
*X
1365 subs_dict
= { X[i]: BX[i] for i in range(len(X)) }
1366 return tuple( a_i
.subs(subs_dict
) for a_i
in a
)
1368 class ConcreteEJA(RationalBasisEJA
):
1370 A class for the Euclidean Jordan algebras that we know by name.
1372 These are the Jordan algebras whose basis, multiplication table,
1373 rank, and so on are known a priori. More to the point, they are
1374 the Euclidean Jordan algebras for which we are able to conjure up
1375 a "random instance."
1379 sage: from mjo.eja.eja_algebra import ConcreteEJA
1383 Our basis is normalized with respect to the algebra's inner
1384 product, unless we specify otherwise::
1386 sage: set_random_seed()
1387 sage: J = ConcreteEJA.random_instance()
1388 sage: all( b.norm() == 1 for b in J.gens() )
1391 Since our basis is orthonormal with respect to the algebra's inner
1392 product, and since we know that this algebra is an EJA, any
1393 left-multiplication operator's matrix will be symmetric because
1394 natural->EJA basis representation is an isometry and within the
1395 EJA the operator is self-adjoint by the Jordan axiom::
1397 sage: set_random_seed()
1398 sage: J = ConcreteEJA.random_instance()
1399 sage: x = J.random_element()
1400 sage: x.operator().is_self_adjoint()
1405 def _max_random_instance_size():
1407 Return an integer "size" that is an upper bound on the size of
1408 this algebra when it is used in a random test
1409 case. Unfortunately, the term "size" is ambiguous -- when
1410 dealing with `R^n` under either the Hadamard or Jordan spin
1411 product, the "size" refers to the dimension `n`. When dealing
1412 with a matrix algebra (real symmetric or complex/quaternion
1413 Hermitian), it refers to the size of the matrix, which is far
1414 less than the dimension of the underlying vector space.
1416 This method must be implemented in each subclass.
1418 raise NotImplementedError
1421 def random_instance(cls
, *args
, **kwargs
):
1423 Return a random instance of this type of algebra.
1425 This method should be implemented in each subclass.
1427 from sage
.misc
.prandom
import choice
1428 eja_class
= choice(cls
.__subclasses
__())
1430 # These all bubble up to the RationalBasisEJA superclass
1431 # constructor, so any (kw)args valid there are also valid
1433 return eja_class
.random_instance(*args
, **kwargs
)
1438 def dimension_over_reals():
1440 The dimension of this matrix's base ring over the reals.
1442 The reals are dimension one over themselves, obviously; that's
1443 just `\mathbb{R}^{1}`. Likewise, the complex numbers `a + bi`
1444 have dimension two. Finally, the quaternions have dimension
1445 four over the reals.
1447 This is used to determine the size of the matrix returned from
1448 :meth:`real_embed`, among other things.
1450 raise NotImplementedError
1453 def real_embed(cls
,M
):
1455 Embed the matrix ``M`` into a space of real matrices.
1457 The matrix ``M`` can have entries in any field at the moment:
1458 the real numbers, complex numbers, or quaternions. And although
1459 they are not a field, we can probably support octonions at some
1460 point, too. This function returns a real matrix that "acts like"
1461 the original with respect to matrix multiplication; i.e.
1463 real_embed(M*N) = real_embed(M)*real_embed(N)
1466 if M
.ncols() != M
.nrows():
1467 raise ValueError("the matrix 'M' must be square")
1472 def real_unembed(cls
,M
):
1474 The inverse of :meth:`real_embed`.
1476 if M
.ncols() != M
.nrows():
1477 raise ValueError("the matrix 'M' must be square")
1478 if not ZZ(M
.nrows()).mod(cls
.dimension_over_reals()).is_zero():
1479 raise ValueError("the matrix 'M' must be a real embedding")
1483 def jordan_product(X
,Y
):
1484 return (X
*Y
+ Y
*X
)/2
1487 def trace_inner_product(cls
,X
,Y
):
1489 Compute the trace inner-product of two real-embeddings.
1493 sage: from mjo.eja.eja_algebra import (RealSymmetricEJA,
1494 ....: ComplexHermitianEJA,
1495 ....: QuaternionHermitianEJA)
1499 This gives the same answer as it would if we computed the trace
1500 from the unembedded (original) matrices::
1502 sage: set_random_seed()
1503 sage: J = RealSymmetricEJA.random_instance()
1504 sage: x,y = J.random_elements(2)
1505 sage: Xe = x.to_matrix()
1506 sage: Ye = y.to_matrix()
1507 sage: X = J.real_unembed(Xe)
1508 sage: Y = J.real_unembed(Ye)
1509 sage: expected = (X*Y).trace()
1510 sage: actual = J.trace_inner_product(Xe,Ye)
1511 sage: actual == expected
1516 sage: set_random_seed()
1517 sage: J = ComplexHermitianEJA.random_instance()
1518 sage: x,y = J.random_elements(2)
1519 sage: Xe = x.to_matrix()
1520 sage: Ye = y.to_matrix()
1521 sage: X = J.real_unembed(Xe)
1522 sage: Y = J.real_unembed(Ye)
1523 sage: expected = (X*Y).trace().real()
1524 sage: actual = J.trace_inner_product(Xe,Ye)
1525 sage: actual == expected
1530 sage: set_random_seed()
1531 sage: J = QuaternionHermitianEJA.random_instance()
1532 sage: x,y = J.random_elements(2)
1533 sage: Xe = x.to_matrix()
1534 sage: Ye = y.to_matrix()
1535 sage: X = J.real_unembed(Xe)
1536 sage: Y = J.real_unembed(Ye)
1537 sage: expected = (X*Y).trace().coefficient_tuple()[0]
1538 sage: actual = J.trace_inner_product(Xe,Ye)
1539 sage: actual == expected
1543 Xu
= cls
.real_unembed(X
)
1544 Yu
= cls
.real_unembed(Y
)
1545 tr
= (Xu
*Yu
).trace()
1548 # Works in QQ, AA, RDF, et cetera.
1550 except AttributeError:
1551 # A quaternion doesn't have a real() method, but does
1552 # have coefficient_tuple() method that returns the
1553 # coefficients of 1, i, j, and k -- in that order.
1554 return tr
.coefficient_tuple()[0]
1557 class RealMatrixEJA(MatrixEJA
):
1559 def dimension_over_reals():
1563 class RealSymmetricEJA(ConcreteEJA
, RealMatrixEJA
):
1565 The rank-n simple EJA consisting of real symmetric n-by-n
1566 matrices, the usual symmetric Jordan product, and the trace inner
1567 product. It has dimension `(n^2 + n)/2` over the reals.
1571 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1575 sage: J = RealSymmetricEJA(2)
1576 sage: e0, e1, e2 = J.gens()
1584 In theory, our "field" can be any subfield of the reals::
1586 sage: RealSymmetricEJA(2, field=RDF)
1587 Euclidean Jordan algebra of dimension 3 over Real Double Field
1588 sage: RealSymmetricEJA(2, field=RR)
1589 Euclidean Jordan algebra of dimension 3 over Real Field with
1590 53 bits of precision
1594 The dimension of this algebra is `(n^2 + n) / 2`::
1596 sage: set_random_seed()
1597 sage: n_max = RealSymmetricEJA._max_random_instance_size()
1598 sage: n = ZZ.random_element(1, n_max)
1599 sage: J = RealSymmetricEJA(n)
1600 sage: J.dimension() == (n^2 + n)/2
1603 The Jordan multiplication is what we think it is::
1605 sage: set_random_seed()
1606 sage: J = RealSymmetricEJA.random_instance()
1607 sage: x,y = J.random_elements(2)
1608 sage: actual = (x*y).to_matrix()
1609 sage: X = x.to_matrix()
1610 sage: Y = y.to_matrix()
1611 sage: expected = (X*Y + Y*X)/2
1612 sage: actual == expected
1614 sage: J(expected) == x*y
1617 We can change the generator prefix::
1619 sage: RealSymmetricEJA(3, prefix='q').gens()
1620 (q0, q1, q2, q3, q4, q5)
1622 We can construct the (trivial) algebra of rank zero::
1624 sage: RealSymmetricEJA(0)
1625 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1629 def _denormalized_basis(cls
, n
):
1631 Return a basis for the space of real symmetric n-by-n matrices.
1635 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1639 sage: set_random_seed()
1640 sage: n = ZZ.random_element(1,5)
1641 sage: B = RealSymmetricEJA._denormalized_basis(n)
1642 sage: all( M.is_symmetric() for M in B)
1646 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1650 for j
in range(i
+1):
1651 Eij
= matrix(ZZ
, n
, lambda k
,l
: k
==i
and l
==j
)
1655 Sij
= Eij
+ Eij
.transpose()
1661 def _max_random_instance_size():
1662 return 4 # Dimension 10
1665 def random_instance(cls
, **kwargs
):
1667 Return a random instance of this type of algebra.
1669 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1670 return cls(n
, **kwargs
)
1672 def __init__(self
, n
, **kwargs
):
1673 # We know this is a valid EJA, but will double-check
1674 # if the user passes check_axioms=True.
1675 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
1677 super(RealSymmetricEJA
, self
).__init
__(self
._denormalized
_basis
(n
),
1678 self
.jordan_product
,
1679 self
.trace_inner_product
,
1682 # TODO: this could be factored out somehow, but is left here
1683 # because the MatrixEJA is not presently a subclass of the
1684 # FDEJA class that defines rank() and one().
1685 self
.rank
.set_cache(n
)
1686 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
1687 self
.one
.set_cache(self(idV
))
1691 class ComplexMatrixEJA(MatrixEJA
):
1692 # A manual dictionary-cache for the complex_extension() method,
1693 # since apparently @classmethods can't also be @cached_methods.
1694 _complex_extension
= {}
1697 def complex_extension(cls
,field
):
1699 The complex field that we embed/unembed, as an extension
1700 of the given ``field``.
1702 if field
in cls
._complex
_extension
:
1703 return cls
._complex
_extension
[field
]
1705 # Sage doesn't know how to adjoin the complex "i" (the root of
1706 # x^2 + 1) to a field in a general way. Here, we just enumerate
1707 # all of the cases that I have cared to support so far.
1709 # Sage doesn't know how to embed AA into QQbar, i.e. how
1710 # to adjoin sqrt(-1) to AA.
1712 elif not field
.is_exact():
1714 F
= field
.complex_field()
1716 # Works for QQ and... maybe some other fields.
1717 R
= PolynomialRing(field
, 'z')
1719 F
= field
.extension(z
**2 + 1, 'I', embedding
=CLF(-1).sqrt())
1721 cls
._complex
_extension
[field
] = F
1725 def dimension_over_reals():
1729 def real_embed(cls
,M
):
1731 Embed the n-by-n complex matrix ``M`` into the space of real
1732 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1733 bi` to the block matrix ``[[a,b],[-b,a]]``.
1737 sage: from mjo.eja.eja_algebra import ComplexMatrixEJA
1741 sage: F = QuadraticField(-1, 'I')
1742 sage: x1 = F(4 - 2*i)
1743 sage: x2 = F(1 + 2*i)
1746 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1747 sage: ComplexMatrixEJA.real_embed(M)
1756 Embedding is a homomorphism (isomorphism, in fact)::
1758 sage: set_random_seed()
1759 sage: n = ZZ.random_element(3)
1760 sage: F = QuadraticField(-1, 'I')
1761 sage: X = random_matrix(F, n)
1762 sage: Y = random_matrix(F, n)
1763 sage: Xe = ComplexMatrixEJA.real_embed(X)
1764 sage: Ye = ComplexMatrixEJA.real_embed(Y)
1765 sage: XYe = ComplexMatrixEJA.real_embed(X*Y)
1770 super(ComplexMatrixEJA
,cls
).real_embed(M
)
1773 # We don't need any adjoined elements...
1774 field
= M
.base_ring().base_ring()
1780 blocks
.append(matrix(field
, 2, [ [ a
, b
],
1783 return matrix
.block(field
, n
, blocks
)
1787 def real_unembed(cls
,M
):
1789 The inverse of _embed_complex_matrix().
1793 sage: from mjo.eja.eja_algebra import ComplexMatrixEJA
1797 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1798 ....: [-2, 1, -4, 3],
1799 ....: [ 9, 10, 11, 12],
1800 ....: [-10, 9, -12, 11] ])
1801 sage: ComplexMatrixEJA.real_unembed(A)
1803 [ 10*I + 9 12*I + 11]
1807 Unembedding is the inverse of embedding::
1809 sage: set_random_seed()
1810 sage: F = QuadraticField(-1, 'I')
1811 sage: M = random_matrix(F, 3)
1812 sage: Me = ComplexMatrixEJA.real_embed(M)
1813 sage: ComplexMatrixEJA.real_unembed(Me) == M
1817 super(ComplexMatrixEJA
,cls
).real_unembed(M
)
1819 d
= cls
.dimension_over_reals()
1820 F
= cls
.complex_extension(M
.base_ring())
1823 # Go top-left to bottom-right (reading order), converting every
1824 # 2-by-2 block we see to a single complex element.
1826 for k
in range(n
/d
):
1827 for j
in range(n
/d
):
1828 submat
= M
[d
*k
:d
*k
+d
,d
*j
:d
*j
+d
]
1829 if submat
[0,0] != submat
[1,1]:
1830 raise ValueError('bad on-diagonal submatrix')
1831 if submat
[0,1] != -submat
[1,0]:
1832 raise ValueError('bad off-diagonal submatrix')
1833 z
= submat
[0,0] + submat
[0,1]*i
1836 return matrix(F
, n
/d
, elements
)
1839 class ComplexHermitianEJA(ConcreteEJA
, ComplexMatrixEJA
):
1841 The rank-n simple EJA consisting of complex Hermitian n-by-n
1842 matrices over the real numbers, the usual symmetric Jordan product,
1843 and the real-part-of-trace inner product. It has dimension `n^2` over
1848 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1852 In theory, our "field" can be any subfield of the reals::
1854 sage: ComplexHermitianEJA(2, field=RDF)
1855 Euclidean Jordan algebra of dimension 4 over Real Double Field
1856 sage: ComplexHermitianEJA(2, field=RR)
1857 Euclidean Jordan algebra of dimension 4 over Real Field with
1858 53 bits of precision
1862 The dimension of this algebra is `n^2`::
1864 sage: set_random_seed()
1865 sage: n_max = ComplexHermitianEJA._max_random_instance_size()
1866 sage: n = ZZ.random_element(1, n_max)
1867 sage: J = ComplexHermitianEJA(n)
1868 sage: J.dimension() == n^2
1871 The Jordan multiplication is what we think it is::
1873 sage: set_random_seed()
1874 sage: J = ComplexHermitianEJA.random_instance()
1875 sage: x,y = J.random_elements(2)
1876 sage: actual = (x*y).to_matrix()
1877 sage: X = x.to_matrix()
1878 sage: Y = y.to_matrix()
1879 sage: expected = (X*Y + Y*X)/2
1880 sage: actual == expected
1882 sage: J(expected) == x*y
1885 We can change the generator prefix::
1887 sage: ComplexHermitianEJA(2, prefix='z').gens()
1890 We can construct the (trivial) algebra of rank zero::
1892 sage: ComplexHermitianEJA(0)
1893 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1898 def _denormalized_basis(cls
, n
):
1900 Returns a basis for the space of complex Hermitian n-by-n matrices.
1902 Why do we embed these? Basically, because all of numerical linear
1903 algebra assumes that you're working with vectors consisting of `n`
1904 entries from a field and scalars from the same field. There's no way
1905 to tell SageMath that (for example) the vectors contain complex
1906 numbers, while the scalar field is real.
1910 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1914 sage: set_random_seed()
1915 sage: n = ZZ.random_element(1,5)
1916 sage: B = ComplexHermitianEJA._denormalized_basis(n)
1917 sage: all( M.is_symmetric() for M in B)
1922 R
= PolynomialRing(field
, 'z')
1924 F
= field
.extension(z
**2 + 1, 'I')
1927 # This is like the symmetric case, but we need to be careful:
1929 # * We want conjugate-symmetry, not just symmetry.
1930 # * The diagonal will (as a result) be real.
1933 Eij
= matrix
.zero(F
,n
)
1935 for j
in range(i
+1):
1939 Sij
= cls
.real_embed(Eij
)
1942 # The second one has a minus because it's conjugated.
1943 Eij
[j
,i
] = 1 # Eij = Eij + Eij.transpose()
1944 Sij_real
= cls
.real_embed(Eij
)
1946 # Eij = I*Eij - I*Eij.transpose()
1949 Sij_imag
= cls
.real_embed(Eij
)
1955 # Since we embedded these, we can drop back to the "field" that we
1956 # started with instead of the complex extension "F".
1957 return tuple( s
.change_ring(field
) for s
in S
)
1960 def __init__(self
, n
, **kwargs
):
1961 # We know this is a valid EJA, but will double-check
1962 # if the user passes check_axioms=True.
1963 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
1965 super(ComplexHermitianEJA
, self
).__init
__(self
._denormalized
_basis
(n
),
1966 self
.jordan_product
,
1967 self
.trace_inner_product
,
1969 # TODO: this could be factored out somehow, but is left here
1970 # because the MatrixEJA is not presently a subclass of the
1971 # FDEJA class that defines rank() and one().
1972 self
.rank
.set_cache(n
)
1973 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
1974 self
.one
.set_cache(self(idV
))
1977 def _max_random_instance_size():
1978 return 3 # Dimension 9
1981 def random_instance(cls
, **kwargs
):
1983 Return a random instance of this type of algebra.
1985 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1986 return cls(n
, **kwargs
)
1988 class QuaternionMatrixEJA(MatrixEJA
):
1990 # A manual dictionary-cache for the quaternion_extension() method,
1991 # since apparently @classmethods can't also be @cached_methods.
1992 _quaternion_extension
= {}
1995 def quaternion_extension(cls
,field
):
1997 The quaternion field that we embed/unembed, as an extension
1998 of the given ``field``.
2000 if field
in cls
._quaternion
_extension
:
2001 return cls
._quaternion
_extension
[field
]
2003 Q
= QuaternionAlgebra(field
,-1,-1)
2005 cls
._quaternion
_extension
[field
] = Q
2009 def dimension_over_reals():
2013 def real_embed(cls
,M
):
2015 Embed the n-by-n quaternion matrix ``M`` into the space of real
2016 matrices of size 4n-by-4n by first sending each quaternion entry `z
2017 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
2018 c+di],[-c + di, a-bi]]`, and then embedding those into a real
2023 sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA
2027 sage: Q = QuaternionAlgebra(QQ,-1,-1)
2028 sage: i,j,k = Q.gens()
2029 sage: x = 1 + 2*i + 3*j + 4*k
2030 sage: M = matrix(Q, 1, [[x]])
2031 sage: QuaternionMatrixEJA.real_embed(M)
2037 Embedding is a homomorphism (isomorphism, in fact)::
2039 sage: set_random_seed()
2040 sage: n = ZZ.random_element(2)
2041 sage: Q = QuaternionAlgebra(QQ,-1,-1)
2042 sage: X = random_matrix(Q, n)
2043 sage: Y = random_matrix(Q, n)
2044 sage: Xe = QuaternionMatrixEJA.real_embed(X)
2045 sage: Ye = QuaternionMatrixEJA.real_embed(Y)
2046 sage: XYe = QuaternionMatrixEJA.real_embed(X*Y)
2051 super(QuaternionMatrixEJA
,cls
).real_embed(M
)
2052 quaternions
= M
.base_ring()
2055 F
= QuadraticField(-1, 'I')
2060 t
= z
.coefficient_tuple()
2065 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
2066 [-c
+ d
*i
, a
- b
*i
]])
2067 realM
= ComplexMatrixEJA
.real_embed(cplxM
)
2068 blocks
.append(realM
)
2070 # We should have real entries by now, so use the realest field
2071 # we've got for the return value.
2072 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
2077 def real_unembed(cls
,M
):
2079 The inverse of _embed_quaternion_matrix().
2083 sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA
2087 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
2088 ....: [-2, 1, -4, 3],
2089 ....: [-3, 4, 1, -2],
2090 ....: [-4, -3, 2, 1]])
2091 sage: QuaternionMatrixEJA.real_unembed(M)
2092 [1 + 2*i + 3*j + 4*k]
2096 Unembedding is the inverse of embedding::
2098 sage: set_random_seed()
2099 sage: Q = QuaternionAlgebra(QQ, -1, -1)
2100 sage: M = random_matrix(Q, 3)
2101 sage: Me = QuaternionMatrixEJA.real_embed(M)
2102 sage: QuaternionMatrixEJA.real_unembed(Me) == M
2106 super(QuaternionMatrixEJA
,cls
).real_unembed(M
)
2108 d
= cls
.dimension_over_reals()
2110 # Use the base ring of the matrix to ensure that its entries can be
2111 # multiplied by elements of the quaternion algebra.
2112 Q
= cls
.quaternion_extension(M
.base_ring())
2115 # Go top-left to bottom-right (reading order), converting every
2116 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
2119 for l
in range(n
/d
):
2120 for m
in range(n
/d
):
2121 submat
= ComplexMatrixEJA
.real_unembed(
2122 M
[d
*l
:d
*l
+d
,d
*m
:d
*m
+d
] )
2123 if submat
[0,0] != submat
[1,1].conjugate():
2124 raise ValueError('bad on-diagonal submatrix')
2125 if submat
[0,1] != -submat
[1,0].conjugate():
2126 raise ValueError('bad off-diagonal submatrix')
2127 z
= submat
[0,0].real()
2128 z
+= submat
[0,0].imag()*i
2129 z
+= submat
[0,1].real()*j
2130 z
+= submat
[0,1].imag()*k
2133 return matrix(Q
, n
/d
, elements
)
2136 class QuaternionHermitianEJA(ConcreteEJA
, QuaternionMatrixEJA
):
2138 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
2139 matrices, the usual symmetric Jordan product, and the
2140 real-part-of-trace inner product. It has dimension `2n^2 - n` over
2145 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2149 In theory, our "field" can be any subfield of the reals::
2151 sage: QuaternionHermitianEJA(2, field=RDF)
2152 Euclidean Jordan algebra of dimension 6 over Real Double Field
2153 sage: QuaternionHermitianEJA(2, field=RR)
2154 Euclidean Jordan algebra of dimension 6 over Real Field with
2155 53 bits of precision
2159 The dimension of this algebra is `2*n^2 - n`::
2161 sage: set_random_seed()
2162 sage: n_max = QuaternionHermitianEJA._max_random_instance_size()
2163 sage: n = ZZ.random_element(1, n_max)
2164 sage: J = QuaternionHermitianEJA(n)
2165 sage: J.dimension() == 2*(n^2) - n
2168 The Jordan multiplication is what we think it is::
2170 sage: set_random_seed()
2171 sage: J = QuaternionHermitianEJA.random_instance()
2172 sage: x,y = J.random_elements(2)
2173 sage: actual = (x*y).to_matrix()
2174 sage: X = x.to_matrix()
2175 sage: Y = y.to_matrix()
2176 sage: expected = (X*Y + Y*X)/2
2177 sage: actual == expected
2179 sage: J(expected) == x*y
2182 We can change the generator prefix::
2184 sage: QuaternionHermitianEJA(2, prefix='a').gens()
2185 (a0, a1, a2, a3, a4, a5)
2187 We can construct the (trivial) algebra of rank zero::
2189 sage: QuaternionHermitianEJA(0)
2190 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2194 def _denormalized_basis(cls
, n
):
2196 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
2198 Why do we embed these? Basically, because all of numerical
2199 linear algebra assumes that you're working with vectors consisting
2200 of `n` entries from a field and scalars from the same field. There's
2201 no way to tell SageMath that (for example) the vectors contain
2202 complex numbers, while the scalar field is real.
2206 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2210 sage: set_random_seed()
2211 sage: n = ZZ.random_element(1,5)
2212 sage: B = QuaternionHermitianEJA._denormalized_basis(n)
2213 sage: all( M.is_symmetric() for M in B )
2218 Q
= QuaternionAlgebra(QQ
,-1,-1)
2221 # This is like the symmetric case, but we need to be careful:
2223 # * We want conjugate-symmetry, not just symmetry.
2224 # * The diagonal will (as a result) be real.
2227 Eij
= matrix
.zero(Q
,n
)
2229 for j
in range(i
+1):
2233 Sij
= cls
.real_embed(Eij
)
2236 # The second, third, and fourth ones have a minus
2237 # because they're conjugated.
2238 # Eij = Eij + Eij.transpose()
2240 Sij_real
= cls
.real_embed(Eij
)
2242 # Eij = I*(Eij - Eij.transpose())
2245 Sij_I
= cls
.real_embed(Eij
)
2247 # Eij = J*(Eij - Eij.transpose())
2250 Sij_J
= cls
.real_embed(Eij
)
2252 # Eij = K*(Eij - Eij.transpose())
2255 Sij_K
= cls
.real_embed(Eij
)
2261 # Since we embedded these, we can drop back to the "field" that we
2262 # started with instead of the quaternion algebra "Q".
2263 return tuple( s
.change_ring(field
) for s
in S
)
2266 def __init__(self
, n
, **kwargs
):
2267 # We know this is a valid EJA, but will double-check
2268 # if the user passes check_axioms=True.
2269 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2271 super(QuaternionHermitianEJA
, self
).__init
__(self
._denormalized
_basis
(n
),
2272 self
.jordan_product
,
2273 self
.trace_inner_product
,
2275 # TODO: this could be factored out somehow, but is left here
2276 # because the MatrixEJA is not presently a subclass of the
2277 # FDEJA class that defines rank() and one().
2278 self
.rank
.set_cache(n
)
2279 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
2280 self
.one
.set_cache(self(idV
))
2284 def _max_random_instance_size():
2286 The maximum rank of a random QuaternionHermitianEJA.
2288 return 2 # Dimension 6
2291 def random_instance(cls
, **kwargs
):
2293 Return a random instance of this type of algebra.
2295 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2296 return cls(n
, **kwargs
)
2299 class HadamardEJA(ConcreteEJA
):
2301 Return the Euclidean Jordan Algebra corresponding to the set
2302 `R^n` under the Hadamard product.
2304 Note: this is nothing more than the Cartesian product of ``n``
2305 copies of the spin algebra. Once Cartesian product algebras
2306 are implemented, this can go.
2310 sage: from mjo.eja.eja_algebra import HadamardEJA
2314 This multiplication table can be verified by hand::
2316 sage: J = HadamardEJA(3)
2317 sage: e0,e1,e2 = J.gens()
2333 We can change the generator prefix::
2335 sage: HadamardEJA(3, prefix='r').gens()
2339 def __init__(self
, n
, **kwargs
):
2341 jordan_product
= lambda x
,y
: x
2342 inner_product
= lambda x
,y
: x
2344 def jordan_product(x
,y
):
2346 return P( xi
*yi
for (xi
,yi
) in zip(x
,y
) )
2348 def inner_product(x
,y
):
2351 # New defaults for keyword arguments. Don't orthonormalize
2352 # because our basis is already orthonormal with respect to our
2353 # inner-product. Don't check the axioms, because we know this
2354 # is a valid EJA... but do double-check if the user passes
2355 # check_axioms=True. Note: we DON'T override the "check_field"
2356 # default here, because the user can pass in a field!
2357 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2358 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2360 column_basis
= tuple( b
.column() for b
in FreeModule(ZZ
, n
).basis() )
2361 super().__init
__(column_basis
, jordan_product
, inner_product
, **kwargs
)
2362 self
.rank
.set_cache(n
)
2365 self
.one
.set_cache( self
.zero() )
2367 self
.one
.set_cache( sum(self
.gens()) )
2370 def _max_random_instance_size():
2372 The maximum dimension of a random HadamardEJA.
2377 def random_instance(cls
, **kwargs
):
2379 Return a random instance of this type of algebra.
2381 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2382 return cls(n
, **kwargs
)
2385 class BilinearFormEJA(ConcreteEJA
):
2387 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2388 with the half-trace inner product and jordan product ``x*y =
2389 (<Bx,y>,y_bar>, x0*y_bar + y0*x_bar)`` where `B = 1 \times B22` is
2390 a symmetric positive-definite "bilinear form" matrix. Its
2391 dimension is the size of `B`, and it has rank two in dimensions
2392 larger than two. It reduces to the ``JordanSpinEJA`` when `B` is
2393 the identity matrix of order ``n``.
2395 We insist that the one-by-one upper-left identity block of `B` be
2396 passed in as well so that we can be passed a matrix of size zero
2397 to construct a trivial algebra.
2401 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
2402 ....: JordanSpinEJA)
2406 When no bilinear form is specified, the identity matrix is used,
2407 and the resulting algebra is the Jordan spin algebra::
2409 sage: B = matrix.identity(AA,3)
2410 sage: J0 = BilinearFormEJA(B)
2411 sage: J1 = JordanSpinEJA(3)
2412 sage: J0.multiplication_table() == J0.multiplication_table()
2415 An error is raised if the matrix `B` does not correspond to a
2416 positive-definite bilinear form::
2418 sage: B = matrix.random(QQ,2,3)
2419 sage: J = BilinearFormEJA(B)
2420 Traceback (most recent call last):
2422 ValueError: bilinear form is not positive-definite
2423 sage: B = matrix.zero(QQ,3)
2424 sage: J = BilinearFormEJA(B)
2425 Traceback (most recent call last):
2427 ValueError: bilinear form is not positive-definite
2431 We can create a zero-dimensional algebra::
2433 sage: B = matrix.identity(AA,0)
2434 sage: J = BilinearFormEJA(B)
2438 We can check the multiplication condition given in the Jordan, von
2439 Neumann, and Wigner paper (and also discussed on my "On the
2440 symmetry..." paper). Note that this relies heavily on the standard
2441 choice of basis, as does anything utilizing the bilinear form
2442 matrix. We opt not to orthonormalize the basis, because if we
2443 did, we would have to normalize the `s_{i}` in a similar manner::
2445 sage: set_random_seed()
2446 sage: n = ZZ.random_element(5)
2447 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2448 sage: B11 = matrix.identity(QQ,1)
2449 sage: B22 = M.transpose()*M
2450 sage: B = block_matrix(2,2,[ [B11,0 ],
2452 sage: J = BilinearFormEJA(B, orthonormalize=False)
2453 sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
2454 sage: V = J.vector_space()
2455 sage: sis = [ J( V([0] + (M.inverse()*ei).list()).column() )
2456 ....: for ei in eis ]
2457 sage: actual = [ sis[i]*sis[j]
2458 ....: for i in range(n-1)
2459 ....: for j in range(n-1) ]
2460 sage: expected = [ J.one() if i == j else J.zero()
2461 ....: for i in range(n-1)
2462 ....: for j in range(n-1) ]
2463 sage: actual == expected
2467 def __init__(self
, B
, **kwargs
):
2468 # The matrix "B" is supplied by the user in most cases,
2469 # so it makes sense to check whether or not its positive-
2470 # definite unless we are specifically asked not to...
2471 if ("check_axioms" not in kwargs
) or kwargs
["check_axioms"]:
2472 if not B
.is_positive_definite():
2473 raise ValueError("bilinear form is not positive-definite")
2475 # However, all of the other data for this EJA is computed
2476 # by us in manner that guarantees the axioms are
2477 # satisfied. So, again, unless we are specifically asked to
2478 # verify things, we'll skip the rest of the checks.
2479 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2481 def inner_product(x
,y
):
2482 return (y
.T
*B
*x
)[0,0]
2484 def jordan_product(x
,y
):
2490 z0
= inner_product(y
,x
)
2491 zbar
= y0
*xbar
+ x0
*ybar
2492 return P([z0
] + zbar
.list())
2495 column_basis
= tuple( b
.column() for b
in FreeModule(ZZ
, n
).basis() )
2496 super(BilinearFormEJA
, self
).__init
__(column_basis
,
2501 # The rank of this algebra is two, unless we're in a
2502 # one-dimensional ambient space (because the rank is bounded
2503 # by the ambient dimension).
2504 self
.rank
.set_cache(min(n
,2))
2507 self
.one
.set_cache( self
.zero() )
2509 self
.one
.set_cache( self
.monomial(0) )
2512 def _max_random_instance_size():
2514 The maximum dimension of a random BilinearFormEJA.
2519 def random_instance(cls
, **kwargs
):
2521 Return a random instance of this algebra.
2523 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2525 B
= matrix
.identity(ZZ
, n
)
2526 return cls(B
, **kwargs
)
2528 B11
= matrix
.identity(ZZ
, 1)
2529 M
= matrix
.random(ZZ
, n
-1)
2530 I
= matrix
.identity(ZZ
, n
-1)
2532 while alpha
.is_zero():
2533 alpha
= ZZ
.random_element().abs()
2534 B22
= M
.transpose()*M
+ alpha
*I
2536 from sage
.matrix
.special
import block_matrix
2537 B
= block_matrix(2,2, [ [B11
, ZZ(0) ],
2540 return cls(B
, **kwargs
)
2543 class JordanSpinEJA(BilinearFormEJA
):
2545 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2546 with the usual inner product and jordan product ``x*y =
2547 (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2552 sage: from mjo.eja.eja_algebra import JordanSpinEJA
2556 This multiplication table can be verified by hand::
2558 sage: J = JordanSpinEJA(4)
2559 sage: e0,e1,e2,e3 = J.gens()
2575 We can change the generator prefix::
2577 sage: JordanSpinEJA(2, prefix='B').gens()
2582 Ensure that we have the usual inner product on `R^n`::
2584 sage: set_random_seed()
2585 sage: J = JordanSpinEJA.random_instance()
2586 sage: x,y = J.random_elements(2)
2587 sage: actual = x.inner_product(y)
2588 sage: expected = x.to_vector().inner_product(y.to_vector())
2589 sage: actual == expected
2593 def __init__(self
, n
, **kwargs
):
2594 # This is a special case of the BilinearFormEJA with the
2595 # identity matrix as its bilinear form.
2596 B
= matrix
.identity(ZZ
, n
)
2598 # Don't orthonormalize because our basis is already
2599 # orthonormal with respect to our inner-product.
2600 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2602 # But also don't pass check_field=False here, because the user
2603 # can pass in a field!
2604 super(JordanSpinEJA
, self
).__init
__(B
, **kwargs
)
2607 def _max_random_instance_size():
2609 The maximum dimension of a random JordanSpinEJA.
2614 def random_instance(cls
, **kwargs
):
2616 Return a random instance of this type of algebra.
2618 Needed here to override the implementation for ``BilinearFormEJA``.
2620 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2621 return cls(n
, **kwargs
)
2624 class TrivialEJA(ConcreteEJA
):
2626 The trivial Euclidean Jordan algebra consisting of only a zero element.
2630 sage: from mjo.eja.eja_algebra import TrivialEJA
2634 sage: J = TrivialEJA()
2641 sage: 7*J.one()*12*J.one()
2643 sage: J.one().inner_product(J.one())
2645 sage: J.one().norm()
2647 sage: J.one().subalgebra_generated_by()
2648 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2653 def __init__(self
, **kwargs
):
2654 jordan_product
= lambda x
,y
: x
2655 inner_product
= lambda x
,y
: 0
2658 # New defaults for keyword arguments
2659 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2660 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2662 super(TrivialEJA
, self
).__init
__(basis
,
2666 # The rank is zero using my definition, namely the dimension of the
2667 # largest subalgebra generated by any element.
2668 self
.rank
.set_cache(0)
2669 self
.one
.set_cache( self
.zero() )
2672 def random_instance(cls
, **kwargs
):
2673 # We don't take a "size" argument so the superclass method is
2674 # inappropriate for us.
2675 return cls(**kwargs
)
2678 class DirectSumEJA(FiniteDimensionalEJA
):
2680 The external (orthogonal) direct sum of two other Euclidean Jordan
2681 algebras. Essentially the Cartesian product of its two factors.
2682 Every Euclidean Jordan algebra decomposes into an orthogonal
2683 direct sum of simple Euclidean Jordan algebras, so no generality
2684 is lost by providing only this construction.
2688 sage: from mjo.eja.eja_algebra import (random_eja,
2690 ....: RealSymmetricEJA,
2695 sage: J1 = HadamardEJA(2)
2696 sage: J2 = RealSymmetricEJA(3)
2697 sage: J = DirectSumEJA(J1,J2)
2702 sage: J.matrix_space()
2703 The Cartesian product of (Full MatrixSpace of 2 by 1 dense matrices
2704 over Algebraic Real Field, Full MatrixSpace of 3 by 3 dense matrices
2705 over Algebraic Real Field)
2709 The external direct sum construction is only valid when the two factors
2710 have the same base ring; an error is raised otherwise::
2712 sage: set_random_seed()
2713 sage: J1 = random_eja(field=AA)
2714 sage: J2 = random_eja(field=QQ,orthonormalize=False)
2715 sage: J = DirectSumEJA(J1,J2)
2716 Traceback (most recent call last):
2718 ValueError: algebras must share the same base field
2721 def __init__(self
, J1
, J2
, **kwargs
):
2722 if J1
.base_ring() != J2
.base_ring():
2723 raise ValueError("algebras must share the same base field")
2724 field
= J1
.base_ring()
2726 M
= J1
.matrix_space().cartesian_product(J2
.matrix_space())
2727 self
._cartprod
_algebra
= J1
.cartesian_product(J2
)
2729 self
._matrix
_basis
= tuple( [M((a
,0)) for a
in J1
.matrix_basis()] +
2730 [M((0,b
)) for b
in J2
.matrix_basis()] )
2732 n
= len(self
._matrix
_basis
)
2734 CombinatorialFreeModule
.__init
__(
2738 category
=self
._cartprod
_algebra
.category(),
2741 self
.rank
.set_cache(J1
.rank() + J2
.rank())
2745 def product(self
,x
,y
):
2749 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
2750 ....: ComplexHermitianEJA,
2755 sage: set_random_seed()
2756 sage: J1 = JordanSpinEJA(3, field=QQ)
2757 sage: J2 = ComplexHermitianEJA(2, field=QQ, orthonormalize=False)
2758 sage: J = DirectSumEJA(J1,J2)
2759 sage: J.random_element()*J.random_element() in J
2763 xv
= self
._cartprod
_algebra
.from_vector(x
.to_vector())
2764 yv
= self
._cartprod
_algebra
.from_vector(y
.to_vector())
2765 return self
.from_vector((xv
*yv
).to_vector())
2768 def cartesian_factors(self
):
2770 Return the pair of this algebra's factors.
2774 sage: from mjo.eja.eja_algebra import (HadamardEJA,
2775 ....: JordanSpinEJA,
2780 sage: J1 = HadamardEJA(2, field=QQ)
2781 sage: J2 = JordanSpinEJA(3, field=QQ)
2782 sage: J = DirectSumEJA(J1,J2)
2783 sage: J.cartesian_factors()
2784 (Euclidean Jordan algebra of dimension 2 over Rational Field,
2785 Euclidean Jordan algebra of dimension 3 over Rational Field)
2788 return self
._cartprod
_algebra
.cartesian_factors()
2791 # def projections(self):
2793 # Return a pair of projections onto this algebra's factors.
2797 # sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
2798 # ....: ComplexHermitianEJA,
2799 # ....: DirectSumEJA)
2803 # sage: J1 = JordanSpinEJA(2)
2804 # sage: J2 = ComplexHermitianEJA(2)
2805 # sage: J = DirectSumEJA(J1,J2)
2806 # sage: (pi_left, pi_right) = J.projections()
2807 # sage: J.one().to_vector()
2808 # (1, 0, 1, 0, 0, 1)
2809 # sage: pi_left(J.one()).to_vector()
2811 # sage: pi_right(J.one()).to_vector()
2815 # (J1,J2) = self.factors()
2816 # m = J1.dimension()
2817 # n = J2.dimension()
2818 # V_basis = self.vector_space().basis()
2819 # # Need to specify the dimensions explicitly so that we don't
2820 # # wind up with a zero-by-zero matrix when we want e.g. a
2821 # # zero-by-two matrix (important for composing things).
2822 # P1 = matrix(self.base_ring(), m, m+n, V_basis[:m])
2823 # P2 = matrix(self.base_ring(), n, m+n, V_basis[m:])
2824 # pi_left = FiniteDimensionalEJAOperator(self,J1,P1)
2825 # pi_right = FiniteDimensionalEJAOperator(self,J2,P2)
2826 # return (pi_left, pi_right)
2828 # def inclusions(self):
2830 # Return the pair of inclusion maps from our factors into us.
2834 # sage: from mjo.eja.eja_algebra import (random_eja,
2835 # ....: JordanSpinEJA,
2836 # ....: RealSymmetricEJA,
2837 # ....: DirectSumEJA)
2841 # sage: J1 = JordanSpinEJA(3)
2842 # sage: J2 = RealSymmetricEJA(2)
2843 # sage: J = DirectSumEJA(J1,J2)
2844 # sage: (iota_left, iota_right) = J.inclusions()
2845 # sage: iota_left(J1.zero()) == J.zero()
2847 # sage: iota_right(J2.zero()) == J.zero()
2849 # sage: J1.one().to_vector()
2851 # sage: iota_left(J1.one()).to_vector()
2852 # (1, 0, 0, 0, 0, 0)
2853 # sage: J2.one().to_vector()
2855 # sage: iota_right(J2.one()).to_vector()
2856 # (0, 0, 0, 1, 0, 1)
2857 # sage: J.one().to_vector()
2858 # (1, 0, 0, 1, 0, 1)
2862 # Composing a projection with the corresponding inclusion should
2863 # produce the identity map, and mismatching them should produce
2866 # sage: set_random_seed()
2867 # sage: J1 = random_eja()
2868 # sage: J2 = random_eja()
2869 # sage: J = DirectSumEJA(J1,J2)
2870 # sage: (iota_left, iota_right) = J.inclusions()
2871 # sage: (pi_left, pi_right) = J.projections()
2872 # sage: pi_left*iota_left == J1.one().operator()
2874 # sage: pi_right*iota_right == J2.one().operator()
2876 # sage: (pi_left*iota_right).is_zero()
2878 # sage: (pi_right*iota_left).is_zero()
2882 # (J1,J2) = self.factors()
2883 # m = J1.dimension()
2884 # n = J2.dimension()
2885 # V_basis = self.vector_space().basis()
2886 # # Need to specify the dimensions explicitly so that we don't
2887 # # wind up with a zero-by-zero matrix when we want e.g. a
2888 # # two-by-zero matrix (important for composing things).
2889 # I1 = matrix.column(self.base_ring(), m, m+n, V_basis[:m])
2890 # I2 = matrix.column(self.base_ring(), n, m+n, V_basis[m:])
2891 # iota_left = FiniteDimensionalEJAOperator(J1,self,I1)
2892 # iota_right = FiniteDimensionalEJAOperator(J2,self,I2)
2893 # return (iota_left, iota_right)
2895 # def inner_product(self, x, y):
2897 # The standard Cartesian inner-product.
2899 # We project ``x`` and ``y`` onto our factors, and add up the
2900 # inner-products from the subalgebras.
2905 # sage: from mjo.eja.eja_algebra import (HadamardEJA,
2906 # ....: QuaternionHermitianEJA,
2907 # ....: DirectSumEJA)
2911 # sage: J1 = HadamardEJA(3,field=QQ)
2912 # sage: J2 = QuaternionHermitianEJA(2,field=QQ,orthonormalize=False)
2913 # sage: J = DirectSumEJA(J1,J2)
2914 # sage: x1 = J1.one()
2916 # sage: y1 = J2.one()
2918 # sage: x1.inner_product(x2)
2920 # sage: y1.inner_product(y2)
2922 # sage: J.one().inner_product(J.one())
2926 # (pi_left, pi_right) = self.projections()
2932 # return (x1.inner_product(y1) + x2.inner_product(y2))
2936 random_eja
= ConcreteEJA
.random_instance