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
71 if not field
.is_subring(RR
):
72 # Note: this does return true for the real algebraic
73 # field, the rationals, and any quadratic field where
74 # we've specified a real embedding.
75 raise ValueError("scalar field is not real")
77 # If the basis given to us wasn't over the field that it's
78 # supposed to be over, fix that. Or, you know, crash.
79 basis
= tuple( b
.change_ring(field
) for b
in basis
)
82 # Check commutativity of the Jordan and inner-products.
83 # This has to be done before we build the multiplication
84 # and inner-product tables/matrices, because we take
85 # advantage of symmetry in the process.
86 if not all( jordan_product(bi
,bj
) == jordan_product(bj
,bi
)
89 raise ValueError("Jordan product is not commutative")
91 if not all( inner_product(bi
,bj
) == inner_product(bj
,bi
)
94 raise ValueError("inner-product is not commutative")
98 category
= MagmaticAlgebras(field
).FiniteDimensional()
99 category
= category
.WithBasis().Unital()
101 # Element subalgebras can take advantage of this.
102 category
= category
.Associative()
104 # Call the superclass constructor so that we can use its from_vector()
105 # method to build our multiplication table.
107 super().__init
__(field
,
113 # Now comes all of the hard work. We'll be constructing an
114 # ambient vector space V that our (vectorized) basis lives in,
115 # as well as a subspace W of V spanned by those (vectorized)
116 # basis elements. The W-coordinates are the coefficients that
117 # we see in things like x = 1*e1 + 2*e2.
122 # Works on both column and square matrices...
123 degree
= len(basis
[0].list())
125 # Build an ambient space that fits our matrix basis when
126 # written out as "long vectors."
127 V
= VectorSpace(field
, degree
)
129 # The matrix that will hole the orthonormal -> unorthonormal
130 # coordinate transformation.
131 self
._deortho
_matrix
= None
134 # Save a copy of the un-orthonormalized basis for later.
135 # Convert it to ambient V (vector) coordinates while we're
136 # at it, because we'd have to do it later anyway.
137 deortho_vector_basis
= tuple( V(b
.list()) for b
in basis
)
139 from mjo
.eja
.eja_utils
import gram_schmidt
140 basis
= tuple(gram_schmidt(basis
, inner_product
))
142 # Save the (possibly orthonormalized) matrix basis for
144 self
._matrix
_basis
= basis
146 # Now create the vector space for the algebra, which will have
147 # its own set of non-ambient coordinates (in terms of the
149 vector_basis
= tuple( V(b
.list()) for b
in basis
)
150 W
= V
.span_of_basis( vector_basis
, check
=check_axioms
)
153 # Now "W" is the vector space of our algebra coordinates. The
154 # variables "X1", "X2",... refer to the entries of vectors in
155 # W. Thus to convert back and forth between the orthonormal
156 # coordinates and the given ones, we need to stick the original
158 U
= V
.span_of_basis( deortho_vector_basis
, check
=check_axioms
)
159 self
._deortho
_matrix
= matrix( U
.coordinate_vector(q
)
160 for q
in vector_basis
)
163 # Now we actually compute the multiplication and inner-product
164 # tables/matrices using the possibly-orthonormalized basis.
165 self
._inner
_product
_matrix
= matrix
.identity(field
, n
)
166 self
._multiplication
_table
= [ [0 for j
in range(i
+1)]
169 # Note: the Jordan and inner-products are defined in terms
170 # of the ambient basis. It's important that their arguments
171 # are in ambient coordinates as well.
174 # ortho basis w.r.t. ambient coords
178 # The jordan product returns a matrixy answer, so we
179 # have to convert it to the algebra coordinates.
180 elt
= jordan_product(q_i
, q_j
)
181 elt
= W
.coordinate_vector(V(elt
.list()))
182 self
._multiplication
_table
[i
][j
] = self
.from_vector(elt
)
184 if not orthonormalize
:
185 # If we're orthonormalizing the basis with respect
186 # to an inner-product, then the inner-product
187 # matrix with respect to the resulting basis is
188 # just going to be the identity.
189 ip
= inner_product(q_i
, q_j
)
190 self
._inner
_product
_matrix
[i
,j
] = ip
191 self
._inner
_product
_matrix
[j
,i
] = ip
193 self
._inner
_product
_matrix
._cache
= {'hermitian': True}
194 self
._inner
_product
_matrix
.set_immutable()
197 if not self
._is
_jordanian
():
198 raise ValueError("Jordan identity does not hold")
199 if not self
._inner
_product
_is
_associative
():
200 raise ValueError("inner product is not associative")
203 def _coerce_map_from_base_ring(self
):
205 Disable the map from the base ring into the algebra.
207 Performing a nonsense conversion like this automatically
208 is counterpedagogical. The fallback is to try the usual
209 element constructor, which should also fail.
213 sage: from mjo.eja.eja_algebra import random_eja
217 sage: set_random_seed()
218 sage: J = random_eja()
220 Traceback (most recent call last):
222 ValueError: not an element of this algebra
228 def product_on_basis(self
, i
, j
):
229 # We only stored the lower-triangular portion of the
230 # multiplication table.
232 return self
._multiplication
_table
[i
][j
]
234 return self
._multiplication
_table
[j
][i
]
236 def inner_product(self
, x
, y
):
238 The inner product associated with this Euclidean Jordan algebra.
240 Defaults to the trace inner product, but can be overridden by
241 subclasses if they are sure that the necessary properties are
246 sage: from mjo.eja.eja_algebra import (random_eja,
248 ....: BilinearFormEJA)
252 Our inner product is "associative," which means the following for
253 a symmetric bilinear form::
255 sage: set_random_seed()
256 sage: J = random_eja()
257 sage: x,y,z = J.random_elements(3)
258 sage: (x*y).inner_product(z) == y.inner_product(x*z)
263 Ensure that this is the usual inner product for the algebras
266 sage: set_random_seed()
267 sage: J = HadamardEJA.random_instance()
268 sage: x,y = J.random_elements(2)
269 sage: actual = x.inner_product(y)
270 sage: expected = x.to_vector().inner_product(y.to_vector())
271 sage: actual == expected
274 Ensure that this is one-half of the trace inner-product in a
275 BilinearFormEJA that isn't just the reals (when ``n`` isn't
276 one). This is in Faraut and Koranyi, and also my "On the
279 sage: set_random_seed()
280 sage: J = BilinearFormEJA.random_instance()
281 sage: n = J.dimension()
282 sage: x = J.random_element()
283 sage: y = J.random_element()
284 sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
287 B
= self
._inner
_product
_matrix
288 return (B
*x
.to_vector()).inner_product(y
.to_vector())
291 def _is_commutative(self
):
293 Whether or not this algebra's multiplication table is commutative.
295 This method should of course always return ``True``, unless
296 this algebra was constructed with ``check_axioms=False`` and
297 passed an invalid multiplication table.
299 return all( self
.product_on_basis(i
,j
) == self
.product_on_basis(i
,j
)
300 for i
in range(self
.dimension())
301 for j
in range(self
.dimension()) )
303 def _is_jordanian(self
):
305 Whether or not this algebra's multiplication table respects the
306 Jordan identity `(x^{2})(xy) = x(x^{2}y)`.
308 We only check one arrangement of `x` and `y`, so for a
309 ``True`` result to be truly true, you should also check
310 :meth:`_is_commutative`. This method should of course always
311 return ``True``, unless this algebra was constructed with
312 ``check_axioms=False`` and passed an invalid multiplication table.
314 return all( (self
.monomial(i
)**2)*(self
.monomial(i
)*self
.monomial(j
))
316 (self
.monomial(i
))*((self
.monomial(i
)**2)*self
.monomial(j
))
317 for i
in range(self
.dimension())
318 for j
in range(self
.dimension()) )
320 def _inner_product_is_associative(self
):
322 Return whether or not this algebra's inner product `B` is
323 associative; that is, whether or not `B(xy,z) = B(x,yz)`.
325 This method should of course always return ``True``, unless
326 this algebra was constructed with ``check_axioms=False`` and
327 passed an invalid Jordan or inner-product.
330 # Used to check whether or not something is zero in an inexact
331 # ring. This number is sufficient to allow the construction of
332 # QuaternionHermitianEJA(2, field=RDF) with check_axioms=True.
335 for i
in range(self
.dimension()):
336 for j
in range(self
.dimension()):
337 for k
in range(self
.dimension()):
341 diff
= (x
*y
).inner_product(z
) - x
.inner_product(y
*z
)
343 if self
.base_ring().is_exact():
347 if diff
.abs() > epsilon
:
352 def _element_constructor_(self
, elt
):
354 Construct an element of this algebra from its vector or matrix
357 This gets called only after the parent element _call_ method
358 fails to find a coercion for the argument.
362 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
364 ....: RealSymmetricEJA)
368 The identity in `S^n` is converted to the identity in the EJA::
370 sage: J = RealSymmetricEJA(3)
371 sage: I = matrix.identity(QQ,3)
372 sage: J(I) == J.one()
375 This skew-symmetric matrix can't be represented in the EJA::
377 sage: J = RealSymmetricEJA(3)
378 sage: A = matrix(QQ,3, lambda i,j: i-j)
380 Traceback (most recent call last):
382 ValueError: not an element of this algebra
386 Ensure that we can convert any element of the two non-matrix
387 simple algebras (whose matrix representations are columns)
388 back and forth faithfully::
390 sage: set_random_seed()
391 sage: J = HadamardEJA.random_instance()
392 sage: x = J.random_element()
393 sage: J(x.to_vector().column()) == x
395 sage: J = JordanSpinEJA.random_instance()
396 sage: x = J.random_element()
397 sage: J(x.to_vector().column()) == x
401 msg
= "not an element of this algebra"
403 # The superclass implementation of random_element()
404 # needs to be able to coerce "0" into the algebra.
406 elif elt
in self
.base_ring():
407 # Ensure that no base ring -> algebra coercion is performed
408 # by this method. There's some stupidity in sage that would
409 # otherwise propagate to this method; for example, sage thinks
410 # that the integer 3 belongs to the space of 2-by-2 matrices.
411 raise ValueError(msg
)
415 except (AttributeError, TypeError):
416 # Try to convert a vector into a column-matrix
419 if elt
not in self
.matrix_space():
420 raise ValueError(msg
)
422 # Thanks for nothing! Matrix spaces aren't vector spaces in
423 # Sage, so we have to figure out its matrix-basis coordinates
424 # ourselves. We use the basis space's ring instead of the
425 # element's ring because the basis space might be an algebraic
426 # closure whereas the base ring of the 3-by-3 identity matrix
427 # could be QQ instead of QQbar.
429 # We pass check=False because the matrix basis is "guaranteed"
430 # to be linearly independent... right? Ha ha.
431 V
= VectorSpace(self
.base_ring(), elt
.nrows()*elt
.ncols())
432 W
= V
.span_of_basis( (_mat2vec(s
) for s
in self
.matrix_basis()),
436 coords
= W
.coordinate_vector(_mat2vec(elt
))
437 except ArithmeticError: # vector is not in free module
438 raise ValueError(msg
)
440 return self
.from_vector(coords
)
444 Return a string representation of ``self``.
448 sage: from mjo.eja.eja_algebra import JordanSpinEJA
452 Ensure that it says what we think it says::
454 sage: JordanSpinEJA(2, field=AA)
455 Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
456 sage: JordanSpinEJA(3, field=RDF)
457 Euclidean Jordan algebra of dimension 3 over Real Double Field
460 fmt
= "Euclidean Jordan algebra of dimension {} over {}"
461 return fmt
.format(self
.dimension(), self
.base_ring())
465 def characteristic_polynomial_of(self
):
467 Return the algebra's "characteristic polynomial of" function,
468 which is itself a multivariate polynomial that, when evaluated
469 at the coordinates of some algebra element, returns that
470 element's characteristic polynomial.
472 The resulting polynomial has `n+1` variables, where `n` is the
473 dimension of this algebra. The first `n` variables correspond to
474 the coordinates of an algebra element: when evaluated at the
475 coordinates of an algebra element with respect to a certain
476 basis, the result is a univariate polynomial (in the one
477 remaining variable ``t``), namely the characteristic polynomial
482 sage: from mjo.eja.eja_algebra import JordanSpinEJA, TrivialEJA
486 The characteristic polynomial in the spin algebra is given in
487 Alizadeh, Example 11.11::
489 sage: J = JordanSpinEJA(3)
490 sage: p = J.characteristic_polynomial_of(); p
491 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
492 sage: xvec = J.one().to_vector()
496 By definition, the characteristic polynomial is a monic
497 degree-zero polynomial in a rank-zero algebra. Note that
498 Cayley-Hamilton is indeed satisfied since the polynomial
499 ``1`` evaluates to the identity element of the algebra on
502 sage: J = TrivialEJA()
503 sage: J.characteristic_polynomial_of()
510 # The list of coefficient polynomials a_0, a_1, a_2, ..., a_(r-1).
511 a
= self
._charpoly
_coefficients
()
513 # We go to a bit of trouble here to reorder the
514 # indeterminates, so that it's easier to evaluate the
515 # characteristic polynomial at x's coordinates and get back
516 # something in terms of t, which is what we want.
517 S
= PolynomialRing(self
.base_ring(),'t')
521 S
= PolynomialRing(S
, R
.variable_names())
524 return (t
**r
+ sum( a
[k
]*(t
**k
) for k
in range(r
) ))
526 def coordinate_polynomial_ring(self
):
528 The multivariate polynomial ring in which this algebra's
529 :meth:`characteristic_polynomial_of` lives.
533 sage: from mjo.eja.eja_algebra import (HadamardEJA,
534 ....: RealSymmetricEJA)
538 sage: J = HadamardEJA(2)
539 sage: J.coordinate_polynomial_ring()
540 Multivariate Polynomial Ring in X1, X2...
541 sage: J = RealSymmetricEJA(3,field=QQ,orthonormalize=False)
542 sage: J.coordinate_polynomial_ring()
543 Multivariate Polynomial Ring in X1, X2, X3, X4, X5, X6...
546 var_names
= tuple( "X%d" % z
for z
in range(1, self
.dimension()+1) )
547 return PolynomialRing(self
.base_ring(), var_names
)
549 def inner_product(self
, x
, y
):
551 The inner product associated with this Euclidean Jordan algebra.
553 Defaults to the trace inner product, but can be overridden by
554 subclasses if they are sure that the necessary properties are
559 sage: from mjo.eja.eja_algebra import (random_eja,
561 ....: BilinearFormEJA)
565 Our inner product is "associative," which means the following for
566 a symmetric bilinear form::
568 sage: set_random_seed()
569 sage: J = random_eja()
570 sage: x,y,z = J.random_elements(3)
571 sage: (x*y).inner_product(z) == y.inner_product(x*z)
576 Ensure that this is the usual inner product for the algebras
579 sage: set_random_seed()
580 sage: J = HadamardEJA.random_instance()
581 sage: x,y = J.random_elements(2)
582 sage: actual = x.inner_product(y)
583 sage: expected = x.to_vector().inner_product(y.to_vector())
584 sage: actual == expected
587 Ensure that this is one-half of the trace inner-product in a
588 BilinearFormEJA that isn't just the reals (when ``n`` isn't
589 one). This is in Faraut and Koranyi, and also my "On the
592 sage: set_random_seed()
593 sage: J = BilinearFormEJA.random_instance()
594 sage: n = J.dimension()
595 sage: x = J.random_element()
596 sage: y = J.random_element()
597 sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
600 B
= self
._inner
_product
_matrix
601 return (B
*x
.to_vector()).inner_product(y
.to_vector())
604 def is_trivial(self
):
606 Return whether or not this algebra is trivial.
608 A trivial algebra contains only the zero element.
612 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
617 sage: J = ComplexHermitianEJA(3)
623 sage: J = TrivialEJA()
628 return self
.dimension() == 0
631 def multiplication_table(self
):
633 Return a visual representation of this algebra's multiplication
634 table (on basis elements).
638 sage: from mjo.eja.eja_algebra import JordanSpinEJA
642 sage: J = JordanSpinEJA(4)
643 sage: J.multiplication_table()
644 +----++----+----+----+----+
645 | * || e0 | e1 | e2 | e3 |
646 +====++====+====+====+====+
647 | e0 || e0 | e1 | e2 | e3 |
648 +----++----+----+----+----+
649 | e1 || e1 | e0 | 0 | 0 |
650 +----++----+----+----+----+
651 | e2 || e2 | 0 | e0 | 0 |
652 +----++----+----+----+----+
653 | e3 || e3 | 0 | 0 | e0 |
654 +----++----+----+----+----+
658 # Prepend the header row.
659 M
= [["*"] + list(self
.gens())]
661 # And to each subsequent row, prepend an entry that belongs to
662 # the left-side "header column."
663 M
+= [ [self
.monomial(i
)] + [ self
.product_on_basis(i
,j
)
667 return table(M
, header_row
=True, header_column
=True, frame
=True)
670 def matrix_basis(self
):
672 Return an (often more natural) representation of this algebras
673 basis as an ordered tuple of matrices.
675 Every finite-dimensional Euclidean Jordan Algebra is a, up to
676 Jordan isomorphism, a direct sum of five simple
677 algebras---four of which comprise Hermitian matrices. And the
678 last type of algebra can of course be thought of as `n`-by-`1`
679 column matrices (ambiguusly called column vectors) to avoid
680 special cases. As a result, matrices (and column vectors) are
681 a natural representation format for Euclidean Jordan algebra
684 But, when we construct an algebra from a basis of matrices,
685 those matrix representations are lost in favor of coordinate
686 vectors *with respect to* that basis. We could eventually
687 convert back if we tried hard enough, but having the original
688 representations handy is valuable enough that we simply store
689 them and return them from this method.
691 Why implement this for non-matrix algebras? Avoiding special
692 cases for the :class:`BilinearFormEJA` pays with simplicity in
693 its own right. But mainly, we would like to be able to assume
694 that elements of a :class:`CartesianProductEJA` can be displayed
695 nicely, without having to have special classes for direct sums
696 one of whose components was a matrix algebra.
700 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
701 ....: RealSymmetricEJA)
705 sage: J = RealSymmetricEJA(2)
707 Finite family {0: e0, 1: e1, 2: e2}
708 sage: J.matrix_basis()
710 [1 0] [ 0 0.7071067811865475?] [0 0]
711 [0 0], [0.7071067811865475? 0], [0 1]
716 sage: J = JordanSpinEJA(2)
718 Finite family {0: e0, 1: e1}
719 sage: J.matrix_basis()
725 return self
._matrix
_basis
728 def matrix_space(self
):
730 Return the matrix space in which this algebra's elements live, if
731 we think of them as matrices (including column vectors of the
734 Generally this will be an `n`-by-`1` column-vector space,
735 except when the algebra is trivial. There it's `n`-by-`n`
736 (where `n` is zero), to ensure that two elements of the matrix
737 space (empty matrices) can be multiplied.
739 Matrix algebras override this with something more useful.
741 if self
.is_trivial():
742 return MatrixSpace(self
.base_ring(), 0)
744 return self
.matrix_basis()[0].parent()
750 Return the unit element of this algebra.
754 sage: from mjo.eja.eja_algebra import (HadamardEJA,
759 We can compute unit element in the Hadamard EJA::
761 sage: J = HadamardEJA(5)
763 e0 + e1 + e2 + e3 + e4
765 The unit element in the Hadamard EJA is inherited in the
766 subalgebras generated by its elements::
768 sage: J = HadamardEJA(5)
770 e0 + e1 + e2 + e3 + e4
771 sage: x = sum(J.gens())
772 sage: A = x.subalgebra_generated_by(orthonormalize=False)
775 sage: A.one().superalgebra_element()
776 e0 + e1 + e2 + e3 + e4
780 The identity element acts like the identity, regardless of
781 whether or not we orthonormalize::
783 sage: set_random_seed()
784 sage: J = random_eja()
785 sage: x = J.random_element()
786 sage: J.one()*x == x and x*J.one() == x
788 sage: A = x.subalgebra_generated_by()
789 sage: y = A.random_element()
790 sage: A.one()*y == y and y*A.one() == y
795 sage: set_random_seed()
796 sage: J = random_eja(field=QQ, orthonormalize=False)
797 sage: x = J.random_element()
798 sage: J.one()*x == x and x*J.one() == x
800 sage: A = x.subalgebra_generated_by(orthonormalize=False)
801 sage: y = A.random_element()
802 sage: A.one()*y == y and y*A.one() == y
805 The matrix of the unit element's operator is the identity,
806 regardless of the base field and whether or not we
809 sage: set_random_seed()
810 sage: J = random_eja()
811 sage: actual = J.one().operator().matrix()
812 sage: expected = matrix.identity(J.base_ring(), J.dimension())
813 sage: actual == expected
815 sage: x = J.random_element()
816 sage: A = x.subalgebra_generated_by()
817 sage: actual = A.one().operator().matrix()
818 sage: expected = matrix.identity(A.base_ring(), A.dimension())
819 sage: actual == expected
824 sage: set_random_seed()
825 sage: J = random_eja(field=QQ, orthonormalize=False)
826 sage: actual = J.one().operator().matrix()
827 sage: expected = matrix.identity(J.base_ring(), J.dimension())
828 sage: actual == expected
830 sage: x = J.random_element()
831 sage: A = x.subalgebra_generated_by(orthonormalize=False)
832 sage: actual = A.one().operator().matrix()
833 sage: expected = matrix.identity(A.base_ring(), A.dimension())
834 sage: actual == expected
837 Ensure that the cached unit element (often precomputed by
838 hand) agrees with the computed one::
840 sage: set_random_seed()
841 sage: J = random_eja()
842 sage: cached = J.one()
843 sage: J.one.clear_cache()
844 sage: J.one() == cached
849 sage: set_random_seed()
850 sage: J = random_eja(field=QQ, orthonormalize=False)
851 sage: cached = J.one()
852 sage: J.one.clear_cache()
853 sage: J.one() == cached
857 # We can brute-force compute the matrices of the operators
858 # that correspond to the basis elements of this algebra.
859 # If some linear combination of those basis elements is the
860 # algebra identity, then the same linear combination of
861 # their matrices has to be the identity matrix.
863 # Of course, matrices aren't vectors in sage, so we have to
864 # appeal to the "long vectors" isometry.
865 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
867 # Now we use basic linear algebra to find the coefficients,
868 # of the matrices-as-vectors-linear-combination, which should
869 # work for the original algebra basis too.
870 A
= matrix(self
.base_ring(), oper_vecs
)
872 # We used the isometry on the left-hand side already, but we
873 # still need to do it for the right-hand side. Recall that we
874 # wanted something that summed to the identity matrix.
875 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
877 # Now if there's an identity element in the algebra, this
878 # should work. We solve on the left to avoid having to
879 # transpose the matrix "A".
880 return self
.from_vector(A
.solve_left(b
))
883 def peirce_decomposition(self
, c
):
885 The Peirce decomposition of this algebra relative to the
888 In the future, this can be extended to a complete system of
889 orthogonal idempotents.
893 - ``c`` -- an idempotent of this algebra.
897 A triple (J0, J5, J1) containing two subalgebras and one subspace
900 - ``J0`` -- the algebra on the eigenspace of ``c.operator()``
901 corresponding to the eigenvalue zero.
903 - ``J5`` -- the eigenspace (NOT a subalgebra) of ``c.operator()``
904 corresponding to the eigenvalue one-half.
906 - ``J1`` -- the algebra on the eigenspace of ``c.operator()``
907 corresponding to the eigenvalue one.
909 These are the only possible eigenspaces for that operator, and this
910 algebra is a direct sum of them. The spaces ``J0`` and ``J1`` are
911 orthogonal, and are subalgebras of this algebra with the appropriate
916 sage: from mjo.eja.eja_algebra import random_eja, RealSymmetricEJA
920 The canonical example comes from the symmetric matrices, which
921 decompose into diagonal and off-diagonal parts::
923 sage: J = RealSymmetricEJA(3)
924 sage: C = matrix(QQ, [ [1,0,0],
928 sage: J0,J5,J1 = J.peirce_decomposition(c)
930 Euclidean Jordan algebra of dimension 1...
932 Vector space of degree 6 and dimension 2...
934 Euclidean Jordan algebra of dimension 3...
935 sage: J0.one().to_matrix()
939 sage: orig_df = AA.options.display_format
940 sage: AA.options.display_format = 'radical'
941 sage: J.from_vector(J5.basis()[0]).to_matrix()
945 sage: J.from_vector(J5.basis()[1]).to_matrix()
949 sage: AA.options.display_format = orig_df
950 sage: J1.one().to_matrix()
957 Every algebra decomposes trivially with respect to its identity
960 sage: set_random_seed()
961 sage: J = random_eja()
962 sage: J0,J5,J1 = J.peirce_decomposition(J.one())
963 sage: J0.dimension() == 0 and J5.dimension() == 0
965 sage: J1.superalgebra() == J and J1.dimension() == J.dimension()
968 The decomposition is into eigenspaces, and its components are
969 therefore necessarily orthogonal. Moreover, the identity
970 elements in the two subalgebras are the projections onto their
971 respective subspaces of the superalgebra's identity element::
973 sage: set_random_seed()
974 sage: J = random_eja()
975 sage: x = J.random_element()
976 sage: if not J.is_trivial():
977 ....: while x.is_nilpotent():
978 ....: x = J.random_element()
979 sage: c = x.subalgebra_idempotent()
980 sage: J0,J5,J1 = J.peirce_decomposition(c)
982 sage: for (w,y,z) in zip(J0.basis(), J5.basis(), J1.basis()):
983 ....: w = w.superalgebra_element()
984 ....: y = J.from_vector(y)
985 ....: z = z.superalgebra_element()
986 ....: ipsum += w.inner_product(y).abs()
987 ....: ipsum += w.inner_product(z).abs()
988 ....: ipsum += y.inner_product(z).abs()
991 sage: J1(c) == J1.one()
993 sage: J0(J.one() - c) == J0.one()
997 if not c
.is_idempotent():
998 raise ValueError("element is not idempotent: %s" % c
)
1000 from mjo
.eja
.eja_subalgebra
import FiniteDimensionalEJASubalgebra
1002 # Default these to what they should be if they turn out to be
1003 # trivial, because eigenspaces_left() won't return eigenvalues
1004 # corresponding to trivial spaces (e.g. it returns only the
1005 # eigenspace corresponding to lambda=1 if you take the
1006 # decomposition relative to the identity element).
1007 trivial
= FiniteDimensionalEJASubalgebra(self
, ())
1008 J0
= trivial
# eigenvalue zero
1009 J5
= VectorSpace(self
.base_ring(), 0) # eigenvalue one-half
1010 J1
= trivial
# eigenvalue one
1012 for (eigval
, eigspace
) in c
.operator().matrix().right_eigenspaces():
1013 if eigval
== ~
(self
.base_ring()(2)):
1016 gens
= tuple( self
.from_vector(b
) for b
in eigspace
.basis() )
1017 subalg
= FiniteDimensionalEJASubalgebra(self
,
1025 raise ValueError("unexpected eigenvalue: %s" % eigval
)
1030 def random_element(self
, thorough
=False):
1032 Return a random element of this algebra.
1034 Our algebra superclass method only returns a linear
1035 combination of at most two basis elements. We instead
1036 want the vector space "random element" method that
1037 returns a more diverse selection.
1041 - ``thorough`` -- (boolean; default False) whether or not we
1042 should generate irrational coefficients for the random
1043 element when our base ring is irrational; this slows the
1044 algebra operations to a crawl, but any truly random method
1048 # For a general base ring... maybe we can trust this to do the
1049 # right thing? Unlikely, but.
1050 V
= self
.vector_space()
1051 v
= V
.random_element()
1053 if self
.base_ring() is AA
:
1054 # The "random element" method of the algebraic reals is
1055 # stupid at the moment, and only returns integers between
1056 # -2 and 2, inclusive:
1058 # https://trac.sagemath.org/ticket/30875
1060 # Instead, we implement our own "random vector" method,
1061 # and then coerce that into the algebra. We use the vector
1062 # space degree here instead of the dimension because a
1063 # subalgebra could (for example) be spanned by only two
1064 # vectors, each with five coordinates. We need to
1065 # generate all five coordinates.
1067 v
*= QQbar
.random_element().real()
1069 v
*= QQ
.random_element()
1071 return self
.from_vector(V
.coordinate_vector(v
))
1073 def random_elements(self
, count
, thorough
=False):
1075 Return ``count`` random elements as a tuple.
1079 - ``thorough`` -- (boolean; default False) whether or not we
1080 should generate irrational coefficients for the random
1081 elements when our base ring is irrational; this slows the
1082 algebra operations to a crawl, but any truly random method
1087 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1091 sage: J = JordanSpinEJA(3)
1092 sage: x,y,z = J.random_elements(3)
1093 sage: all( [ x in J, y in J, z in J ])
1095 sage: len( J.random_elements(10) ) == 10
1099 return tuple( self
.random_element(thorough
)
1100 for idx
in range(count
) )
1104 def _charpoly_coefficients(self
):
1106 The `r` polynomial coefficients of the "characteristic polynomial
1111 sage: from mjo.eja.eja_algebra import random_eja
1115 The theory shows that these are all homogeneous polynomials of
1118 sage: set_random_seed()
1119 sage: J = random_eja()
1120 sage: all(p.is_homogeneous() for p in J._charpoly_coefficients())
1124 n
= self
.dimension()
1125 R
= self
.coordinate_polynomial_ring()
1127 F
= R
.fraction_field()
1130 # From a result in my book, these are the entries of the
1131 # basis representation of L_x.
1132 return sum( vars[k
]*self
.monomial(k
).operator().matrix()[i
,j
]
1135 L_x
= matrix(F
, n
, n
, L_x_i_j
)
1138 if self
.rank
.is_in_cache():
1140 # There's no need to pad the system with redundant
1141 # columns if we *know* they'll be redundant.
1144 # Compute an extra power in case the rank is equal to
1145 # the dimension (otherwise, we would stop at x^(r-1)).
1146 x_powers
= [ (L_x
**k
)*self
.one().to_vector()
1147 for k
in range(n
+1) ]
1148 A
= matrix
.column(F
, x_powers
[:n
])
1149 AE
= A
.extended_echelon_form()
1156 # The theory says that only the first "r" coefficients are
1157 # nonzero, and they actually live in the original polynomial
1158 # ring and not the fraction field. We negate them because in
1159 # the actual characteristic polynomial, they get moved to the
1160 # other side where x^r lives. We don't bother to trim A_rref
1161 # down to a square matrix and solve the resulting system,
1162 # because the upper-left r-by-r portion of A_rref is
1163 # guaranteed to be the identity matrix, so e.g.
1165 # A_rref.solve_right(Y)
1167 # would just be returning Y.
1168 return (-E
*b
)[:r
].change_ring(R
)
1173 Return the rank of this EJA.
1175 This is a cached method because we know the rank a priori for
1176 all of the algebras we can construct. Thus we can avoid the
1177 expensive ``_charpoly_coefficients()`` call unless we truly
1178 need to compute the whole characteristic polynomial.
1182 sage: from mjo.eja.eja_algebra import (HadamardEJA,
1183 ....: JordanSpinEJA,
1184 ....: RealSymmetricEJA,
1185 ....: ComplexHermitianEJA,
1186 ....: QuaternionHermitianEJA,
1191 The rank of the Jordan spin algebra is always two::
1193 sage: JordanSpinEJA(2).rank()
1195 sage: JordanSpinEJA(3).rank()
1197 sage: JordanSpinEJA(4).rank()
1200 The rank of the `n`-by-`n` Hermitian real, complex, or
1201 quaternion matrices is `n`::
1203 sage: RealSymmetricEJA(4).rank()
1205 sage: ComplexHermitianEJA(3).rank()
1207 sage: QuaternionHermitianEJA(2).rank()
1212 Ensure that every EJA that we know how to construct has a
1213 positive integer rank, unless the algebra is trivial in
1214 which case its rank will be zero::
1216 sage: set_random_seed()
1217 sage: J = random_eja()
1221 sage: r > 0 or (r == 0 and J.is_trivial())
1224 Ensure that computing the rank actually works, since the ranks
1225 of all simple algebras are known and will be cached by default::
1227 sage: set_random_seed() # long time
1228 sage: J = random_eja() # long time
1229 sage: cached = J.rank() # long time
1230 sage: J.rank.clear_cache() # long time
1231 sage: J.rank() == cached # long time
1235 return len(self
._charpoly
_coefficients
())
1238 def vector_space(self
):
1240 Return the vector space that underlies this algebra.
1244 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1248 sage: J = RealSymmetricEJA(2)
1249 sage: J.vector_space()
1250 Vector space of dimension 3 over...
1253 return self
.zero().to_vector().parent().ambient_vector_space()
1256 Element
= FiniteDimensionalEJAElement
1258 class RationalBasisEJA(FiniteDimensionalEJA
):
1260 New class for algebras whose supplied basis elements have all rational entries.
1264 sage: from mjo.eja.eja_algebra import BilinearFormEJA
1268 The supplied basis is orthonormalized by default::
1270 sage: B = matrix(QQ, [[1, 0, 0], [0, 25, -32], [0, -32, 41]])
1271 sage: J = BilinearFormEJA(B)
1272 sage: J.matrix_basis()
1289 # Abuse the check_field parameter to check that the entries of
1290 # out basis (in ambient coordinates) are in the field QQ.
1291 if not all( all(b_i
in QQ
for b_i
in b
.list()) for b
in basis
):
1292 raise TypeError("basis not rational")
1294 self
._rational
_algebra
= None
1296 # There's no point in constructing the extra algebra if this
1297 # one is already rational.
1299 # Note: the same Jordan and inner-products work here,
1300 # because they are necessarily defined with respect to
1301 # ambient coordinates and not any particular basis.
1302 self
._rational
_algebra
= FiniteDimensionalEJA(
1307 orthonormalize
=False,
1311 super().__init
__(basis
,
1315 check_field
=check_field
,
1319 def _charpoly_coefficients(self
):
1323 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
1324 ....: JordanSpinEJA)
1328 The base ring of the resulting polynomial coefficients is what
1329 it should be, and not the rationals (unless the algebra was
1330 already over the rationals)::
1332 sage: J = JordanSpinEJA(3)
1333 sage: J._charpoly_coefficients()
1334 (X1^2 - X2^2 - X3^2, -2*X1)
1335 sage: a0 = J._charpoly_coefficients()[0]
1337 Algebraic Real Field
1338 sage: a0.base_ring()
1339 Algebraic Real Field
1342 if self
._rational
_algebra
is None:
1343 # There's no need to construct *another* algebra over the
1344 # rationals if this one is already over the
1345 # rationals. Likewise, if we never orthonormalized our
1346 # basis, we might as well just use the given one.
1347 return super()._charpoly
_coefficients
()
1349 # Do the computation over the rationals. The answer will be
1350 # the same, because all we've done is a change of basis.
1351 # Then, change back from QQ to our real base ring
1352 a
= ( a_i
.change_ring(self
.base_ring())
1353 for a_i
in self
._rational
_algebra
._charpoly
_coefficients
() )
1355 if self
._deortho
_matrix
is None:
1356 # This can happen if our base ring was, say, AA and we
1357 # chose not to (or didn't need to) orthonormalize. It's
1358 # still faster to do the computations over QQ even if
1359 # the numbers in the boxes stay the same.
1362 # Otherwise, convert the coordinate variables back to the
1363 # deorthonormalized ones.
1364 R
= self
.coordinate_polynomial_ring()
1365 from sage
.modules
.free_module_element
import vector
1366 X
= vector(R
, R
.gens())
1367 BX
= self
._deortho
_matrix
*X
1369 subs_dict
= { X[i]: BX[i] for i in range(len(X)) }
1370 return tuple( a_i
.subs(subs_dict
) for a_i
in a
)
1372 class ConcreteEJA(RationalBasisEJA
):
1374 A class for the Euclidean Jordan algebras that we know by name.
1376 These are the Jordan algebras whose basis, multiplication table,
1377 rank, and so on are known a priori. More to the point, they are
1378 the Euclidean Jordan algebras for which we are able to conjure up
1379 a "random instance."
1383 sage: from mjo.eja.eja_algebra import ConcreteEJA
1387 Our basis is normalized with respect to the algebra's inner
1388 product, unless we specify otherwise::
1390 sage: set_random_seed()
1391 sage: J = ConcreteEJA.random_instance()
1392 sage: all( b.norm() == 1 for b in J.gens() )
1395 Since our basis is orthonormal with respect to the algebra's inner
1396 product, and since we know that this algebra is an EJA, any
1397 left-multiplication operator's matrix will be symmetric because
1398 natural->EJA basis representation is an isometry and within the
1399 EJA the operator is self-adjoint by the Jordan axiom::
1401 sage: set_random_seed()
1402 sage: J = ConcreteEJA.random_instance()
1403 sage: x = J.random_element()
1404 sage: x.operator().is_self_adjoint()
1409 def _max_random_instance_size():
1411 Return an integer "size" that is an upper bound on the size of
1412 this algebra when it is used in a random test
1413 case. Unfortunately, the term "size" is ambiguous -- when
1414 dealing with `R^n` under either the Hadamard or Jordan spin
1415 product, the "size" refers to the dimension `n`. When dealing
1416 with a matrix algebra (real symmetric or complex/quaternion
1417 Hermitian), it refers to the size of the matrix, which is far
1418 less than the dimension of the underlying vector space.
1420 This method must be implemented in each subclass.
1422 raise NotImplementedError
1425 def random_instance(cls
, *args
, **kwargs
):
1427 Return a random instance of this type of algebra.
1429 This method should be implemented in each subclass.
1431 from sage
.misc
.prandom
import choice
1432 eja_class
= choice(cls
.__subclasses
__())
1434 # These all bubble up to the RationalBasisEJA superclass
1435 # constructor, so any (kw)args valid there are also valid
1437 return eja_class
.random_instance(*args
, **kwargs
)
1442 def dimension_over_reals():
1444 The dimension of this matrix's base ring over the reals.
1446 The reals are dimension one over themselves, obviously; that's
1447 just `\mathbb{R}^{1}`. Likewise, the complex numbers `a + bi`
1448 have dimension two. Finally, the quaternions have dimension
1449 four over the reals.
1451 This is used to determine the size of the matrix returned from
1452 :meth:`real_embed`, among other things.
1454 raise NotImplementedError
1457 def real_embed(cls
,M
):
1459 Embed the matrix ``M`` into a space of real matrices.
1461 The matrix ``M`` can have entries in any field at the moment:
1462 the real numbers, complex numbers, or quaternions. And although
1463 they are not a field, we can probably support octonions at some
1464 point, too. This function returns a real matrix that "acts like"
1465 the original with respect to matrix multiplication; i.e.
1467 real_embed(M*N) = real_embed(M)*real_embed(N)
1470 if M
.ncols() != M
.nrows():
1471 raise ValueError("the matrix 'M' must be square")
1476 def real_unembed(cls
,M
):
1478 The inverse of :meth:`real_embed`.
1480 if M
.ncols() != M
.nrows():
1481 raise ValueError("the matrix 'M' must be square")
1482 if not ZZ(M
.nrows()).mod(cls
.dimension_over_reals()).is_zero():
1483 raise ValueError("the matrix 'M' must be a real embedding")
1487 def jordan_product(X
,Y
):
1488 return (X
*Y
+ Y
*X
)/2
1491 def trace_inner_product(cls
,X
,Y
):
1493 Compute the trace inner-product of two real-embeddings.
1497 sage: from mjo.eja.eja_algebra import (RealSymmetricEJA,
1498 ....: ComplexHermitianEJA,
1499 ....: QuaternionHermitianEJA)
1503 This gives the same answer as it would if we computed the trace
1504 from the unembedded (original) matrices::
1506 sage: set_random_seed()
1507 sage: J = RealSymmetricEJA.random_instance()
1508 sage: x,y = J.random_elements(2)
1509 sage: Xe = x.to_matrix()
1510 sage: Ye = y.to_matrix()
1511 sage: X = J.real_unembed(Xe)
1512 sage: Y = J.real_unembed(Ye)
1513 sage: expected = (X*Y).trace()
1514 sage: actual = J.trace_inner_product(Xe,Ye)
1515 sage: actual == expected
1520 sage: set_random_seed()
1521 sage: J = ComplexHermitianEJA.random_instance()
1522 sage: x,y = J.random_elements(2)
1523 sage: Xe = x.to_matrix()
1524 sage: Ye = y.to_matrix()
1525 sage: X = J.real_unembed(Xe)
1526 sage: Y = J.real_unembed(Ye)
1527 sage: expected = (X*Y).trace().real()
1528 sage: actual = J.trace_inner_product(Xe,Ye)
1529 sage: actual == expected
1534 sage: set_random_seed()
1535 sage: J = QuaternionHermitianEJA.random_instance()
1536 sage: x,y = J.random_elements(2)
1537 sage: Xe = x.to_matrix()
1538 sage: Ye = y.to_matrix()
1539 sage: X = J.real_unembed(Xe)
1540 sage: Y = J.real_unembed(Ye)
1541 sage: expected = (X*Y).trace().coefficient_tuple()[0]
1542 sage: actual = J.trace_inner_product(Xe,Ye)
1543 sage: actual == expected
1547 Xu
= cls
.real_unembed(X
)
1548 Yu
= cls
.real_unembed(Y
)
1549 tr
= (Xu
*Yu
).trace()
1552 # Works in QQ, AA, RDF, et cetera.
1554 except AttributeError:
1555 # A quaternion doesn't have a real() method, but does
1556 # have coefficient_tuple() method that returns the
1557 # coefficients of 1, i, j, and k -- in that order.
1558 return tr
.coefficient_tuple()[0]
1561 class RealMatrixEJA(MatrixEJA
):
1563 def dimension_over_reals():
1567 class RealSymmetricEJA(ConcreteEJA
, RealMatrixEJA
):
1569 The rank-n simple EJA consisting of real symmetric n-by-n
1570 matrices, the usual symmetric Jordan product, and the trace inner
1571 product. It has dimension `(n^2 + n)/2` over the reals.
1575 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1579 sage: J = RealSymmetricEJA(2)
1580 sage: e0, e1, e2 = J.gens()
1588 In theory, our "field" can be any subfield of the reals::
1590 sage: RealSymmetricEJA(2, field=RDF)
1591 Euclidean Jordan algebra of dimension 3 over Real Double Field
1592 sage: RealSymmetricEJA(2, field=RR)
1593 Euclidean Jordan algebra of dimension 3 over Real Field with
1594 53 bits of precision
1598 The dimension of this algebra is `(n^2 + n) / 2`::
1600 sage: set_random_seed()
1601 sage: n_max = RealSymmetricEJA._max_random_instance_size()
1602 sage: n = ZZ.random_element(1, n_max)
1603 sage: J = RealSymmetricEJA(n)
1604 sage: J.dimension() == (n^2 + n)/2
1607 The Jordan multiplication is what we think it is::
1609 sage: set_random_seed()
1610 sage: J = RealSymmetricEJA.random_instance()
1611 sage: x,y = J.random_elements(2)
1612 sage: actual = (x*y).to_matrix()
1613 sage: X = x.to_matrix()
1614 sage: Y = y.to_matrix()
1615 sage: expected = (X*Y + Y*X)/2
1616 sage: actual == expected
1618 sage: J(expected) == x*y
1621 We can change the generator prefix::
1623 sage: RealSymmetricEJA(3, prefix='q').gens()
1624 (q0, q1, q2, q3, q4, q5)
1626 We can construct the (trivial) algebra of rank zero::
1628 sage: RealSymmetricEJA(0)
1629 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1633 def _denormalized_basis(cls
, n
):
1635 Return a basis for the space of real symmetric n-by-n matrices.
1639 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1643 sage: set_random_seed()
1644 sage: n = ZZ.random_element(1,5)
1645 sage: B = RealSymmetricEJA._denormalized_basis(n)
1646 sage: all( M.is_symmetric() for M in B)
1650 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1654 for j
in range(i
+1):
1655 Eij
= matrix(ZZ
, n
, lambda k
,l
: k
==i
and l
==j
)
1659 Sij
= Eij
+ Eij
.transpose()
1665 def _max_random_instance_size():
1666 return 4 # Dimension 10
1669 def random_instance(cls
, **kwargs
):
1671 Return a random instance of this type of algebra.
1673 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1674 return cls(n
, **kwargs
)
1676 def __init__(self
, n
, **kwargs
):
1677 # We know this is a valid EJA, but will double-check
1678 # if the user passes check_axioms=True.
1679 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
1681 super(RealSymmetricEJA
, self
).__init
__(self
._denormalized
_basis
(n
),
1682 self
.jordan_product
,
1683 self
.trace_inner_product
,
1686 # TODO: this could be factored out somehow, but is left here
1687 # because the MatrixEJA is not presently a subclass of the
1688 # FDEJA class that defines rank() and one().
1689 self
.rank
.set_cache(n
)
1690 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
1691 self
.one
.set_cache(self(idV
))
1695 class ComplexMatrixEJA(MatrixEJA
):
1696 # A manual dictionary-cache for the complex_extension() method,
1697 # since apparently @classmethods can't also be @cached_methods.
1698 _complex_extension
= {}
1701 def complex_extension(cls
,field
):
1703 The complex field that we embed/unembed, as an extension
1704 of the given ``field``.
1706 if field
in cls
._complex
_extension
:
1707 return cls
._complex
_extension
[field
]
1709 # Sage doesn't know how to adjoin the complex "i" (the root of
1710 # x^2 + 1) to a field in a general way. Here, we just enumerate
1711 # all of the cases that I have cared to support so far.
1713 # Sage doesn't know how to embed AA into QQbar, i.e. how
1714 # to adjoin sqrt(-1) to AA.
1716 elif not field
.is_exact():
1718 F
= field
.complex_field()
1720 # Works for QQ and... maybe some other fields.
1721 R
= PolynomialRing(field
, 'z')
1723 F
= field
.extension(z
**2 + 1, 'I', embedding
=CLF(-1).sqrt())
1725 cls
._complex
_extension
[field
] = F
1729 def dimension_over_reals():
1733 def real_embed(cls
,M
):
1735 Embed the n-by-n complex matrix ``M`` into the space of real
1736 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1737 bi` to the block matrix ``[[a,b],[-b,a]]``.
1741 sage: from mjo.eja.eja_algebra import ComplexMatrixEJA
1745 sage: F = QuadraticField(-1, 'I')
1746 sage: x1 = F(4 - 2*i)
1747 sage: x2 = F(1 + 2*i)
1750 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1751 sage: ComplexMatrixEJA.real_embed(M)
1760 Embedding is a homomorphism (isomorphism, in fact)::
1762 sage: set_random_seed()
1763 sage: n = ZZ.random_element(3)
1764 sage: F = QuadraticField(-1, 'I')
1765 sage: X = random_matrix(F, n)
1766 sage: Y = random_matrix(F, n)
1767 sage: Xe = ComplexMatrixEJA.real_embed(X)
1768 sage: Ye = ComplexMatrixEJA.real_embed(Y)
1769 sage: XYe = ComplexMatrixEJA.real_embed(X*Y)
1774 super(ComplexMatrixEJA
,cls
).real_embed(M
)
1777 # We don't need any adjoined elements...
1778 field
= M
.base_ring().base_ring()
1784 blocks
.append(matrix(field
, 2, [ [ a
, b
],
1787 return matrix
.block(field
, n
, blocks
)
1791 def real_unembed(cls
,M
):
1793 The inverse of _embed_complex_matrix().
1797 sage: from mjo.eja.eja_algebra import ComplexMatrixEJA
1801 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1802 ....: [-2, 1, -4, 3],
1803 ....: [ 9, 10, 11, 12],
1804 ....: [-10, 9, -12, 11] ])
1805 sage: ComplexMatrixEJA.real_unembed(A)
1807 [ 10*I + 9 12*I + 11]
1811 Unembedding is the inverse of embedding::
1813 sage: set_random_seed()
1814 sage: F = QuadraticField(-1, 'I')
1815 sage: M = random_matrix(F, 3)
1816 sage: Me = ComplexMatrixEJA.real_embed(M)
1817 sage: ComplexMatrixEJA.real_unembed(Me) == M
1821 super(ComplexMatrixEJA
,cls
).real_unembed(M
)
1823 d
= cls
.dimension_over_reals()
1824 F
= cls
.complex_extension(M
.base_ring())
1827 # Go top-left to bottom-right (reading order), converting every
1828 # 2-by-2 block we see to a single complex element.
1830 for k
in range(n
/d
):
1831 for j
in range(n
/d
):
1832 submat
= M
[d
*k
:d
*k
+d
,d
*j
:d
*j
+d
]
1833 if submat
[0,0] != submat
[1,1]:
1834 raise ValueError('bad on-diagonal submatrix')
1835 if submat
[0,1] != -submat
[1,0]:
1836 raise ValueError('bad off-diagonal submatrix')
1837 z
= submat
[0,0] + submat
[0,1]*i
1840 return matrix(F
, n
/d
, elements
)
1843 class ComplexHermitianEJA(ConcreteEJA
, ComplexMatrixEJA
):
1845 The rank-n simple EJA consisting of complex Hermitian n-by-n
1846 matrices over the real numbers, the usual symmetric Jordan product,
1847 and the real-part-of-trace inner product. It has dimension `n^2` over
1852 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1856 In theory, our "field" can be any subfield of the reals::
1858 sage: ComplexHermitianEJA(2, field=RDF)
1859 Euclidean Jordan algebra of dimension 4 over Real Double Field
1860 sage: ComplexHermitianEJA(2, field=RR)
1861 Euclidean Jordan algebra of dimension 4 over Real Field with
1862 53 bits of precision
1866 The dimension of this algebra is `n^2`::
1868 sage: set_random_seed()
1869 sage: n_max = ComplexHermitianEJA._max_random_instance_size()
1870 sage: n = ZZ.random_element(1, n_max)
1871 sage: J = ComplexHermitianEJA(n)
1872 sage: J.dimension() == n^2
1875 The Jordan multiplication is what we think it is::
1877 sage: set_random_seed()
1878 sage: J = ComplexHermitianEJA.random_instance()
1879 sage: x,y = J.random_elements(2)
1880 sage: actual = (x*y).to_matrix()
1881 sage: X = x.to_matrix()
1882 sage: Y = y.to_matrix()
1883 sage: expected = (X*Y + Y*X)/2
1884 sage: actual == expected
1886 sage: J(expected) == x*y
1889 We can change the generator prefix::
1891 sage: ComplexHermitianEJA(2, prefix='z').gens()
1894 We can construct the (trivial) algebra of rank zero::
1896 sage: ComplexHermitianEJA(0)
1897 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1902 def _denormalized_basis(cls
, n
):
1904 Returns a basis for the space of complex Hermitian n-by-n matrices.
1906 Why do we embed these? Basically, because all of numerical linear
1907 algebra assumes that you're working with vectors consisting of `n`
1908 entries from a field and scalars from the same field. There's no way
1909 to tell SageMath that (for example) the vectors contain complex
1910 numbers, while the scalar field is real.
1914 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1918 sage: set_random_seed()
1919 sage: n = ZZ.random_element(1,5)
1920 sage: B = ComplexHermitianEJA._denormalized_basis(n)
1921 sage: all( M.is_symmetric() for M in B)
1926 R
= PolynomialRing(field
, 'z')
1928 F
= field
.extension(z
**2 + 1, 'I')
1931 # This is like the symmetric case, but we need to be careful:
1933 # * We want conjugate-symmetry, not just symmetry.
1934 # * The diagonal will (as a result) be real.
1937 Eij
= matrix
.zero(F
,n
)
1939 for j
in range(i
+1):
1943 Sij
= cls
.real_embed(Eij
)
1946 # The second one has a minus because it's conjugated.
1947 Eij
[j
,i
] = 1 # Eij = Eij + Eij.transpose()
1948 Sij_real
= cls
.real_embed(Eij
)
1950 # Eij = I*Eij - I*Eij.transpose()
1953 Sij_imag
= cls
.real_embed(Eij
)
1959 # Since we embedded these, we can drop back to the "field" that we
1960 # started with instead of the complex extension "F".
1961 return tuple( s
.change_ring(field
) for s
in S
)
1964 def __init__(self
, n
, **kwargs
):
1965 # We know this is a valid EJA, but will double-check
1966 # if the user passes check_axioms=True.
1967 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
1969 super(ComplexHermitianEJA
, self
).__init
__(self
._denormalized
_basis
(n
),
1970 self
.jordan_product
,
1971 self
.trace_inner_product
,
1973 # TODO: this could be factored out somehow, but is left here
1974 # because the MatrixEJA is not presently a subclass of the
1975 # FDEJA class that defines rank() and one().
1976 self
.rank
.set_cache(n
)
1977 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
1978 self
.one
.set_cache(self(idV
))
1981 def _max_random_instance_size():
1982 return 3 # Dimension 9
1985 def random_instance(cls
, **kwargs
):
1987 Return a random instance of this type of algebra.
1989 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1990 return cls(n
, **kwargs
)
1992 class QuaternionMatrixEJA(MatrixEJA
):
1994 # A manual dictionary-cache for the quaternion_extension() method,
1995 # since apparently @classmethods can't also be @cached_methods.
1996 _quaternion_extension
= {}
1999 def quaternion_extension(cls
,field
):
2001 The quaternion field that we embed/unembed, as an extension
2002 of the given ``field``.
2004 if field
in cls
._quaternion
_extension
:
2005 return cls
._quaternion
_extension
[field
]
2007 Q
= QuaternionAlgebra(field
,-1,-1)
2009 cls
._quaternion
_extension
[field
] = Q
2013 def dimension_over_reals():
2017 def real_embed(cls
,M
):
2019 Embed the n-by-n quaternion matrix ``M`` into the space of real
2020 matrices of size 4n-by-4n by first sending each quaternion entry `z
2021 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
2022 c+di],[-c + di, a-bi]]`, and then embedding those into a real
2027 sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA
2031 sage: Q = QuaternionAlgebra(QQ,-1,-1)
2032 sage: i,j,k = Q.gens()
2033 sage: x = 1 + 2*i + 3*j + 4*k
2034 sage: M = matrix(Q, 1, [[x]])
2035 sage: QuaternionMatrixEJA.real_embed(M)
2041 Embedding is a homomorphism (isomorphism, in fact)::
2043 sage: set_random_seed()
2044 sage: n = ZZ.random_element(2)
2045 sage: Q = QuaternionAlgebra(QQ,-1,-1)
2046 sage: X = random_matrix(Q, n)
2047 sage: Y = random_matrix(Q, n)
2048 sage: Xe = QuaternionMatrixEJA.real_embed(X)
2049 sage: Ye = QuaternionMatrixEJA.real_embed(Y)
2050 sage: XYe = QuaternionMatrixEJA.real_embed(X*Y)
2055 super(QuaternionMatrixEJA
,cls
).real_embed(M
)
2056 quaternions
= M
.base_ring()
2059 F
= QuadraticField(-1, 'I')
2064 t
= z
.coefficient_tuple()
2069 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
2070 [-c
+ d
*i
, a
- b
*i
]])
2071 realM
= ComplexMatrixEJA
.real_embed(cplxM
)
2072 blocks
.append(realM
)
2074 # We should have real entries by now, so use the realest field
2075 # we've got for the return value.
2076 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
2081 def real_unembed(cls
,M
):
2083 The inverse of _embed_quaternion_matrix().
2087 sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA
2091 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
2092 ....: [-2, 1, -4, 3],
2093 ....: [-3, 4, 1, -2],
2094 ....: [-4, -3, 2, 1]])
2095 sage: QuaternionMatrixEJA.real_unembed(M)
2096 [1 + 2*i + 3*j + 4*k]
2100 Unembedding is the inverse of embedding::
2102 sage: set_random_seed()
2103 sage: Q = QuaternionAlgebra(QQ, -1, -1)
2104 sage: M = random_matrix(Q, 3)
2105 sage: Me = QuaternionMatrixEJA.real_embed(M)
2106 sage: QuaternionMatrixEJA.real_unembed(Me) == M
2110 super(QuaternionMatrixEJA
,cls
).real_unembed(M
)
2112 d
= cls
.dimension_over_reals()
2114 # Use the base ring of the matrix to ensure that its entries can be
2115 # multiplied by elements of the quaternion algebra.
2116 Q
= cls
.quaternion_extension(M
.base_ring())
2119 # Go top-left to bottom-right (reading order), converting every
2120 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
2123 for l
in range(n
/d
):
2124 for m
in range(n
/d
):
2125 submat
= ComplexMatrixEJA
.real_unembed(
2126 M
[d
*l
:d
*l
+d
,d
*m
:d
*m
+d
] )
2127 if submat
[0,0] != submat
[1,1].conjugate():
2128 raise ValueError('bad on-diagonal submatrix')
2129 if submat
[0,1] != -submat
[1,0].conjugate():
2130 raise ValueError('bad off-diagonal submatrix')
2131 z
= submat
[0,0].real()
2132 z
+= submat
[0,0].imag()*i
2133 z
+= submat
[0,1].real()*j
2134 z
+= submat
[0,1].imag()*k
2137 return matrix(Q
, n
/d
, elements
)
2140 class QuaternionHermitianEJA(ConcreteEJA
, QuaternionMatrixEJA
):
2142 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
2143 matrices, the usual symmetric Jordan product, and the
2144 real-part-of-trace inner product. It has dimension `2n^2 - n` over
2149 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2153 In theory, our "field" can be any subfield of the reals::
2155 sage: QuaternionHermitianEJA(2, field=RDF)
2156 Euclidean Jordan algebra of dimension 6 over Real Double Field
2157 sage: QuaternionHermitianEJA(2, field=RR)
2158 Euclidean Jordan algebra of dimension 6 over Real Field with
2159 53 bits of precision
2163 The dimension of this algebra is `2*n^2 - n`::
2165 sage: set_random_seed()
2166 sage: n_max = QuaternionHermitianEJA._max_random_instance_size()
2167 sage: n = ZZ.random_element(1, n_max)
2168 sage: J = QuaternionHermitianEJA(n)
2169 sage: J.dimension() == 2*(n^2) - n
2172 The Jordan multiplication is what we think it is::
2174 sage: set_random_seed()
2175 sage: J = QuaternionHermitianEJA.random_instance()
2176 sage: x,y = J.random_elements(2)
2177 sage: actual = (x*y).to_matrix()
2178 sage: X = x.to_matrix()
2179 sage: Y = y.to_matrix()
2180 sage: expected = (X*Y + Y*X)/2
2181 sage: actual == expected
2183 sage: J(expected) == x*y
2186 We can change the generator prefix::
2188 sage: QuaternionHermitianEJA(2, prefix='a').gens()
2189 (a0, a1, a2, a3, a4, a5)
2191 We can construct the (trivial) algebra of rank zero::
2193 sage: QuaternionHermitianEJA(0)
2194 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2198 def _denormalized_basis(cls
, n
):
2200 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
2202 Why do we embed these? Basically, because all of numerical
2203 linear algebra assumes that you're working with vectors consisting
2204 of `n` entries from a field and scalars from the same field. There's
2205 no way to tell SageMath that (for example) the vectors contain
2206 complex numbers, while the scalar field is real.
2210 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2214 sage: set_random_seed()
2215 sage: n = ZZ.random_element(1,5)
2216 sage: B = QuaternionHermitianEJA._denormalized_basis(n)
2217 sage: all( M.is_symmetric() for M in B )
2222 Q
= QuaternionAlgebra(QQ
,-1,-1)
2225 # This is like the symmetric case, but we need to be careful:
2227 # * We want conjugate-symmetry, not just symmetry.
2228 # * The diagonal will (as a result) be real.
2231 Eij
= matrix
.zero(Q
,n
)
2233 for j
in range(i
+1):
2237 Sij
= cls
.real_embed(Eij
)
2240 # The second, third, and fourth ones have a minus
2241 # because they're conjugated.
2242 # Eij = Eij + Eij.transpose()
2244 Sij_real
= cls
.real_embed(Eij
)
2246 # Eij = I*(Eij - Eij.transpose())
2249 Sij_I
= cls
.real_embed(Eij
)
2251 # Eij = J*(Eij - Eij.transpose())
2254 Sij_J
= cls
.real_embed(Eij
)
2256 # Eij = K*(Eij - Eij.transpose())
2259 Sij_K
= cls
.real_embed(Eij
)
2265 # Since we embedded these, we can drop back to the "field" that we
2266 # started with instead of the quaternion algebra "Q".
2267 return tuple( s
.change_ring(field
) for s
in S
)
2270 def __init__(self
, n
, **kwargs
):
2271 # We know this is a valid EJA, but will double-check
2272 # if the user passes check_axioms=True.
2273 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2275 super(QuaternionHermitianEJA
, self
).__init
__(self
._denormalized
_basis
(n
),
2276 self
.jordan_product
,
2277 self
.trace_inner_product
,
2279 # TODO: this could be factored out somehow, but is left here
2280 # because the MatrixEJA is not presently a subclass of the
2281 # FDEJA class that defines rank() and one().
2282 self
.rank
.set_cache(n
)
2283 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
2284 self
.one
.set_cache(self(idV
))
2288 def _max_random_instance_size():
2290 The maximum rank of a random QuaternionHermitianEJA.
2292 return 2 # Dimension 6
2295 def random_instance(cls
, **kwargs
):
2297 Return a random instance of this type of algebra.
2299 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2300 return cls(n
, **kwargs
)
2303 class HadamardEJA(ConcreteEJA
):
2305 Return the Euclidean Jordan Algebra corresponding to the set
2306 `R^n` under the Hadamard product.
2308 Note: this is nothing more than the Cartesian product of ``n``
2309 copies of the spin algebra. Once Cartesian product algebras
2310 are implemented, this can go.
2314 sage: from mjo.eja.eja_algebra import HadamardEJA
2318 This multiplication table can be verified by hand::
2320 sage: J = HadamardEJA(3)
2321 sage: e0,e1,e2 = J.gens()
2337 We can change the generator prefix::
2339 sage: HadamardEJA(3, prefix='r').gens()
2343 def __init__(self
, n
, **kwargs
):
2345 jordan_product
= lambda x
,y
: x
2346 inner_product
= lambda x
,y
: x
2348 def jordan_product(x
,y
):
2350 return P( xi
*yi
for (xi
,yi
) in zip(x
,y
) )
2352 def inner_product(x
,y
):
2355 # New defaults for keyword arguments. Don't orthonormalize
2356 # because our basis is already orthonormal with respect to our
2357 # inner-product. Don't check the axioms, because we know this
2358 # is a valid EJA... but do double-check if the user passes
2359 # check_axioms=True. Note: we DON'T override the "check_field"
2360 # default here, because the user can pass in a field!
2361 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2362 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2364 column_basis
= tuple( b
.column() for b
in FreeModule(ZZ
, n
).basis() )
2365 super().__init
__(column_basis
, jordan_product
, inner_product
, **kwargs
)
2366 self
.rank
.set_cache(n
)
2369 self
.one
.set_cache( self
.zero() )
2371 self
.one
.set_cache( sum(self
.gens()) )
2374 def _max_random_instance_size():
2376 The maximum dimension of a random HadamardEJA.
2381 def random_instance(cls
, **kwargs
):
2383 Return a random instance of this type of algebra.
2385 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2386 return cls(n
, **kwargs
)
2389 class BilinearFormEJA(ConcreteEJA
):
2391 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2392 with the half-trace inner product and jordan product ``x*y =
2393 (<Bx,y>,y_bar>, x0*y_bar + y0*x_bar)`` where `B = 1 \times B22` is
2394 a symmetric positive-definite "bilinear form" matrix. Its
2395 dimension is the size of `B`, and it has rank two in dimensions
2396 larger than two. It reduces to the ``JordanSpinEJA`` when `B` is
2397 the identity matrix of order ``n``.
2399 We insist that the one-by-one upper-left identity block of `B` be
2400 passed in as well so that we can be passed a matrix of size zero
2401 to construct a trivial algebra.
2405 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
2406 ....: JordanSpinEJA)
2410 When no bilinear form is specified, the identity matrix is used,
2411 and the resulting algebra is the Jordan spin algebra::
2413 sage: B = matrix.identity(AA,3)
2414 sage: J0 = BilinearFormEJA(B)
2415 sage: J1 = JordanSpinEJA(3)
2416 sage: J0.multiplication_table() == J0.multiplication_table()
2419 An error is raised if the matrix `B` does not correspond to a
2420 positive-definite bilinear form::
2422 sage: B = matrix.random(QQ,2,3)
2423 sage: J = BilinearFormEJA(B)
2424 Traceback (most recent call last):
2426 ValueError: bilinear form is not positive-definite
2427 sage: B = matrix.zero(QQ,3)
2428 sage: J = BilinearFormEJA(B)
2429 Traceback (most recent call last):
2431 ValueError: bilinear form is not positive-definite
2435 We can create a zero-dimensional algebra::
2437 sage: B = matrix.identity(AA,0)
2438 sage: J = BilinearFormEJA(B)
2442 We can check the multiplication condition given in the Jordan, von
2443 Neumann, and Wigner paper (and also discussed on my "On the
2444 symmetry..." paper). Note that this relies heavily on the standard
2445 choice of basis, as does anything utilizing the bilinear form
2446 matrix. We opt not to orthonormalize the basis, because if we
2447 did, we would have to normalize the `s_{i}` in a similar manner::
2449 sage: set_random_seed()
2450 sage: n = ZZ.random_element(5)
2451 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2452 sage: B11 = matrix.identity(QQ,1)
2453 sage: B22 = M.transpose()*M
2454 sage: B = block_matrix(2,2,[ [B11,0 ],
2456 sage: J = BilinearFormEJA(B, orthonormalize=False)
2457 sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
2458 sage: V = J.vector_space()
2459 sage: sis = [ J( V([0] + (M.inverse()*ei).list()).column() )
2460 ....: for ei in eis ]
2461 sage: actual = [ sis[i]*sis[j]
2462 ....: for i in range(n-1)
2463 ....: for j in range(n-1) ]
2464 sage: expected = [ J.one() if i == j else J.zero()
2465 ....: for i in range(n-1)
2466 ....: for j in range(n-1) ]
2467 sage: actual == expected
2471 def __init__(self
, B
, **kwargs
):
2472 # The matrix "B" is supplied by the user in most cases,
2473 # so it makes sense to check whether or not its positive-
2474 # definite unless we are specifically asked not to...
2475 if ("check_axioms" not in kwargs
) or kwargs
["check_axioms"]:
2476 if not B
.is_positive_definite():
2477 raise ValueError("bilinear form is not positive-definite")
2479 # However, all of the other data for this EJA is computed
2480 # by us in manner that guarantees the axioms are
2481 # satisfied. So, again, unless we are specifically asked to
2482 # verify things, we'll skip the rest of the checks.
2483 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2485 def inner_product(x
,y
):
2486 return (y
.T
*B
*x
)[0,0]
2488 def jordan_product(x
,y
):
2494 z0
= inner_product(y
,x
)
2495 zbar
= y0
*xbar
+ x0
*ybar
2496 return P([z0
] + zbar
.list())
2499 column_basis
= tuple( b
.column() for b
in FreeModule(ZZ
, n
).basis() )
2500 super(BilinearFormEJA
, self
).__init
__(column_basis
,
2505 # The rank of this algebra is two, unless we're in a
2506 # one-dimensional ambient space (because the rank is bounded
2507 # by the ambient dimension).
2508 self
.rank
.set_cache(min(n
,2))
2511 self
.one
.set_cache( self
.zero() )
2513 self
.one
.set_cache( self
.monomial(0) )
2516 def _max_random_instance_size():
2518 The maximum dimension of a random BilinearFormEJA.
2523 def random_instance(cls
, **kwargs
):
2525 Return a random instance of this algebra.
2527 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2529 B
= matrix
.identity(ZZ
, n
)
2530 return cls(B
, **kwargs
)
2532 B11
= matrix
.identity(ZZ
, 1)
2533 M
= matrix
.random(ZZ
, n
-1)
2534 I
= matrix
.identity(ZZ
, n
-1)
2536 while alpha
.is_zero():
2537 alpha
= ZZ
.random_element().abs()
2538 B22
= M
.transpose()*M
+ alpha
*I
2540 from sage
.matrix
.special
import block_matrix
2541 B
= block_matrix(2,2, [ [B11
, ZZ(0) ],
2544 return cls(B
, **kwargs
)
2547 class JordanSpinEJA(BilinearFormEJA
):
2549 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2550 with the usual inner product and jordan product ``x*y =
2551 (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2556 sage: from mjo.eja.eja_algebra import JordanSpinEJA
2560 This multiplication table can be verified by hand::
2562 sage: J = JordanSpinEJA(4)
2563 sage: e0,e1,e2,e3 = J.gens()
2579 We can change the generator prefix::
2581 sage: JordanSpinEJA(2, prefix='B').gens()
2586 Ensure that we have the usual inner product on `R^n`::
2588 sage: set_random_seed()
2589 sage: J = JordanSpinEJA.random_instance()
2590 sage: x,y = J.random_elements(2)
2591 sage: actual = x.inner_product(y)
2592 sage: expected = x.to_vector().inner_product(y.to_vector())
2593 sage: actual == expected
2597 def __init__(self
, n
, **kwargs
):
2598 # This is a special case of the BilinearFormEJA with the
2599 # identity matrix as its bilinear form.
2600 B
= matrix
.identity(ZZ
, n
)
2602 # Don't orthonormalize because our basis is already
2603 # orthonormal with respect to our inner-product.
2604 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2606 # But also don't pass check_field=False here, because the user
2607 # can pass in a field!
2608 super(JordanSpinEJA
, self
).__init
__(B
, **kwargs
)
2611 def _max_random_instance_size():
2613 The maximum dimension of a random JordanSpinEJA.
2618 def random_instance(cls
, **kwargs
):
2620 Return a random instance of this type of algebra.
2622 Needed here to override the implementation for ``BilinearFormEJA``.
2624 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2625 return cls(n
, **kwargs
)
2628 class TrivialEJA(ConcreteEJA
):
2630 The trivial Euclidean Jordan algebra consisting of only a zero element.
2634 sage: from mjo.eja.eja_algebra import TrivialEJA
2638 sage: J = TrivialEJA()
2645 sage: 7*J.one()*12*J.one()
2647 sage: J.one().inner_product(J.one())
2649 sage: J.one().norm()
2651 sage: J.one().subalgebra_generated_by()
2652 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2657 def __init__(self
, **kwargs
):
2658 jordan_product
= lambda x
,y
: x
2659 inner_product
= lambda x
,y
: 0
2662 # New defaults for keyword arguments
2663 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2664 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2666 super(TrivialEJA
, self
).__init
__(basis
,
2670 # The rank is zero using my definition, namely the dimension of the
2671 # largest subalgebra generated by any element.
2672 self
.rank
.set_cache(0)
2673 self
.one
.set_cache( self
.zero() )
2676 def random_instance(cls
, **kwargs
):
2677 # We don't take a "size" argument so the superclass method is
2678 # inappropriate for us.
2679 return cls(**kwargs
)
2682 class CartesianProductEJA(CombinatorialFreeModule_CartesianProduct
,
2683 FiniteDimensionalEJA
):
2685 The external (orthogonal) direct sum of two or more Euclidean
2686 Jordan algebras. Every Euclidean Jordan algebra decomposes into an
2687 orthogonal direct sum of simple Euclidean Jordan algebras which is
2688 then isometric to a Cartesian product, so no generality is lost by
2689 providing only this construction.
2693 sage: from mjo.eja.eja_algebra import (CartesianProductEJA,
2695 ....: JordanSpinEJA,
2696 ....: RealSymmetricEJA)
2700 The Jordan product is inherited from our factors and implemented by
2701 our CombinatorialFreeModule Cartesian product superclass::
2703 sage: set_random_seed()
2704 sage: J1 = HadamardEJA(2)
2705 sage: J2 = RealSymmetricEJA(2)
2706 sage: J = cartesian_product([J1,J2])
2707 sage: x,y = J.random_elements(2)
2711 The ability to retrieve the original factors is implemented by our
2712 CombinatorialFreeModule Cartesian product superclass::
2714 sage: J1 = HadamardEJA(2, field=QQ)
2715 sage: J2 = JordanSpinEJA(3, field=QQ)
2716 sage: J = cartesian_product([J1,J2])
2717 sage: J.cartesian_factors()
2718 (Euclidean Jordan algebra of dimension 2 over Rational Field,
2719 Euclidean Jordan algebra of dimension 3 over Rational Field)
2721 You can provide more than two factors::
2723 sage: J1 = HadamardEJA(2)
2724 sage: J2 = JordanSpinEJA(3)
2725 sage: J3 = RealSymmetricEJA(3)
2726 sage: cartesian_product([J1,J2,J3])
2727 Euclidean Jordan algebra of dimension 2 over Algebraic Real
2728 Field (+) Euclidean Jordan algebra of dimension 3 over Algebraic
2729 Real Field (+) Euclidean Jordan algebra of dimension 6 over
2730 Algebraic Real Field
2734 All factors must share the same base field::
2736 sage: J1 = HadamardEJA(2, field=QQ)
2737 sage: J2 = RealSymmetricEJA(2)
2738 sage: CartesianProductEJA((J1,J2))
2739 Traceback (most recent call last):
2741 ValueError: all factors must share the same base field
2744 def __init__(self
, modules
, **kwargs
):
2745 CombinatorialFreeModule_CartesianProduct
.__init
__(self
,
2748 field
= modules
[0].base_ring()
2749 if not all( J
.base_ring() == field
for J
in modules
):
2750 raise ValueError("all factors must share the same base field")
2752 basis
= tuple( b
.to_vector().column() for b
in self
.basis() )
2754 # Define jordan/inner products that operate on the basis.
2755 def jordan_product(x_mat
,y_mat
):
2756 x
= self
.from_vector(_mat2vec(x_mat
))
2757 y
= self
.from_vector(_mat2vec(y_mat
))
2758 return self
.cartesian_jordan_product(x
,y
).to_vector().column()
2760 def inner_product(x_mat
, y_mat
):
2761 x
= self
.from_vector(_mat2vec(x_mat
))
2762 y
= self
.from_vector(_mat2vec(y_mat
))
2763 return self
.cartesian_inner_product(x
,y
)
2765 # Use whatever category the superclass came up with. Usually
2766 # some join of the EJA and Cartesian product
2767 # categories. There's no need to check the field since it
2768 # already came from an EJA. Likewise the axioms are guaranteed
2769 # to be satisfied. We can't orthonormalize by default because
2770 # there's no way to pass "orthonormalize=False" to
2771 # cartesian_product(...) when the base ring is QQ and
2772 # orthonormalizing would give us irrational entries.
2774 # TODO: create a separate constructor that is capable of
2775 # orthonormalizing and is only used by the cartesian_product()
2777 FiniteDimensionalEJA
.__init
__(self
,
2782 orthonormalize
=False,
2785 category
=self
.category())
2788 # Initialize the FDEJA class, too. Does this override the
2789 # initialization that we did for the
2790 # CombinatorialFreeModule_CartesianProduct class? If not, we
2791 # will probably have to duplicate some of the work (i.e. one
2792 # of the constructors). Since the CartesianProduct one is
2793 # smaller, that makes the most sense to copy/paste if it comes
2797 self
.rank
.set_cache(sum(J
.rank() for J
in modules
))
2800 def cartesian_projection(self
, i
):
2804 sage: from mjo.eja.eja_algebra import (random_eja,
2805 ....: JordanSpinEJA,
2807 ....: RealSymmetricEJA,
2808 ....: ComplexHermitianEJA)
2812 The projection morphisms are Euclidean Jordan algebra
2815 sage: J1 = HadamardEJA(2)
2816 sage: J2 = RealSymmetricEJA(2)
2817 sage: J = cartesian_product([J1,J2])
2818 sage: J.cartesian_projection(0)
2819 Linear operator between finite-dimensional Euclidean Jordan
2820 algebras represented by the matrix:
2823 Domain: Euclidean Jordan algebra of dimension 2 over Algebraic
2824 Real Field (+) Euclidean Jordan algebra of dimension 3 over
2825 Algebraic Real Field
2826 Codomain: Euclidean Jordan algebra of dimension 2 over Algebraic
2828 sage: J.cartesian_projection(1)
2829 Linear operator between finite-dimensional Euclidean Jordan
2830 algebras represented by the matrix:
2834 Domain: Euclidean Jordan algebra of dimension 2 over Algebraic
2835 Real Field (+) Euclidean Jordan algebra of dimension 3 over
2836 Algebraic Real Field
2837 Codomain: Euclidean Jordan algebra of dimension 3 over Algebraic
2840 The projections work the way you'd expect on the vector
2841 representation of an element::
2843 sage: J1 = JordanSpinEJA(2)
2844 sage: J2 = ComplexHermitianEJA(2)
2845 sage: J = cartesian_product([J1,J2])
2846 sage: pi_left = J.cartesian_projection(0)
2847 sage: pi_right = J.cartesian_projection(1)
2848 sage: pi_left(J.one()).to_vector()
2850 sage: pi_right(J.one()).to_vector()
2852 sage: J.one().to_vector()
2857 The answer never changes::
2859 sage: set_random_seed()
2860 sage: J1 = random_eja()
2861 sage: J2 = random_eja()
2862 sage: J = cartesian_product([J1,J2])
2863 sage: P0 = J.cartesian_projection(0)
2864 sage: P1 = J.cartesian_projection(0)
2869 Ji
= self
.cartesian_factors()[i
]
2870 # Requires the fix on Trac 31421/31422 to work!
2871 Pi
= super().cartesian_projection(i
)
2872 return FiniteDimensionalEJAOperator(self
,Ji
,Pi
.matrix())
2875 def cartesian_embedding(self
, i
):
2879 sage: from mjo.eja.eja_algebra import (random_eja,
2880 ....: JordanSpinEJA,
2882 ....: RealSymmetricEJA)
2886 The embedding morphisms are Euclidean Jordan algebra
2889 sage: J1 = HadamardEJA(2)
2890 sage: J2 = RealSymmetricEJA(2)
2891 sage: J = cartesian_product([J1,J2])
2892 sage: J.cartesian_embedding(0)
2893 Linear operator between finite-dimensional Euclidean Jordan
2894 algebras represented by the matrix:
2900 Domain: Euclidean Jordan algebra of dimension 2 over
2901 Algebraic Real Field
2902 Codomain: Euclidean Jordan algebra of dimension 2 over
2903 Algebraic Real Field (+) Euclidean Jordan algebra of
2904 dimension 3 over Algebraic Real Field
2905 sage: J.cartesian_embedding(1)
2906 Linear operator between finite-dimensional Euclidean Jordan
2907 algebras represented by the matrix:
2913 Domain: Euclidean Jordan algebra of dimension 3 over
2914 Algebraic Real Field
2915 Codomain: Euclidean Jordan algebra of dimension 2 over
2916 Algebraic Real Field (+) Euclidean Jordan algebra of
2917 dimension 3 over Algebraic Real Field
2919 The embeddings work the way you'd expect on the vector
2920 representation of an element::
2922 sage: J1 = JordanSpinEJA(3)
2923 sage: J2 = RealSymmetricEJA(2)
2924 sage: J = cartesian_product([J1,J2])
2925 sage: iota_left = J.cartesian_embedding(0)
2926 sage: iota_right = J.cartesian_embedding(1)
2927 sage: iota_left(J1.zero()) == J.zero()
2929 sage: iota_right(J2.zero()) == J.zero()
2931 sage: J1.one().to_vector()
2933 sage: iota_left(J1.one()).to_vector()
2935 sage: J2.one().to_vector()
2937 sage: iota_right(J2.one()).to_vector()
2939 sage: J.one().to_vector()
2944 The answer never changes::
2946 sage: set_random_seed()
2947 sage: J1 = random_eja()
2948 sage: J2 = random_eja()
2949 sage: J = cartesian_product([J1,J2])
2950 sage: E0 = J.cartesian_embedding(0)
2951 sage: E1 = J.cartesian_embedding(0)
2955 Composing a projection with the corresponding inclusion should
2956 produce the identity map, and mismatching them should produce
2959 sage: set_random_seed()
2960 sage: J1 = random_eja()
2961 sage: J2 = random_eja()
2962 sage: J = cartesian_product([J1,J2])
2963 sage: iota_left = J.cartesian_embedding(0)
2964 sage: iota_right = J.cartesian_embedding(1)
2965 sage: pi_left = J.cartesian_projection(0)
2966 sage: pi_right = J.cartesian_projection(1)
2967 sage: pi_left*iota_left == J1.one().operator()
2969 sage: pi_right*iota_right == J2.one().operator()
2971 sage: (pi_left*iota_right).is_zero()
2973 sage: (pi_right*iota_left).is_zero()
2977 Ji
= self
.cartesian_factors()[i
]
2978 # Requires the fix on Trac 31421/31422 to work!
2979 Ei
= super().cartesian_embedding(i
)
2980 return FiniteDimensionalEJAOperator(Ji
,self
,Ei
.matrix())
2983 def cartesian_jordan_product(self
, x
, y
):
2985 The componentwise Jordan product.
2987 We project ``x`` and ``y`` onto our factors, and add up the
2988 Jordan products from the subalgebras. This may still be useful
2989 after (if) the default Jordan product in the Cartesian product
2990 algebra is overridden.
2994 sage: from mjo.eja.eja_algebra import (HadamardEJA,
2995 ....: JordanSpinEJA)
2999 sage: J1 = HadamardEJA(3)
3000 sage: J2 = JordanSpinEJA(3)
3001 sage: J = cartesian_product([J1,J2])
3002 sage: x1 = J1.from_vector(vector(QQ,(1,2,1)))
3003 sage: y1 = J1.from_vector(vector(QQ,(1,0,2)))
3004 sage: x2 = J2.from_vector(vector(QQ,(1,2,3)))
3005 sage: y2 = J2.from_vector(vector(QQ,(1,1,1)))
3006 sage: z1 = J.from_vector(vector(QQ,(1,2,1,1,2,3)))
3007 sage: z2 = J.from_vector(vector(QQ,(1,0,2,1,1,1)))
3008 sage: (x1*y1).to_vector()
3010 sage: (x2*y2).to_vector()
3012 sage: J.cartesian_jordan_product(z1,z2).to_vector()
3016 m
= len(self
.cartesian_factors())
3017 projections
= ( self
.cartesian_projection(i
) for i
in range(m
) )
3018 products
= ( P(x
)*P(y
) for P
in projections
)
3019 return self
._cartesian
_product
_of
_elements
(tuple(products
))
3021 def cartesian_inner_product(self
, x
, y
):
3023 The standard componentwise Cartesian inner-product.
3025 We project ``x`` and ``y`` onto our factors, and add up the
3026 inner-products from the subalgebras. This may still be useful
3027 after (if) the default inner product in the Cartesian product
3028 algebra is overridden.
3032 sage: from mjo.eja.eja_algebra import (HadamardEJA,
3033 ....: QuaternionHermitianEJA)
3037 sage: J1 = HadamardEJA(3,field=QQ)
3038 sage: J2 = QuaternionHermitianEJA(2,field=QQ,orthonormalize=False)
3039 sage: J = cartesian_product([J1,J2])
3044 sage: x1.inner_product(x2)
3046 sage: y1.inner_product(y2)
3048 sage: z1 = J._cartesian_product_of_elements((x1,y1))
3049 sage: z2 = J._cartesian_product_of_elements((x2,y2))
3050 sage: J.cartesian_inner_product(z1,z2)
3054 m
= len(self
.cartesian_factors())
3055 projections
= ( self
.cartesian_projection(i
) for i
in range(m
) )
3056 return sum( P(x
).inner_product(P(y
)) for P
in projections
)
3059 Element
= FiniteDimensionalEJAElement
3062 FiniteDimensionalEJA
.CartesianProduct
= CartesianProductEJA
3063 random_eja
= ConcreteEJA
.random_instance