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 (CartesianProductEJAElement
,
35 FiniteDimensionalEJAElement
)
36 from mjo
.eja
.eja_operator
import FiniteDimensionalEJAOperator
37 from mjo
.eja
.eja_utils
import _mat2vec
39 class FiniteDimensionalEJA(CombinatorialFreeModule
):
41 A finite-dimensional Euclidean Jordan algebra.
45 - basis -- a tuple of basis elements in their matrix form.
47 - jordan_product -- function of two elements (in matrix form)
48 that returns their jordan product in this algebra; this will
49 be applied to ``basis`` to compute a multiplication table for
52 - inner_product -- function of two elements (in matrix form) that
53 returns their inner product. This will be applied to ``basis`` to
54 compute an inner-product table (basically a matrix) for this algebra.
57 Element
= FiniteDimensionalEJAElement
72 if not field
.is_subring(RR
):
73 # Note: this does return true for the real algebraic
74 # field, the rationals, and any quadratic field where
75 # we've specified a real embedding.
76 raise ValueError("scalar field is not real")
78 # If the basis given to us wasn't over the field that it's
79 # supposed to be over, fix that. Or, you know, crash.
80 basis
= tuple( b
.change_ring(field
) for b
in basis
)
83 # Check commutativity of the Jordan and inner-products.
84 # This has to be done before we build the multiplication
85 # and inner-product tables/matrices, because we take
86 # advantage of symmetry in the process.
87 if not all( jordan_product(bi
,bj
) == jordan_product(bj
,bi
)
90 raise ValueError("Jordan product is not commutative")
92 if not all( inner_product(bi
,bj
) == inner_product(bj
,bi
)
95 raise ValueError("inner-product is not commutative")
99 category
= MagmaticAlgebras(field
).FiniteDimensional()
100 category
= category
.WithBasis().Unital()
102 # Element subalgebras can take advantage of this.
103 category
= category
.Associative()
105 # Call the superclass constructor so that we can use its from_vector()
106 # method to build our multiplication table.
108 super().__init
__(field
,
114 # Now comes all of the hard work. We'll be constructing an
115 # ambient vector space V that our (vectorized) basis lives in,
116 # as well as a subspace W of V spanned by those (vectorized)
117 # basis elements. The W-coordinates are the coefficients that
118 # we see in things like x = 1*e1 + 2*e2.
123 # Works on both column and square matrices...
124 degree
= len(basis
[0].list())
126 # Build an ambient space that fits our matrix basis when
127 # written out as "long vectors."
128 V
= VectorSpace(field
, degree
)
130 # The matrix that will hole the orthonormal -> unorthonormal
131 # coordinate transformation.
132 self
._deortho
_matrix
= None
135 # Save a copy of the un-orthonormalized basis for later.
136 # Convert it to ambient V (vector) coordinates while we're
137 # at it, because we'd have to do it later anyway.
138 deortho_vector_basis
= tuple( V(b
.list()) for b
in basis
)
140 from mjo
.eja
.eja_utils
import gram_schmidt
141 basis
= tuple(gram_schmidt(basis
, inner_product
))
143 # Save the (possibly orthonormalized) matrix basis for
145 self
._matrix
_basis
= basis
147 # Now create the vector space for the algebra, which will have
148 # its own set of non-ambient coordinates (in terms of the
150 vector_basis
= tuple( V(b
.list()) for b
in basis
)
151 W
= V
.span_of_basis( vector_basis
, check
=check_axioms
)
154 # Now "W" is the vector space of our algebra coordinates. The
155 # variables "X1", "X2",... refer to the entries of vectors in
156 # W. Thus to convert back and forth between the orthonormal
157 # coordinates and the given ones, we need to stick the original
159 U
= V
.span_of_basis( deortho_vector_basis
, check
=check_axioms
)
160 self
._deortho
_matrix
= matrix( U
.coordinate_vector(q
)
161 for q
in vector_basis
)
164 # Now we actually compute the multiplication and inner-product
165 # tables/matrices using the possibly-orthonormalized basis.
166 self
._inner
_product
_matrix
= matrix
.identity(field
, n
)
167 self
._multiplication
_table
= [ [0 for j
in range(i
+1)]
170 # Note: the Jordan and inner-products are defined in terms
171 # of the ambient basis. It's important that their arguments
172 # are in ambient coordinates as well.
175 # ortho basis w.r.t. ambient coords
179 # The jordan product returns a matrixy answer, so we
180 # have to convert it to the algebra coordinates.
181 elt
= jordan_product(q_i
, q_j
)
182 elt
= W
.coordinate_vector(V(elt
.list()))
183 self
._multiplication
_table
[i
][j
] = self
.from_vector(elt
)
185 if not orthonormalize
:
186 # If we're orthonormalizing the basis with respect
187 # to an inner-product, then the inner-product
188 # matrix with respect to the resulting basis is
189 # just going to be the identity.
190 ip
= inner_product(q_i
, q_j
)
191 self
._inner
_product
_matrix
[i
,j
] = ip
192 self
._inner
_product
_matrix
[j
,i
] = ip
194 self
._inner
_product
_matrix
._cache
= {'hermitian': True}
195 self
._inner
_product
_matrix
.set_immutable()
198 if not self
._is
_jordanian
():
199 raise ValueError("Jordan identity does not hold")
200 if not self
._inner
_product
_is
_associative
():
201 raise ValueError("inner product is not associative")
204 def _coerce_map_from_base_ring(self
):
206 Disable the map from the base ring into the algebra.
208 Performing a nonsense conversion like this automatically
209 is counterpedagogical. The fallback is to try the usual
210 element constructor, which should also fail.
214 sage: from mjo.eja.eja_algebra import random_eja
218 sage: set_random_seed()
219 sage: J = random_eja()
221 Traceback (most recent call last):
223 ValueError: not an element of this algebra
229 def product_on_basis(self
, i
, j
):
230 # We only stored the lower-triangular portion of the
231 # multiplication table.
233 return self
._multiplication
_table
[i
][j
]
235 return self
._multiplication
_table
[j
][i
]
237 def inner_product(self
, x
, y
):
239 The inner product associated with this Euclidean Jordan algebra.
241 Defaults to the trace inner product, but can be overridden by
242 subclasses if they are sure that the necessary properties are
247 sage: from mjo.eja.eja_algebra import (random_eja,
249 ....: BilinearFormEJA)
253 Our inner product is "associative," which means the following for
254 a symmetric bilinear form::
256 sage: set_random_seed()
257 sage: J = random_eja()
258 sage: x,y,z = J.random_elements(3)
259 sage: (x*y).inner_product(z) == y.inner_product(x*z)
264 Ensure that this is the usual inner product for the algebras
267 sage: set_random_seed()
268 sage: J = HadamardEJA.random_instance()
269 sage: x,y = J.random_elements(2)
270 sage: actual = x.inner_product(y)
271 sage: expected = x.to_vector().inner_product(y.to_vector())
272 sage: actual == expected
275 Ensure that this is one-half of the trace inner-product in a
276 BilinearFormEJA that isn't just the reals (when ``n`` isn't
277 one). This is in Faraut and Koranyi, and also my "On the
280 sage: set_random_seed()
281 sage: J = BilinearFormEJA.random_instance()
282 sage: n = J.dimension()
283 sage: x = J.random_element()
284 sage: y = J.random_element()
285 sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
288 B
= self
._inner
_product
_matrix
289 return (B
*x
.to_vector()).inner_product(y
.to_vector())
292 def _is_commutative(self
):
294 Whether or not this algebra's multiplication table is commutative.
296 This method should of course always return ``True``, unless
297 this algebra was constructed with ``check_axioms=False`` and
298 passed an invalid multiplication table.
300 return all( self
.product_on_basis(i
,j
) == self
.product_on_basis(i
,j
)
301 for i
in range(self
.dimension())
302 for j
in range(self
.dimension()) )
304 def _is_jordanian(self
):
306 Whether or not this algebra's multiplication table respects the
307 Jordan identity `(x^{2})(xy) = x(x^{2}y)`.
309 We only check one arrangement of `x` and `y`, so for a
310 ``True`` result to be truly true, you should also check
311 :meth:`_is_commutative`. This method should of course always
312 return ``True``, unless this algebra was constructed with
313 ``check_axioms=False`` and passed an invalid multiplication table.
315 return all( (self
.gens()[i
]**2)*(self
.gens()[i
]*self
.gens()[j
])
317 (self
.gens()[i
])*((self
.gens()[i
]**2)*self
.gens()[j
])
318 for i
in range(self
.dimension())
319 for j
in range(self
.dimension()) )
321 def _inner_product_is_associative(self
):
323 Return whether or not this algebra's inner product `B` is
324 associative; that is, whether or not `B(xy,z) = B(x,yz)`.
326 This method should of course always return ``True``, unless
327 this algebra was constructed with ``check_axioms=False`` and
328 passed an invalid Jordan or inner-product.
331 # Used to check whether or not something is zero in an inexact
332 # ring. This number is sufficient to allow the construction of
333 # QuaternionHermitianEJA(2, field=RDF) with check_axioms=True.
336 for i
in range(self
.dimension()):
337 for j
in range(self
.dimension()):
338 for k
in range(self
.dimension()):
342 diff
= (x
*y
).inner_product(z
) - x
.inner_product(y
*z
)
344 if self
.base_ring().is_exact():
348 if diff
.abs() > epsilon
:
353 def _element_constructor_(self
, elt
):
355 Construct an element of this algebra from its vector or matrix
358 This gets called only after the parent element _call_ method
359 fails to find a coercion for the argument.
363 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
365 ....: RealSymmetricEJA)
369 The identity in `S^n` is converted to the identity in the EJA::
371 sage: J = RealSymmetricEJA(3)
372 sage: I = matrix.identity(QQ,3)
373 sage: J(I) == J.one()
376 This skew-symmetric matrix can't be represented in the EJA::
378 sage: J = RealSymmetricEJA(3)
379 sage: A = matrix(QQ,3, lambda i,j: i-j)
381 Traceback (most recent call last):
383 ValueError: not an element of this algebra
387 Ensure that we can convert any element of the two non-matrix
388 simple algebras (whose matrix representations are columns)
389 back and forth faithfully::
391 sage: set_random_seed()
392 sage: J = HadamardEJA.random_instance()
393 sage: x = J.random_element()
394 sage: J(x.to_vector().column()) == x
396 sage: J = JordanSpinEJA.random_instance()
397 sage: x = J.random_element()
398 sage: J(x.to_vector().column()) == x
402 msg
= "not an element of this algebra"
404 # The superclass implementation of random_element()
405 # needs to be able to coerce "0" into the algebra.
407 elif elt
in self
.base_ring():
408 # Ensure that no base ring -> algebra coercion is performed
409 # by this method. There's some stupidity in sage that would
410 # otherwise propagate to this method; for example, sage thinks
411 # that the integer 3 belongs to the space of 2-by-2 matrices.
412 raise ValueError(msg
)
416 except (AttributeError, TypeError):
417 # Try to convert a vector into a column-matrix
420 if elt
not in self
.matrix_space():
421 raise ValueError(msg
)
423 # Thanks for nothing! Matrix spaces aren't vector spaces in
424 # Sage, so we have to figure out its matrix-basis coordinates
425 # ourselves. We use the basis space's ring instead of the
426 # element's ring because the basis space might be an algebraic
427 # closure whereas the base ring of the 3-by-3 identity matrix
428 # could be QQ instead of QQbar.
430 # We pass check=False because the matrix basis is "guaranteed"
431 # to be linearly independent... right? Ha ha.
432 V
= VectorSpace(self
.base_ring(), elt
.nrows()*elt
.ncols())
433 W
= V
.span_of_basis( (_mat2vec(s
) for s
in self
.matrix_basis()),
437 coords
= W
.coordinate_vector(_mat2vec(elt
))
438 except ArithmeticError: # vector is not in free module
439 raise ValueError(msg
)
441 return self
.from_vector(coords
)
445 Return a string representation of ``self``.
449 sage: from mjo.eja.eja_algebra import JordanSpinEJA
453 Ensure that it says what we think it says::
455 sage: JordanSpinEJA(2, field=AA)
456 Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
457 sage: JordanSpinEJA(3, field=RDF)
458 Euclidean Jordan algebra of dimension 3 over Real Double Field
461 fmt
= "Euclidean Jordan algebra of dimension {} over {}"
462 return fmt
.format(self
.dimension(), self
.base_ring())
466 def characteristic_polynomial_of(self
):
468 Return the algebra's "characteristic polynomial of" function,
469 which is itself a multivariate polynomial that, when evaluated
470 at the coordinates of some algebra element, returns that
471 element's characteristic polynomial.
473 The resulting polynomial has `n+1` variables, where `n` is the
474 dimension of this algebra. The first `n` variables correspond to
475 the coordinates of an algebra element: when evaluated at the
476 coordinates of an algebra element with respect to a certain
477 basis, the result is a univariate polynomial (in the one
478 remaining variable ``t``), namely the characteristic polynomial
483 sage: from mjo.eja.eja_algebra import JordanSpinEJA, TrivialEJA
487 The characteristic polynomial in the spin algebra is given in
488 Alizadeh, Example 11.11::
490 sage: J = JordanSpinEJA(3)
491 sage: p = J.characteristic_polynomial_of(); p
492 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
493 sage: xvec = J.one().to_vector()
497 By definition, the characteristic polynomial is a monic
498 degree-zero polynomial in a rank-zero algebra. Note that
499 Cayley-Hamilton is indeed satisfied since the polynomial
500 ``1`` evaluates to the identity element of the algebra on
503 sage: J = TrivialEJA()
504 sage: J.characteristic_polynomial_of()
511 # The list of coefficient polynomials a_0, a_1, a_2, ..., a_(r-1).
512 a
= self
._charpoly
_coefficients
()
514 # We go to a bit of trouble here to reorder the
515 # indeterminates, so that it's easier to evaluate the
516 # characteristic polynomial at x's coordinates and get back
517 # something in terms of t, which is what we want.
518 S
= PolynomialRing(self
.base_ring(),'t')
522 S
= PolynomialRing(S
, R
.variable_names())
525 return (t
**r
+ sum( a
[k
]*(t
**k
) for k
in range(r
) ))
527 def coordinate_polynomial_ring(self
):
529 The multivariate polynomial ring in which this algebra's
530 :meth:`characteristic_polynomial_of` lives.
534 sage: from mjo.eja.eja_algebra import (HadamardEJA,
535 ....: RealSymmetricEJA)
539 sage: J = HadamardEJA(2)
540 sage: J.coordinate_polynomial_ring()
541 Multivariate Polynomial Ring in X1, X2...
542 sage: J = RealSymmetricEJA(3,field=QQ,orthonormalize=False)
543 sage: J.coordinate_polynomial_ring()
544 Multivariate Polynomial Ring in X1, X2, X3, X4, X5, X6...
547 var_names
= tuple( "X%d" % z
for z
in range(1, self
.dimension()+1) )
548 return PolynomialRing(self
.base_ring(), var_names
)
550 def inner_product(self
, x
, y
):
552 The inner product associated with this Euclidean Jordan algebra.
554 Defaults to the trace inner product, but can be overridden by
555 subclasses if they are sure that the necessary properties are
560 sage: from mjo.eja.eja_algebra import (random_eja,
562 ....: BilinearFormEJA)
566 Our inner product is "associative," which means the following for
567 a symmetric bilinear form::
569 sage: set_random_seed()
570 sage: J = random_eja()
571 sage: x,y,z = J.random_elements(3)
572 sage: (x*y).inner_product(z) == y.inner_product(x*z)
577 Ensure that this is the usual inner product for the algebras
580 sage: set_random_seed()
581 sage: J = HadamardEJA.random_instance()
582 sage: x,y = J.random_elements(2)
583 sage: actual = x.inner_product(y)
584 sage: expected = x.to_vector().inner_product(y.to_vector())
585 sage: actual == expected
588 Ensure that this is one-half of the trace inner-product in a
589 BilinearFormEJA that isn't just the reals (when ``n`` isn't
590 one). This is in Faraut and Koranyi, and also my "On the
593 sage: set_random_seed()
594 sage: J = BilinearFormEJA.random_instance()
595 sage: n = J.dimension()
596 sage: x = J.random_element()
597 sage: y = J.random_element()
598 sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
601 B
= self
._inner
_product
_matrix
602 return (B
*x
.to_vector()).inner_product(y
.to_vector())
605 def is_trivial(self
):
607 Return whether or not this algebra is trivial.
609 A trivial algebra contains only the zero element.
613 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
618 sage: J = ComplexHermitianEJA(3)
624 sage: J = TrivialEJA()
629 return self
.dimension() == 0
632 def multiplication_table(self
):
634 Return a visual representation of this algebra's multiplication
635 table (on basis elements).
639 sage: from mjo.eja.eja_algebra import JordanSpinEJA
643 sage: J = JordanSpinEJA(4)
644 sage: J.multiplication_table()
645 +----++----+----+----+----+
646 | * || e0 | e1 | e2 | e3 |
647 +====++====+====+====+====+
648 | e0 || e0 | e1 | e2 | e3 |
649 +----++----+----+----+----+
650 | e1 || e1 | e0 | 0 | 0 |
651 +----++----+----+----+----+
652 | e2 || e2 | 0 | e0 | 0 |
653 +----++----+----+----+----+
654 | e3 || e3 | 0 | 0 | e0 |
655 +----++----+----+----+----+
659 # Prepend the header row.
660 M
= [["*"] + list(self
.gens())]
662 # And to each subsequent row, prepend an entry that belongs to
663 # the left-side "header column."
664 M
+= [ [self
.gens()[i
]] + [ self
.product_on_basis(i
,j
)
668 return table(M
, header_row
=True, header_column
=True, frame
=True)
671 def matrix_basis(self
):
673 Return an (often more natural) representation of this algebras
674 basis as an ordered tuple of matrices.
676 Every finite-dimensional Euclidean Jordan Algebra is a, up to
677 Jordan isomorphism, a direct sum of five simple
678 algebras---four of which comprise Hermitian matrices. And the
679 last type of algebra can of course be thought of as `n`-by-`1`
680 column matrices (ambiguusly called column vectors) to avoid
681 special cases. As a result, matrices (and column vectors) are
682 a natural representation format for Euclidean Jordan algebra
685 But, when we construct an algebra from a basis of matrices,
686 those matrix representations are lost in favor of coordinate
687 vectors *with respect to* that basis. We could eventually
688 convert back if we tried hard enough, but having the original
689 representations handy is valuable enough that we simply store
690 them and return them from this method.
692 Why implement this for non-matrix algebras? Avoiding special
693 cases for the :class:`BilinearFormEJA` pays with simplicity in
694 its own right. But mainly, we would like to be able to assume
695 that elements of a :class:`CartesianProductEJA` can be displayed
696 nicely, without having to have special classes for direct sums
697 one of whose components was a matrix algebra.
701 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
702 ....: RealSymmetricEJA)
706 sage: J = RealSymmetricEJA(2)
708 Finite family {0: e0, 1: e1, 2: e2}
709 sage: J.matrix_basis()
711 [1 0] [ 0 0.7071067811865475?] [0 0]
712 [0 0], [0.7071067811865475? 0], [0 1]
717 sage: J = JordanSpinEJA(2)
719 Finite family {0: e0, 1: e1}
720 sage: J.matrix_basis()
726 return self
._matrix
_basis
729 def matrix_space(self
):
731 Return the matrix space in which this algebra's elements live, if
732 we think of them as matrices (including column vectors of the
735 Generally this will be an `n`-by-`1` column-vector space,
736 except when the algebra is trivial. There it's `n`-by-`n`
737 (where `n` is zero), to ensure that two elements of the matrix
738 space (empty matrices) can be multiplied.
740 Matrix algebras override this with something more useful.
742 if self
.is_trivial():
743 return MatrixSpace(self
.base_ring(), 0)
745 return self
.matrix_basis()[0].parent()
751 Return the unit element of this algebra.
755 sage: from mjo.eja.eja_algebra import (HadamardEJA,
760 We can compute unit element in the Hadamard EJA::
762 sage: J = HadamardEJA(5)
764 e0 + e1 + e2 + e3 + e4
766 The unit element in the Hadamard EJA is inherited in the
767 subalgebras generated by its elements::
769 sage: J = HadamardEJA(5)
771 e0 + e1 + e2 + e3 + e4
772 sage: x = sum(J.gens())
773 sage: A = x.subalgebra_generated_by(orthonormalize=False)
776 sage: A.one().superalgebra_element()
777 e0 + e1 + e2 + e3 + e4
781 The identity element acts like the identity, regardless of
782 whether or not we orthonormalize::
784 sage: set_random_seed()
785 sage: J = random_eja()
786 sage: x = J.random_element()
787 sage: J.one()*x == x and x*J.one() == x
789 sage: A = x.subalgebra_generated_by()
790 sage: y = A.random_element()
791 sage: A.one()*y == y and y*A.one() == y
796 sage: set_random_seed()
797 sage: J = random_eja(field=QQ, orthonormalize=False)
798 sage: x = J.random_element()
799 sage: J.one()*x == x and x*J.one() == x
801 sage: A = x.subalgebra_generated_by(orthonormalize=False)
802 sage: y = A.random_element()
803 sage: A.one()*y == y and y*A.one() == y
806 The matrix of the unit element's operator is the identity,
807 regardless of the base field and whether or not we
810 sage: set_random_seed()
811 sage: J = random_eja()
812 sage: actual = J.one().operator().matrix()
813 sage: expected = matrix.identity(J.base_ring(), J.dimension())
814 sage: actual == expected
816 sage: x = J.random_element()
817 sage: A = x.subalgebra_generated_by()
818 sage: actual = A.one().operator().matrix()
819 sage: expected = matrix.identity(A.base_ring(), A.dimension())
820 sage: actual == expected
825 sage: set_random_seed()
826 sage: J = random_eja(field=QQ, orthonormalize=False)
827 sage: actual = J.one().operator().matrix()
828 sage: expected = matrix.identity(J.base_ring(), J.dimension())
829 sage: actual == expected
831 sage: x = J.random_element()
832 sage: A = x.subalgebra_generated_by(orthonormalize=False)
833 sage: actual = A.one().operator().matrix()
834 sage: expected = matrix.identity(A.base_ring(), A.dimension())
835 sage: actual == expected
838 Ensure that the cached unit element (often precomputed by
839 hand) agrees with the computed one::
841 sage: set_random_seed()
842 sage: J = random_eja()
843 sage: cached = J.one()
844 sage: J.one.clear_cache()
845 sage: J.one() == cached
850 sage: set_random_seed()
851 sage: J = random_eja(field=QQ, orthonormalize=False)
852 sage: cached = J.one()
853 sage: J.one.clear_cache()
854 sage: J.one() == cached
858 # We can brute-force compute the matrices of the operators
859 # that correspond to the basis elements of this algebra.
860 # If some linear combination of those basis elements is the
861 # algebra identity, then the same linear combination of
862 # their matrices has to be the identity matrix.
864 # Of course, matrices aren't vectors in sage, so we have to
865 # appeal to the "long vectors" isometry.
866 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
868 # Now we use basic linear algebra to find the coefficients,
869 # of the matrices-as-vectors-linear-combination, which should
870 # work for the original algebra basis too.
871 A
= matrix(self
.base_ring(), oper_vecs
)
873 # We used the isometry on the left-hand side already, but we
874 # still need to do it for the right-hand side. Recall that we
875 # wanted something that summed to the identity matrix.
876 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
878 # Now if there's an identity element in the algebra, this
879 # should work. We solve on the left to avoid having to
880 # transpose the matrix "A".
881 return self
.from_vector(A
.solve_left(b
))
884 def peirce_decomposition(self
, c
):
886 The Peirce decomposition of this algebra relative to the
889 In the future, this can be extended to a complete system of
890 orthogonal idempotents.
894 - ``c`` -- an idempotent of this algebra.
898 A triple (J0, J5, J1) containing two subalgebras and one subspace
901 - ``J0`` -- the algebra on the eigenspace of ``c.operator()``
902 corresponding to the eigenvalue zero.
904 - ``J5`` -- the eigenspace (NOT a subalgebra) of ``c.operator()``
905 corresponding to the eigenvalue one-half.
907 - ``J1`` -- the algebra on the eigenspace of ``c.operator()``
908 corresponding to the eigenvalue one.
910 These are the only possible eigenspaces for that operator, and this
911 algebra is a direct sum of them. The spaces ``J0`` and ``J1`` are
912 orthogonal, and are subalgebras of this algebra with the appropriate
917 sage: from mjo.eja.eja_algebra import random_eja, RealSymmetricEJA
921 The canonical example comes from the symmetric matrices, which
922 decompose into diagonal and off-diagonal parts::
924 sage: J = RealSymmetricEJA(3)
925 sage: C = matrix(QQ, [ [1,0,0],
929 sage: J0,J5,J1 = J.peirce_decomposition(c)
931 Euclidean Jordan algebra of dimension 1...
933 Vector space of degree 6 and dimension 2...
935 Euclidean Jordan algebra of dimension 3...
936 sage: J0.one().to_matrix()
940 sage: orig_df = AA.options.display_format
941 sage: AA.options.display_format = 'radical'
942 sage: J.from_vector(J5.basis()[0]).to_matrix()
946 sage: J.from_vector(J5.basis()[1]).to_matrix()
950 sage: AA.options.display_format = orig_df
951 sage: J1.one().to_matrix()
958 Every algebra decomposes trivially with respect to its identity
961 sage: set_random_seed()
962 sage: J = random_eja()
963 sage: J0,J5,J1 = J.peirce_decomposition(J.one())
964 sage: J0.dimension() == 0 and J5.dimension() == 0
966 sage: J1.superalgebra() == J and J1.dimension() == J.dimension()
969 The decomposition is into eigenspaces, and its components are
970 therefore necessarily orthogonal. Moreover, the identity
971 elements in the two subalgebras are the projections onto their
972 respective subspaces of the superalgebra's identity element::
974 sage: set_random_seed()
975 sage: J = random_eja()
976 sage: x = J.random_element()
977 sage: if not J.is_trivial():
978 ....: while x.is_nilpotent():
979 ....: x = J.random_element()
980 sage: c = x.subalgebra_idempotent()
981 sage: J0,J5,J1 = J.peirce_decomposition(c)
983 sage: for (w,y,z) in zip(J0.basis(), J5.basis(), J1.basis()):
984 ....: w = w.superalgebra_element()
985 ....: y = J.from_vector(y)
986 ....: z = z.superalgebra_element()
987 ....: ipsum += w.inner_product(y).abs()
988 ....: ipsum += w.inner_product(z).abs()
989 ....: ipsum += y.inner_product(z).abs()
992 sage: J1(c) == J1.one()
994 sage: J0(J.one() - c) == J0.one()
998 if not c
.is_idempotent():
999 raise ValueError("element is not idempotent: %s" % c
)
1001 from mjo
.eja
.eja_subalgebra
import FiniteDimensionalEJASubalgebra
1003 # Default these to what they should be if they turn out to be
1004 # trivial, because eigenspaces_left() won't return eigenvalues
1005 # corresponding to trivial spaces (e.g. it returns only the
1006 # eigenspace corresponding to lambda=1 if you take the
1007 # decomposition relative to the identity element).
1008 trivial
= FiniteDimensionalEJASubalgebra(self
, ())
1009 J0
= trivial
# eigenvalue zero
1010 J5
= VectorSpace(self
.base_ring(), 0) # eigenvalue one-half
1011 J1
= trivial
# eigenvalue one
1013 for (eigval
, eigspace
) in c
.operator().matrix().right_eigenspaces():
1014 if eigval
== ~
(self
.base_ring()(2)):
1017 gens
= tuple( self
.from_vector(b
) for b
in eigspace
.basis() )
1018 subalg
= FiniteDimensionalEJASubalgebra(self
,
1026 raise ValueError("unexpected eigenvalue: %s" % eigval
)
1031 def random_element(self
, thorough
=False):
1033 Return a random element of this algebra.
1035 Our algebra superclass method only returns a linear
1036 combination of at most two basis elements. We instead
1037 want the vector space "random element" method that
1038 returns a more diverse selection.
1042 - ``thorough`` -- (boolean; default False) whether or not we
1043 should generate irrational coefficients for the random
1044 element when our base ring is irrational; this slows the
1045 algebra operations to a crawl, but any truly random method
1049 # For a general base ring... maybe we can trust this to do the
1050 # right thing? Unlikely, but.
1051 V
= self
.vector_space()
1052 v
= V
.random_element()
1054 if self
.base_ring() is AA
:
1055 # The "random element" method of the algebraic reals is
1056 # stupid at the moment, and only returns integers between
1057 # -2 and 2, inclusive:
1059 # https://trac.sagemath.org/ticket/30875
1061 # Instead, we implement our own "random vector" method,
1062 # and then coerce that into the algebra. We use the vector
1063 # space degree here instead of the dimension because a
1064 # subalgebra could (for example) be spanned by only two
1065 # vectors, each with five coordinates. We need to
1066 # generate all five coordinates.
1068 v
*= QQbar
.random_element().real()
1070 v
*= QQ
.random_element()
1072 return self
.from_vector(V
.coordinate_vector(v
))
1074 def random_elements(self
, count
, thorough
=False):
1076 Return ``count`` random elements as a tuple.
1080 - ``thorough`` -- (boolean; default False) whether or not we
1081 should generate irrational coefficients for the random
1082 elements when our base ring is irrational; this slows the
1083 algebra operations to a crawl, but any truly random method
1088 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1092 sage: J = JordanSpinEJA(3)
1093 sage: x,y,z = J.random_elements(3)
1094 sage: all( [ x in J, y in J, z in J ])
1096 sage: len( J.random_elements(10) ) == 10
1100 return tuple( self
.random_element(thorough
)
1101 for idx
in range(count
) )
1105 def _charpoly_coefficients(self
):
1107 The `r` polynomial coefficients of the "characteristic polynomial
1112 sage: from mjo.eja.eja_algebra import random_eja
1116 The theory shows that these are all homogeneous polynomials of
1119 sage: set_random_seed()
1120 sage: J = random_eja()
1121 sage: all(p.is_homogeneous() for p in J._charpoly_coefficients())
1125 n
= self
.dimension()
1126 R
= self
.coordinate_polynomial_ring()
1128 F
= R
.fraction_field()
1131 # From a result in my book, these are the entries of the
1132 # basis representation of L_x.
1133 return sum( vars[k
]*self
.gens()[k
].operator().matrix()[i
,j
]
1136 L_x
= matrix(F
, n
, n
, L_x_i_j
)
1139 if self
.rank
.is_in_cache():
1141 # There's no need to pad the system with redundant
1142 # columns if we *know* they'll be redundant.
1145 # Compute an extra power in case the rank is equal to
1146 # the dimension (otherwise, we would stop at x^(r-1)).
1147 x_powers
= [ (L_x
**k
)*self
.one().to_vector()
1148 for k
in range(n
+1) ]
1149 A
= matrix
.column(F
, x_powers
[:n
])
1150 AE
= A
.extended_echelon_form()
1157 # The theory says that only the first "r" coefficients are
1158 # nonzero, and they actually live in the original polynomial
1159 # ring and not the fraction field. We negate them because in
1160 # the actual characteristic polynomial, they get moved to the
1161 # other side where x^r lives. We don't bother to trim A_rref
1162 # down to a square matrix and solve the resulting system,
1163 # because the upper-left r-by-r portion of A_rref is
1164 # guaranteed to be the identity matrix, so e.g.
1166 # A_rref.solve_right(Y)
1168 # would just be returning Y.
1169 return (-E
*b
)[:r
].change_ring(R
)
1174 Return the rank of this EJA.
1176 This is a cached method because we know the rank a priori for
1177 all of the algebras we can construct. Thus we can avoid the
1178 expensive ``_charpoly_coefficients()`` call unless we truly
1179 need to compute the whole characteristic polynomial.
1183 sage: from mjo.eja.eja_algebra import (HadamardEJA,
1184 ....: JordanSpinEJA,
1185 ....: RealSymmetricEJA,
1186 ....: ComplexHermitianEJA,
1187 ....: QuaternionHermitianEJA,
1192 The rank of the Jordan spin algebra is always two::
1194 sage: JordanSpinEJA(2).rank()
1196 sage: JordanSpinEJA(3).rank()
1198 sage: JordanSpinEJA(4).rank()
1201 The rank of the `n`-by-`n` Hermitian real, complex, or
1202 quaternion matrices is `n`::
1204 sage: RealSymmetricEJA(4).rank()
1206 sage: ComplexHermitianEJA(3).rank()
1208 sage: QuaternionHermitianEJA(2).rank()
1213 Ensure that every EJA that we know how to construct has a
1214 positive integer rank, unless the algebra is trivial in
1215 which case its rank will be zero::
1217 sage: set_random_seed()
1218 sage: J = random_eja()
1222 sage: r > 0 or (r == 0 and J.is_trivial())
1225 Ensure that computing the rank actually works, since the ranks
1226 of all simple algebras are known and will be cached by default::
1228 sage: set_random_seed() # long time
1229 sage: J = random_eja() # long time
1230 sage: cached = J.rank() # long time
1231 sage: J.rank.clear_cache() # long time
1232 sage: J.rank() == cached # long time
1236 return len(self
._charpoly
_coefficients
())
1239 def vector_space(self
):
1241 Return the vector space that underlies this algebra.
1245 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1249 sage: J = RealSymmetricEJA(2)
1250 sage: J.vector_space()
1251 Vector space of dimension 3 over...
1254 return self
.zero().to_vector().parent().ambient_vector_space()
1257 Element
= FiniteDimensionalEJAElement
1259 class RationalBasisEJA(FiniteDimensionalEJA
):
1261 New class for algebras whose supplied basis elements have all rational entries.
1265 sage: from mjo.eja.eja_algebra import BilinearFormEJA
1269 The supplied basis is orthonormalized by default::
1271 sage: B = matrix(QQ, [[1, 0, 0], [0, 25, -32], [0, -32, 41]])
1272 sage: J = BilinearFormEJA(B)
1273 sage: J.matrix_basis()
1290 # Abuse the check_field parameter to check that the entries of
1291 # out basis (in ambient coordinates) are in the field QQ.
1292 if not all( all(b_i
in QQ
for b_i
in b
.list()) for b
in basis
):
1293 raise TypeError("basis not rational")
1295 self
._rational
_algebra
= None
1297 # There's no point in constructing the extra algebra if this
1298 # one is already rational.
1300 # Note: the same Jordan and inner-products work here,
1301 # because they are necessarily defined with respect to
1302 # ambient coordinates and not any particular basis.
1303 self
._rational
_algebra
= FiniteDimensionalEJA(
1308 orthonormalize
=False,
1312 super().__init
__(basis
,
1316 check_field
=check_field
,
1320 def _charpoly_coefficients(self
):
1324 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
1325 ....: JordanSpinEJA)
1329 The base ring of the resulting polynomial coefficients is what
1330 it should be, and not the rationals (unless the algebra was
1331 already over the rationals)::
1333 sage: J = JordanSpinEJA(3)
1334 sage: J._charpoly_coefficients()
1335 (X1^2 - X2^2 - X3^2, -2*X1)
1336 sage: a0 = J._charpoly_coefficients()[0]
1338 Algebraic Real Field
1339 sage: a0.base_ring()
1340 Algebraic Real Field
1343 if self
._rational
_algebra
is None:
1344 # There's no need to construct *another* algebra over the
1345 # rationals if this one is already over the
1346 # rationals. Likewise, if we never orthonormalized our
1347 # basis, we might as well just use the given one.
1348 return super()._charpoly
_coefficients
()
1350 # Do the computation over the rationals. The answer will be
1351 # the same, because all we've done is a change of basis.
1352 # Then, change back from QQ to our real base ring
1353 a
= ( a_i
.change_ring(self
.base_ring())
1354 for a_i
in self
._rational
_algebra
._charpoly
_coefficients
() )
1356 if self
._deortho
_matrix
is None:
1357 # This can happen if our base ring was, say, AA and we
1358 # chose not to (or didn't need to) orthonormalize. It's
1359 # still faster to do the computations over QQ even if
1360 # the numbers in the boxes stay the same.
1363 # Otherwise, convert the coordinate variables back to the
1364 # deorthonormalized ones.
1365 R
= self
.coordinate_polynomial_ring()
1366 from sage
.modules
.free_module_element
import vector
1367 X
= vector(R
, R
.gens())
1368 BX
= self
._deortho
_matrix
*X
1370 subs_dict
= { X[i]: BX[i] for i in range(len(X)) }
1371 return tuple( a_i
.subs(subs_dict
) for a_i
in a
)
1373 class ConcreteEJA(RationalBasisEJA
):
1375 A class for the Euclidean Jordan algebras that we know by name.
1377 These are the Jordan algebras whose basis, multiplication table,
1378 rank, and so on are known a priori. More to the point, they are
1379 the Euclidean Jordan algebras for which we are able to conjure up
1380 a "random instance."
1384 sage: from mjo.eja.eja_algebra import ConcreteEJA
1388 Our basis is normalized with respect to the algebra's inner
1389 product, unless we specify otherwise::
1391 sage: set_random_seed()
1392 sage: J = ConcreteEJA.random_instance()
1393 sage: all( b.norm() == 1 for b in J.gens() )
1396 Since our basis is orthonormal with respect to the algebra's inner
1397 product, and since we know that this algebra is an EJA, any
1398 left-multiplication operator's matrix will be symmetric because
1399 natural->EJA basis representation is an isometry and within the
1400 EJA the operator is self-adjoint by the Jordan axiom::
1402 sage: set_random_seed()
1403 sage: J = ConcreteEJA.random_instance()
1404 sage: x = J.random_element()
1405 sage: x.operator().is_self_adjoint()
1410 def _max_random_instance_size():
1412 Return an integer "size" that is an upper bound on the size of
1413 this algebra when it is used in a random test
1414 case. Unfortunately, the term "size" is ambiguous -- when
1415 dealing with `R^n` under either the Hadamard or Jordan spin
1416 product, the "size" refers to the dimension `n`. When dealing
1417 with a matrix algebra (real symmetric or complex/quaternion
1418 Hermitian), it refers to the size of the matrix, which is far
1419 less than the dimension of the underlying vector space.
1421 This method must be implemented in each subclass.
1423 raise NotImplementedError
1426 def random_instance(cls
, *args
, **kwargs
):
1428 Return a random instance of this type of algebra.
1430 This method should be implemented in each subclass.
1432 from sage
.misc
.prandom
import choice
1433 eja_class
= choice(cls
.__subclasses
__())
1435 # These all bubble up to the RationalBasisEJA superclass
1436 # constructor, so any (kw)args valid there are also valid
1438 return eja_class
.random_instance(*args
, **kwargs
)
1443 def dimension_over_reals():
1445 The dimension of this matrix's base ring over the reals.
1447 The reals are dimension one over themselves, obviously; that's
1448 just `\mathbb{R}^{1}`. Likewise, the complex numbers `a + bi`
1449 have dimension two. Finally, the quaternions have dimension
1450 four over the reals.
1452 This is used to determine the size of the matrix returned from
1453 :meth:`real_embed`, among other things.
1455 raise NotImplementedError
1458 def real_embed(cls
,M
):
1460 Embed the matrix ``M`` into a space of real matrices.
1462 The matrix ``M`` can have entries in any field at the moment:
1463 the real numbers, complex numbers, or quaternions. And although
1464 they are not a field, we can probably support octonions at some
1465 point, too. This function returns a real matrix that "acts like"
1466 the original with respect to matrix multiplication; i.e.
1468 real_embed(M*N) = real_embed(M)*real_embed(N)
1471 if M
.ncols() != M
.nrows():
1472 raise ValueError("the matrix 'M' must be square")
1477 def real_unembed(cls
,M
):
1479 The inverse of :meth:`real_embed`.
1481 if M
.ncols() != M
.nrows():
1482 raise ValueError("the matrix 'M' must be square")
1483 if not ZZ(M
.nrows()).mod(cls
.dimension_over_reals()).is_zero():
1484 raise ValueError("the matrix 'M' must be a real embedding")
1488 def jordan_product(X
,Y
):
1489 return (X
*Y
+ Y
*X
)/2
1492 def trace_inner_product(cls
,X
,Y
):
1494 Compute the trace inner-product of two real-embeddings.
1498 sage: from mjo.eja.eja_algebra import (RealSymmetricEJA,
1499 ....: ComplexHermitianEJA,
1500 ....: QuaternionHermitianEJA)
1504 This gives the same answer as it would if we computed the trace
1505 from the unembedded (original) matrices::
1507 sage: set_random_seed()
1508 sage: J = RealSymmetricEJA.random_instance()
1509 sage: x,y = J.random_elements(2)
1510 sage: Xe = x.to_matrix()
1511 sage: Ye = y.to_matrix()
1512 sage: X = J.real_unembed(Xe)
1513 sage: Y = J.real_unembed(Ye)
1514 sage: expected = (X*Y).trace()
1515 sage: actual = J.trace_inner_product(Xe,Ye)
1516 sage: actual == expected
1521 sage: set_random_seed()
1522 sage: J = ComplexHermitianEJA.random_instance()
1523 sage: x,y = J.random_elements(2)
1524 sage: Xe = x.to_matrix()
1525 sage: Ye = y.to_matrix()
1526 sage: X = J.real_unembed(Xe)
1527 sage: Y = J.real_unembed(Ye)
1528 sage: expected = (X*Y).trace().real()
1529 sage: actual = J.trace_inner_product(Xe,Ye)
1530 sage: actual == expected
1535 sage: set_random_seed()
1536 sage: J = QuaternionHermitianEJA.random_instance()
1537 sage: x,y = J.random_elements(2)
1538 sage: Xe = x.to_matrix()
1539 sage: Ye = y.to_matrix()
1540 sage: X = J.real_unembed(Xe)
1541 sage: Y = J.real_unembed(Ye)
1542 sage: expected = (X*Y).trace().coefficient_tuple()[0]
1543 sage: actual = J.trace_inner_product(Xe,Ye)
1544 sage: actual == expected
1548 Xu
= cls
.real_unembed(X
)
1549 Yu
= cls
.real_unembed(Y
)
1550 tr
= (Xu
*Yu
).trace()
1553 # Works in QQ, AA, RDF, et cetera.
1555 except AttributeError:
1556 # A quaternion doesn't have a real() method, but does
1557 # have coefficient_tuple() method that returns the
1558 # coefficients of 1, i, j, and k -- in that order.
1559 return tr
.coefficient_tuple()[0]
1562 class RealMatrixEJA(MatrixEJA
):
1564 def dimension_over_reals():
1568 class RealSymmetricEJA(ConcreteEJA
, RealMatrixEJA
):
1570 The rank-n simple EJA consisting of real symmetric n-by-n
1571 matrices, the usual symmetric Jordan product, and the trace inner
1572 product. It has dimension `(n^2 + n)/2` over the reals.
1576 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1580 sage: J = RealSymmetricEJA(2)
1581 sage: e0, e1, e2 = J.gens()
1589 In theory, our "field" can be any subfield of the reals::
1591 sage: RealSymmetricEJA(2, field=RDF)
1592 Euclidean Jordan algebra of dimension 3 over Real Double Field
1593 sage: RealSymmetricEJA(2, field=RR)
1594 Euclidean Jordan algebra of dimension 3 over Real Field with
1595 53 bits of precision
1599 The dimension of this algebra is `(n^2 + n) / 2`::
1601 sage: set_random_seed()
1602 sage: n_max = RealSymmetricEJA._max_random_instance_size()
1603 sage: n = ZZ.random_element(1, n_max)
1604 sage: J = RealSymmetricEJA(n)
1605 sage: J.dimension() == (n^2 + n)/2
1608 The Jordan multiplication is what we think it is::
1610 sage: set_random_seed()
1611 sage: J = RealSymmetricEJA.random_instance()
1612 sage: x,y = J.random_elements(2)
1613 sage: actual = (x*y).to_matrix()
1614 sage: X = x.to_matrix()
1615 sage: Y = y.to_matrix()
1616 sage: expected = (X*Y + Y*X)/2
1617 sage: actual == expected
1619 sage: J(expected) == x*y
1622 We can change the generator prefix::
1624 sage: RealSymmetricEJA(3, prefix='q').gens()
1625 (q0, q1, q2, q3, q4, q5)
1627 We can construct the (trivial) algebra of rank zero::
1629 sage: RealSymmetricEJA(0)
1630 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1634 def _denormalized_basis(cls
, n
):
1636 Return a basis for the space of real symmetric n-by-n matrices.
1640 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1644 sage: set_random_seed()
1645 sage: n = ZZ.random_element(1,5)
1646 sage: B = RealSymmetricEJA._denormalized_basis(n)
1647 sage: all( M.is_symmetric() for M in B)
1651 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1655 for j
in range(i
+1):
1656 Eij
= matrix(ZZ
, n
, lambda k
,l
: k
==i
and l
==j
)
1660 Sij
= Eij
+ Eij
.transpose()
1666 def _max_random_instance_size():
1667 return 4 # Dimension 10
1670 def random_instance(cls
, **kwargs
):
1672 Return a random instance of this type of algebra.
1674 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1675 return cls(n
, **kwargs
)
1677 def __init__(self
, n
, **kwargs
):
1678 # We know this is a valid EJA, but will double-check
1679 # if the user passes check_axioms=True.
1680 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
1682 super(RealSymmetricEJA
, self
).__init
__(self
._denormalized
_basis
(n
),
1683 self
.jordan_product
,
1684 self
.trace_inner_product
,
1687 # TODO: this could be factored out somehow, but is left here
1688 # because the MatrixEJA is not presently a subclass of the
1689 # FDEJA class that defines rank() and one().
1690 self
.rank
.set_cache(n
)
1691 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
1692 self
.one
.set_cache(self(idV
))
1696 class ComplexMatrixEJA(MatrixEJA
):
1697 # A manual dictionary-cache for the complex_extension() method,
1698 # since apparently @classmethods can't also be @cached_methods.
1699 _complex_extension
= {}
1702 def complex_extension(cls
,field
):
1704 The complex field that we embed/unembed, as an extension
1705 of the given ``field``.
1707 if field
in cls
._complex
_extension
:
1708 return cls
._complex
_extension
[field
]
1710 # Sage doesn't know how to adjoin the complex "i" (the root of
1711 # x^2 + 1) to a field in a general way. Here, we just enumerate
1712 # all of the cases that I have cared to support so far.
1714 # Sage doesn't know how to embed AA into QQbar, i.e. how
1715 # to adjoin sqrt(-1) to AA.
1717 elif not field
.is_exact():
1719 F
= field
.complex_field()
1721 # Works for QQ and... maybe some other fields.
1722 R
= PolynomialRing(field
, 'z')
1724 F
= field
.extension(z
**2 + 1, 'I', embedding
=CLF(-1).sqrt())
1726 cls
._complex
_extension
[field
] = F
1730 def dimension_over_reals():
1734 def real_embed(cls
,M
):
1736 Embed the n-by-n complex matrix ``M`` into the space of real
1737 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1738 bi` to the block matrix ``[[a,b],[-b,a]]``.
1742 sage: from mjo.eja.eja_algebra import ComplexMatrixEJA
1746 sage: F = QuadraticField(-1, 'I')
1747 sage: x1 = F(4 - 2*i)
1748 sage: x2 = F(1 + 2*i)
1751 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1752 sage: ComplexMatrixEJA.real_embed(M)
1761 Embedding is a homomorphism (isomorphism, in fact)::
1763 sage: set_random_seed()
1764 sage: n = ZZ.random_element(3)
1765 sage: F = QuadraticField(-1, 'I')
1766 sage: X = random_matrix(F, n)
1767 sage: Y = random_matrix(F, n)
1768 sage: Xe = ComplexMatrixEJA.real_embed(X)
1769 sage: Ye = ComplexMatrixEJA.real_embed(Y)
1770 sage: XYe = ComplexMatrixEJA.real_embed(X*Y)
1775 super(ComplexMatrixEJA
,cls
).real_embed(M
)
1778 # We don't need any adjoined elements...
1779 field
= M
.base_ring().base_ring()
1785 blocks
.append(matrix(field
, 2, [ [ a
, b
],
1788 return matrix
.block(field
, n
, blocks
)
1792 def real_unembed(cls
,M
):
1794 The inverse of _embed_complex_matrix().
1798 sage: from mjo.eja.eja_algebra import ComplexMatrixEJA
1802 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1803 ....: [-2, 1, -4, 3],
1804 ....: [ 9, 10, 11, 12],
1805 ....: [-10, 9, -12, 11] ])
1806 sage: ComplexMatrixEJA.real_unembed(A)
1808 [ 10*I + 9 12*I + 11]
1812 Unembedding is the inverse of embedding::
1814 sage: set_random_seed()
1815 sage: F = QuadraticField(-1, 'I')
1816 sage: M = random_matrix(F, 3)
1817 sage: Me = ComplexMatrixEJA.real_embed(M)
1818 sage: ComplexMatrixEJA.real_unembed(Me) == M
1822 super(ComplexMatrixEJA
,cls
).real_unembed(M
)
1824 d
= cls
.dimension_over_reals()
1825 F
= cls
.complex_extension(M
.base_ring())
1828 # Go top-left to bottom-right (reading order), converting every
1829 # 2-by-2 block we see to a single complex element.
1831 for k
in range(n
/d
):
1832 for j
in range(n
/d
):
1833 submat
= M
[d
*k
:d
*k
+d
,d
*j
:d
*j
+d
]
1834 if submat
[0,0] != submat
[1,1]:
1835 raise ValueError('bad on-diagonal submatrix')
1836 if submat
[0,1] != -submat
[1,0]:
1837 raise ValueError('bad off-diagonal submatrix')
1838 z
= submat
[0,0] + submat
[0,1]*i
1841 return matrix(F
, n
/d
, elements
)
1844 class ComplexHermitianEJA(ConcreteEJA
, ComplexMatrixEJA
):
1846 The rank-n simple EJA consisting of complex Hermitian n-by-n
1847 matrices over the real numbers, the usual symmetric Jordan product,
1848 and the real-part-of-trace inner product. It has dimension `n^2` over
1853 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1857 In theory, our "field" can be any subfield of the reals::
1859 sage: ComplexHermitianEJA(2, field=RDF)
1860 Euclidean Jordan algebra of dimension 4 over Real Double Field
1861 sage: ComplexHermitianEJA(2, field=RR)
1862 Euclidean Jordan algebra of dimension 4 over Real Field with
1863 53 bits of precision
1867 The dimension of this algebra is `n^2`::
1869 sage: set_random_seed()
1870 sage: n_max = ComplexHermitianEJA._max_random_instance_size()
1871 sage: n = ZZ.random_element(1, n_max)
1872 sage: J = ComplexHermitianEJA(n)
1873 sage: J.dimension() == n^2
1876 The Jordan multiplication is what we think it is::
1878 sage: set_random_seed()
1879 sage: J = ComplexHermitianEJA.random_instance()
1880 sage: x,y = J.random_elements(2)
1881 sage: actual = (x*y).to_matrix()
1882 sage: X = x.to_matrix()
1883 sage: Y = y.to_matrix()
1884 sage: expected = (X*Y + Y*X)/2
1885 sage: actual == expected
1887 sage: J(expected) == x*y
1890 We can change the generator prefix::
1892 sage: ComplexHermitianEJA(2, prefix='z').gens()
1895 We can construct the (trivial) algebra of rank zero::
1897 sage: ComplexHermitianEJA(0)
1898 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1903 def _denormalized_basis(cls
, n
):
1905 Returns a basis for the space of complex Hermitian n-by-n matrices.
1907 Why do we embed these? Basically, because all of numerical linear
1908 algebra assumes that you're working with vectors consisting of `n`
1909 entries from a field and scalars from the same field. There's no way
1910 to tell SageMath that (for example) the vectors contain complex
1911 numbers, while the scalar field is real.
1915 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1919 sage: set_random_seed()
1920 sage: n = ZZ.random_element(1,5)
1921 sage: B = ComplexHermitianEJA._denormalized_basis(n)
1922 sage: all( M.is_symmetric() for M in B)
1927 R
= PolynomialRing(field
, 'z')
1929 F
= field
.extension(z
**2 + 1, 'I')
1932 # This is like the symmetric case, but we need to be careful:
1934 # * We want conjugate-symmetry, not just symmetry.
1935 # * The diagonal will (as a result) be real.
1938 Eij
= matrix
.zero(F
,n
)
1940 for j
in range(i
+1):
1944 Sij
= cls
.real_embed(Eij
)
1947 # The second one has a minus because it's conjugated.
1948 Eij
[j
,i
] = 1 # Eij = Eij + Eij.transpose()
1949 Sij_real
= cls
.real_embed(Eij
)
1951 # Eij = I*Eij - I*Eij.transpose()
1954 Sij_imag
= cls
.real_embed(Eij
)
1960 # Since we embedded these, we can drop back to the "field" that we
1961 # started with instead of the complex extension "F".
1962 return tuple( s
.change_ring(field
) for s
in S
)
1965 def __init__(self
, n
, **kwargs
):
1966 # We know this is a valid EJA, but will double-check
1967 # if the user passes check_axioms=True.
1968 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
1970 super(ComplexHermitianEJA
, self
).__init
__(self
._denormalized
_basis
(n
),
1971 self
.jordan_product
,
1972 self
.trace_inner_product
,
1974 # TODO: this could be factored out somehow, but is left here
1975 # because the MatrixEJA is not presently a subclass of the
1976 # FDEJA class that defines rank() and one().
1977 self
.rank
.set_cache(n
)
1978 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
1979 self
.one
.set_cache(self(idV
))
1982 def _max_random_instance_size():
1983 return 3 # Dimension 9
1986 def random_instance(cls
, **kwargs
):
1988 Return a random instance of this type of algebra.
1990 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1991 return cls(n
, **kwargs
)
1993 class QuaternionMatrixEJA(MatrixEJA
):
1995 # A manual dictionary-cache for the quaternion_extension() method,
1996 # since apparently @classmethods can't also be @cached_methods.
1997 _quaternion_extension
= {}
2000 def quaternion_extension(cls
,field
):
2002 The quaternion field that we embed/unembed, as an extension
2003 of the given ``field``.
2005 if field
in cls
._quaternion
_extension
:
2006 return cls
._quaternion
_extension
[field
]
2008 Q
= QuaternionAlgebra(field
,-1,-1)
2010 cls
._quaternion
_extension
[field
] = Q
2014 def dimension_over_reals():
2018 def real_embed(cls
,M
):
2020 Embed the n-by-n quaternion matrix ``M`` into the space of real
2021 matrices of size 4n-by-4n by first sending each quaternion entry `z
2022 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
2023 c+di],[-c + di, a-bi]]`, and then embedding those into a real
2028 sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA
2032 sage: Q = QuaternionAlgebra(QQ,-1,-1)
2033 sage: i,j,k = Q.gens()
2034 sage: x = 1 + 2*i + 3*j + 4*k
2035 sage: M = matrix(Q, 1, [[x]])
2036 sage: QuaternionMatrixEJA.real_embed(M)
2042 Embedding is a homomorphism (isomorphism, in fact)::
2044 sage: set_random_seed()
2045 sage: n = ZZ.random_element(2)
2046 sage: Q = QuaternionAlgebra(QQ,-1,-1)
2047 sage: X = random_matrix(Q, n)
2048 sage: Y = random_matrix(Q, n)
2049 sage: Xe = QuaternionMatrixEJA.real_embed(X)
2050 sage: Ye = QuaternionMatrixEJA.real_embed(Y)
2051 sage: XYe = QuaternionMatrixEJA.real_embed(X*Y)
2056 super(QuaternionMatrixEJA
,cls
).real_embed(M
)
2057 quaternions
= M
.base_ring()
2060 F
= QuadraticField(-1, 'I')
2065 t
= z
.coefficient_tuple()
2070 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
2071 [-c
+ d
*i
, a
- b
*i
]])
2072 realM
= ComplexMatrixEJA
.real_embed(cplxM
)
2073 blocks
.append(realM
)
2075 # We should have real entries by now, so use the realest field
2076 # we've got for the return value.
2077 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
2082 def real_unembed(cls
,M
):
2084 The inverse of _embed_quaternion_matrix().
2088 sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA
2092 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
2093 ....: [-2, 1, -4, 3],
2094 ....: [-3, 4, 1, -2],
2095 ....: [-4, -3, 2, 1]])
2096 sage: QuaternionMatrixEJA.real_unembed(M)
2097 [1 + 2*i + 3*j + 4*k]
2101 Unembedding is the inverse of embedding::
2103 sage: set_random_seed()
2104 sage: Q = QuaternionAlgebra(QQ, -1, -1)
2105 sage: M = random_matrix(Q, 3)
2106 sage: Me = QuaternionMatrixEJA.real_embed(M)
2107 sage: QuaternionMatrixEJA.real_unembed(Me) == M
2111 super(QuaternionMatrixEJA
,cls
).real_unembed(M
)
2113 d
= cls
.dimension_over_reals()
2115 # Use the base ring of the matrix to ensure that its entries can be
2116 # multiplied by elements of the quaternion algebra.
2117 Q
= cls
.quaternion_extension(M
.base_ring())
2120 # Go top-left to bottom-right (reading order), converting every
2121 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
2124 for l
in range(n
/d
):
2125 for m
in range(n
/d
):
2126 submat
= ComplexMatrixEJA
.real_unembed(
2127 M
[d
*l
:d
*l
+d
,d
*m
:d
*m
+d
] )
2128 if submat
[0,0] != submat
[1,1].conjugate():
2129 raise ValueError('bad on-diagonal submatrix')
2130 if submat
[0,1] != -submat
[1,0].conjugate():
2131 raise ValueError('bad off-diagonal submatrix')
2132 z
= submat
[0,0].real()
2133 z
+= submat
[0,0].imag()*i
2134 z
+= submat
[0,1].real()*j
2135 z
+= submat
[0,1].imag()*k
2138 return matrix(Q
, n
/d
, elements
)
2141 class QuaternionHermitianEJA(ConcreteEJA
, QuaternionMatrixEJA
):
2143 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
2144 matrices, the usual symmetric Jordan product, and the
2145 real-part-of-trace inner product. It has dimension `2n^2 - n` over
2150 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2154 In theory, our "field" can be any subfield of the reals::
2156 sage: QuaternionHermitianEJA(2, field=RDF)
2157 Euclidean Jordan algebra of dimension 6 over Real Double Field
2158 sage: QuaternionHermitianEJA(2, field=RR)
2159 Euclidean Jordan algebra of dimension 6 over Real Field with
2160 53 bits of precision
2164 The dimension of this algebra is `2*n^2 - n`::
2166 sage: set_random_seed()
2167 sage: n_max = QuaternionHermitianEJA._max_random_instance_size()
2168 sage: n = ZZ.random_element(1, n_max)
2169 sage: J = QuaternionHermitianEJA(n)
2170 sage: J.dimension() == 2*(n^2) - n
2173 The Jordan multiplication is what we think it is::
2175 sage: set_random_seed()
2176 sage: J = QuaternionHermitianEJA.random_instance()
2177 sage: x,y = J.random_elements(2)
2178 sage: actual = (x*y).to_matrix()
2179 sage: X = x.to_matrix()
2180 sage: Y = y.to_matrix()
2181 sage: expected = (X*Y + Y*X)/2
2182 sage: actual == expected
2184 sage: J(expected) == x*y
2187 We can change the generator prefix::
2189 sage: QuaternionHermitianEJA(2, prefix='a').gens()
2190 (a0, a1, a2, a3, a4, a5)
2192 We can construct the (trivial) algebra of rank zero::
2194 sage: QuaternionHermitianEJA(0)
2195 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2199 def _denormalized_basis(cls
, n
):
2201 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
2203 Why do we embed these? Basically, because all of numerical
2204 linear algebra assumes that you're working with vectors consisting
2205 of `n` entries from a field and scalars from the same field. There's
2206 no way to tell SageMath that (for example) the vectors contain
2207 complex numbers, while the scalar field is real.
2211 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2215 sage: set_random_seed()
2216 sage: n = ZZ.random_element(1,5)
2217 sage: B = QuaternionHermitianEJA._denormalized_basis(n)
2218 sage: all( M.is_symmetric() for M in B )
2223 Q
= QuaternionAlgebra(QQ
,-1,-1)
2226 # This is like the symmetric case, but we need to be careful:
2228 # * We want conjugate-symmetry, not just symmetry.
2229 # * The diagonal will (as a result) be real.
2232 Eij
= matrix
.zero(Q
,n
)
2234 for j
in range(i
+1):
2238 Sij
= cls
.real_embed(Eij
)
2241 # The second, third, and fourth ones have a minus
2242 # because they're conjugated.
2243 # Eij = Eij + Eij.transpose()
2245 Sij_real
= cls
.real_embed(Eij
)
2247 # Eij = I*(Eij - Eij.transpose())
2250 Sij_I
= cls
.real_embed(Eij
)
2252 # Eij = J*(Eij - Eij.transpose())
2255 Sij_J
= cls
.real_embed(Eij
)
2257 # Eij = K*(Eij - Eij.transpose())
2260 Sij_K
= cls
.real_embed(Eij
)
2266 # Since we embedded these, we can drop back to the "field" that we
2267 # started with instead of the quaternion algebra "Q".
2268 return tuple( s
.change_ring(field
) for s
in S
)
2271 def __init__(self
, n
, **kwargs
):
2272 # We know this is a valid EJA, but will double-check
2273 # if the user passes check_axioms=True.
2274 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2276 super(QuaternionHermitianEJA
, self
).__init
__(self
._denormalized
_basis
(n
),
2277 self
.jordan_product
,
2278 self
.trace_inner_product
,
2280 # TODO: this could be factored out somehow, but is left here
2281 # because the MatrixEJA is not presently a subclass of the
2282 # FDEJA class that defines rank() and one().
2283 self
.rank
.set_cache(n
)
2284 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
2285 self
.one
.set_cache(self(idV
))
2289 def _max_random_instance_size():
2291 The maximum rank of a random QuaternionHermitianEJA.
2293 return 2 # Dimension 6
2296 def random_instance(cls
, **kwargs
):
2298 Return a random instance of this type of algebra.
2300 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2301 return cls(n
, **kwargs
)
2304 class HadamardEJA(ConcreteEJA
):
2306 Return the Euclidean Jordan Algebra corresponding to the set
2307 `R^n` under the Hadamard product.
2309 Note: this is nothing more than the Cartesian product of ``n``
2310 copies of the spin algebra. Once Cartesian product algebras
2311 are implemented, this can go.
2315 sage: from mjo.eja.eja_algebra import HadamardEJA
2319 This multiplication table can be verified by hand::
2321 sage: J = HadamardEJA(3)
2322 sage: e0,e1,e2 = J.gens()
2338 We can change the generator prefix::
2340 sage: HadamardEJA(3, prefix='r').gens()
2344 def __init__(self
, n
, **kwargs
):
2346 jordan_product
= lambda x
,y
: x
2347 inner_product
= lambda x
,y
: x
2349 def jordan_product(x
,y
):
2351 return P( xi
*yi
for (xi
,yi
) in zip(x
,y
) )
2353 def inner_product(x
,y
):
2356 # New defaults for keyword arguments. Don't orthonormalize
2357 # because our basis is already orthonormal with respect to our
2358 # inner-product. Don't check the axioms, because we know this
2359 # is a valid EJA... but do double-check if the user passes
2360 # check_axioms=True. Note: we DON'T override the "check_field"
2361 # default here, because the user can pass in a field!
2362 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2363 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2365 column_basis
= tuple( b
.column() for b
in FreeModule(ZZ
, n
).basis() )
2366 super().__init
__(column_basis
, jordan_product
, inner_product
, **kwargs
)
2367 self
.rank
.set_cache(n
)
2370 self
.one
.set_cache( self
.zero() )
2372 self
.one
.set_cache( sum(self
.gens()) )
2375 def _max_random_instance_size():
2377 The maximum dimension of a random HadamardEJA.
2382 def random_instance(cls
, **kwargs
):
2384 Return a random instance of this type of algebra.
2386 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2387 return cls(n
, **kwargs
)
2390 class BilinearFormEJA(ConcreteEJA
):
2392 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2393 with the half-trace inner product and jordan product ``x*y =
2394 (<Bx,y>,y_bar>, x0*y_bar + y0*x_bar)`` where `B = 1 \times B22` is
2395 a symmetric positive-definite "bilinear form" matrix. Its
2396 dimension is the size of `B`, and it has rank two in dimensions
2397 larger than two. It reduces to the ``JordanSpinEJA`` when `B` is
2398 the identity matrix of order ``n``.
2400 We insist that the one-by-one upper-left identity block of `B` be
2401 passed in as well so that we can be passed a matrix of size zero
2402 to construct a trivial algebra.
2406 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
2407 ....: JordanSpinEJA)
2411 When no bilinear form is specified, the identity matrix is used,
2412 and the resulting algebra is the Jordan spin algebra::
2414 sage: B = matrix.identity(AA,3)
2415 sage: J0 = BilinearFormEJA(B)
2416 sage: J1 = JordanSpinEJA(3)
2417 sage: J0.multiplication_table() == J0.multiplication_table()
2420 An error is raised if the matrix `B` does not correspond to a
2421 positive-definite bilinear form::
2423 sage: B = matrix.random(QQ,2,3)
2424 sage: J = BilinearFormEJA(B)
2425 Traceback (most recent call last):
2427 ValueError: bilinear form is not positive-definite
2428 sage: B = matrix.zero(QQ,3)
2429 sage: J = BilinearFormEJA(B)
2430 Traceback (most recent call last):
2432 ValueError: bilinear form is not positive-definite
2436 We can create a zero-dimensional algebra::
2438 sage: B = matrix.identity(AA,0)
2439 sage: J = BilinearFormEJA(B)
2443 We can check the multiplication condition given in the Jordan, von
2444 Neumann, and Wigner paper (and also discussed on my "On the
2445 symmetry..." paper). Note that this relies heavily on the standard
2446 choice of basis, as does anything utilizing the bilinear form
2447 matrix. We opt not to orthonormalize the basis, because if we
2448 did, we would have to normalize the `s_{i}` in a similar manner::
2450 sage: set_random_seed()
2451 sage: n = ZZ.random_element(5)
2452 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2453 sage: B11 = matrix.identity(QQ,1)
2454 sage: B22 = M.transpose()*M
2455 sage: B = block_matrix(2,2,[ [B11,0 ],
2457 sage: J = BilinearFormEJA(B, orthonormalize=False)
2458 sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
2459 sage: V = J.vector_space()
2460 sage: sis = [ J( V([0] + (M.inverse()*ei).list()).column() )
2461 ....: for ei in eis ]
2462 sage: actual = [ sis[i]*sis[j]
2463 ....: for i in range(n-1)
2464 ....: for j in range(n-1) ]
2465 sage: expected = [ J.one() if i == j else J.zero()
2466 ....: for i in range(n-1)
2467 ....: for j in range(n-1) ]
2468 sage: actual == expected
2472 def __init__(self
, B
, **kwargs
):
2473 # The matrix "B" is supplied by the user in most cases,
2474 # so it makes sense to check whether or not its positive-
2475 # definite unless we are specifically asked not to...
2476 if ("check_axioms" not in kwargs
) or kwargs
["check_axioms"]:
2477 if not B
.is_positive_definite():
2478 raise ValueError("bilinear form is not positive-definite")
2480 # However, all of the other data for this EJA is computed
2481 # by us in manner that guarantees the axioms are
2482 # satisfied. So, again, unless we are specifically asked to
2483 # verify things, we'll skip the rest of the checks.
2484 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2486 def inner_product(x
,y
):
2487 return (y
.T
*B
*x
)[0,0]
2489 def jordan_product(x
,y
):
2495 z0
= inner_product(y
,x
)
2496 zbar
= y0
*xbar
+ x0
*ybar
2497 return P([z0
] + zbar
.list())
2500 column_basis
= tuple( b
.column() for b
in FreeModule(ZZ
, n
).basis() )
2501 super(BilinearFormEJA
, self
).__init
__(column_basis
,
2506 # The rank of this algebra is two, unless we're in a
2507 # one-dimensional ambient space (because the rank is bounded
2508 # by the ambient dimension).
2509 self
.rank
.set_cache(min(n
,2))
2512 self
.one
.set_cache( self
.zero() )
2514 self
.one
.set_cache( self
.monomial(0) )
2517 def _max_random_instance_size():
2519 The maximum dimension of a random BilinearFormEJA.
2524 def random_instance(cls
, **kwargs
):
2526 Return a random instance of this algebra.
2528 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2530 B
= matrix
.identity(ZZ
, n
)
2531 return cls(B
, **kwargs
)
2533 B11
= matrix
.identity(ZZ
, 1)
2534 M
= matrix
.random(ZZ
, n
-1)
2535 I
= matrix
.identity(ZZ
, n
-1)
2537 while alpha
.is_zero():
2538 alpha
= ZZ
.random_element().abs()
2539 B22
= M
.transpose()*M
+ alpha
*I
2541 from sage
.matrix
.special
import block_matrix
2542 B
= block_matrix(2,2, [ [B11
, ZZ(0) ],
2545 return cls(B
, **kwargs
)
2548 class JordanSpinEJA(BilinearFormEJA
):
2550 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2551 with the usual inner product and jordan product ``x*y =
2552 (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2557 sage: from mjo.eja.eja_algebra import JordanSpinEJA
2561 This multiplication table can be verified by hand::
2563 sage: J = JordanSpinEJA(4)
2564 sage: e0,e1,e2,e3 = J.gens()
2580 We can change the generator prefix::
2582 sage: JordanSpinEJA(2, prefix='B').gens()
2587 Ensure that we have the usual inner product on `R^n`::
2589 sage: set_random_seed()
2590 sage: J = JordanSpinEJA.random_instance()
2591 sage: x,y = J.random_elements(2)
2592 sage: actual = x.inner_product(y)
2593 sage: expected = x.to_vector().inner_product(y.to_vector())
2594 sage: actual == expected
2598 def __init__(self
, n
, **kwargs
):
2599 # This is a special case of the BilinearFormEJA with the
2600 # identity matrix as its bilinear form.
2601 B
= matrix
.identity(ZZ
, n
)
2603 # Don't orthonormalize because our basis is already
2604 # orthonormal with respect to our inner-product.
2605 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2607 # But also don't pass check_field=False here, because the user
2608 # can pass in a field!
2609 super(JordanSpinEJA
, self
).__init
__(B
, **kwargs
)
2612 def _max_random_instance_size():
2614 The maximum dimension of a random JordanSpinEJA.
2619 def random_instance(cls
, **kwargs
):
2621 Return a random instance of this type of algebra.
2623 Needed here to override the implementation for ``BilinearFormEJA``.
2625 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2626 return cls(n
, **kwargs
)
2629 class TrivialEJA(ConcreteEJA
):
2631 The trivial Euclidean Jordan algebra consisting of only a zero element.
2635 sage: from mjo.eja.eja_algebra import TrivialEJA
2639 sage: J = TrivialEJA()
2646 sage: 7*J.one()*12*J.one()
2648 sage: J.one().inner_product(J.one())
2650 sage: J.one().norm()
2652 sage: J.one().subalgebra_generated_by()
2653 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2658 def __init__(self
, **kwargs
):
2659 jordan_product
= lambda x
,y
: x
2660 inner_product
= lambda x
,y
: 0
2663 # New defaults for keyword arguments
2664 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2665 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2667 super(TrivialEJA
, self
).__init
__(basis
,
2671 # The rank is zero using my definition, namely the dimension of the
2672 # largest subalgebra generated by any element.
2673 self
.rank
.set_cache(0)
2674 self
.one
.set_cache( self
.zero() )
2677 def random_instance(cls
, **kwargs
):
2678 # We don't take a "size" argument so the superclass method is
2679 # inappropriate for us.
2680 return cls(**kwargs
)
2683 class CartesianProductEJA(CombinatorialFreeModule_CartesianProduct
,
2684 FiniteDimensionalEJA
):
2686 The external (orthogonal) direct sum of two or more Euclidean
2687 Jordan algebras. Every Euclidean Jordan algebra decomposes into an
2688 orthogonal direct sum of simple Euclidean Jordan algebras which is
2689 then isometric to a Cartesian product, so no generality is lost by
2690 providing only this construction.
2694 sage: from mjo.eja.eja_algebra import (random_eja,
2695 ....: CartesianProductEJA,
2697 ....: JordanSpinEJA,
2698 ....: RealSymmetricEJA)
2702 The Jordan product is inherited from our factors and implemented by
2703 our CombinatorialFreeModule Cartesian product superclass::
2705 sage: set_random_seed()
2706 sage: J1 = HadamardEJA(2)
2707 sage: J2 = RealSymmetricEJA(2)
2708 sage: J = cartesian_product([J1,J2])
2709 sage: x,y = J.random_elements(2)
2713 The ability to retrieve the original factors is implemented by our
2714 CombinatorialFreeModule Cartesian product superclass::
2716 sage: J1 = HadamardEJA(2, field=QQ)
2717 sage: J2 = JordanSpinEJA(3, field=QQ)
2718 sage: J = cartesian_product([J1,J2])
2719 sage: J.cartesian_factors()
2720 (Euclidean Jordan algebra of dimension 2 over Rational Field,
2721 Euclidean Jordan algebra of dimension 3 over Rational Field)
2723 You can provide more than two factors::
2725 sage: J1 = HadamardEJA(2)
2726 sage: J2 = JordanSpinEJA(3)
2727 sage: J3 = RealSymmetricEJA(3)
2728 sage: cartesian_product([J1,J2,J3])
2729 Euclidean Jordan algebra of dimension 2 over Algebraic Real
2730 Field (+) Euclidean Jordan algebra of dimension 3 over Algebraic
2731 Real Field (+) Euclidean Jordan algebra of dimension 6 over
2732 Algebraic Real Field
2734 Rank is additive on a Cartesian product::
2736 sage: J1 = HadamardEJA(1)
2737 sage: J2 = RealSymmetricEJA(2)
2738 sage: J = cartesian_product([J1,J2])
2739 sage: J1.rank.clear_cache()
2740 sage: J2.rank.clear_cache()
2741 sage: J.rank.clear_cache()
2744 sage: J.rank() == J1.rank() + J2.rank()
2747 The same rank computation works over the rationals, with whatever
2750 sage: J1 = HadamardEJA(1, field=QQ, orthonormalize=False)
2751 sage: J2 = RealSymmetricEJA(2, field=QQ, orthonormalize=False)
2752 sage: J = cartesian_product([J1,J2])
2753 sage: J1.rank.clear_cache()
2754 sage: J2.rank.clear_cache()
2755 sage: J.rank.clear_cache()
2758 sage: J.rank() == J1.rank() + J2.rank()
2763 All factors must share the same base field::
2765 sage: J1 = HadamardEJA(2, field=QQ)
2766 sage: J2 = RealSymmetricEJA(2)
2767 sage: CartesianProductEJA((J1,J2))
2768 Traceback (most recent call last):
2770 ValueError: all factors must share the same base field
2772 The "cached" Jordan and inner products are the componentwise
2775 sage: set_random_seed()
2776 sage: J1 = random_eja()
2777 sage: J2 = random_eja()
2778 sage: J = cartesian_product([J1,J2])
2779 sage: x,y = J.random_elements(2)
2780 sage: x*y == J.cartesian_jordan_product(x,y)
2782 sage: x.inner_product(y) == J.cartesian_inner_product(x,y)
2785 The cached unit element is the same one that would be computed::
2787 sage: set_random_seed() # long time
2788 sage: J1 = random_eja() # long time
2789 sage: J2 = random_eja() # long time
2790 sage: J = cartesian_product([J1,J2]) # long time
2791 sage: actual = J.one() # long time
2792 sage: J.one.clear_cache() # long time
2793 sage: expected = J.one() # long time
2794 sage: actual == expected # long time
2798 def __init__(self
, modules
, **kwargs
):
2799 CombinatorialFreeModule_CartesianProduct
.__init
__(self
,
2802 field
= modules
[0].base_ring()
2803 if not all( J
.base_ring() == field
for J
in modules
):
2804 raise ValueError("all factors must share the same base field")
2806 basis
= tuple( b
.to_vector().column() for b
in self
.basis() )
2808 # Define jordan/inner products that operate on the basis.
2809 def jordan_product(x_mat
,y_mat
):
2810 x
= self
.from_vector(_mat2vec(x_mat
))
2811 y
= self
.from_vector(_mat2vec(y_mat
))
2812 return self
.cartesian_jordan_product(x
,y
).to_vector().column()
2814 def inner_product(x_mat
, y_mat
):
2815 x
= self
.from_vector(_mat2vec(x_mat
))
2816 y
= self
.from_vector(_mat2vec(y_mat
))
2817 return self
.cartesian_inner_product(x
,y
)
2819 # Use whatever category the superclass came up with. Usually
2820 # some join of the EJA and Cartesian product
2821 # categories. There's no need to check the field since it
2822 # already came from an EJA.
2824 # If you want the basis to be orthonormalized, orthonormalize
2826 FiniteDimensionalEJA
.__init
__(self
,
2831 orthonormalize
=False,
2834 category
=self
.category())
2836 ones
= tuple(J
.one() for J
in modules
)
2837 self
.one
.set_cache(self
._cartesian
_product
_of
_elements
(ones
))
2838 self
.rank
.set_cache(sum(J
.rank() for J
in modules
))
2840 # Now that everything else is ready, we clobber our computed
2841 # matrix basis with the "correct" one consisting of ordered
2842 # tuples. Since we didn't orthonormalize our basis, we can
2843 # create these from the basis that was handed to us; that is,
2844 # we don't need to use the one that the earlier __init__()
2845 # method came up with.
2846 m
= len(self
.cartesian_factors())
2847 MS
= self
.matrix_space()
2848 self
._matrix
_basis
= tuple(
2849 MS(tuple( self
.cartesian_projection(i
)(b
).to_matrix()
2850 for i
in range(m
) ))
2851 for b
in self
.basis()
2854 def matrix_space(self
):
2856 Return the space that our matrix basis lives in as a Cartesian
2861 sage: from mjo.eja.eja_algebra import (HadamardEJA,
2862 ....: RealSymmetricEJA)
2866 sage: J1 = HadamardEJA(1)
2867 sage: J2 = RealSymmetricEJA(2)
2868 sage: J = cartesian_product([J1,J2])
2869 sage: J.matrix_space()
2870 The Cartesian product of (Full MatrixSpace of 1 by 1 dense
2871 matrices over Algebraic Real Field, Full MatrixSpace of 2
2872 by 2 dense matrices over Algebraic Real Field)
2875 from sage
.categories
.cartesian_product
import cartesian_product
2876 return cartesian_product( [J
.matrix_space()
2877 for J
in self
.cartesian_factors()] )
2880 def cartesian_projection(self
, i
):
2884 sage: from mjo.eja.eja_algebra import (random_eja,
2885 ....: JordanSpinEJA,
2887 ....: RealSymmetricEJA,
2888 ....: ComplexHermitianEJA)
2892 The projection morphisms are Euclidean Jordan algebra
2895 sage: J1 = HadamardEJA(2)
2896 sage: J2 = RealSymmetricEJA(2)
2897 sage: J = cartesian_product([J1,J2])
2898 sage: J.cartesian_projection(0)
2899 Linear operator between finite-dimensional Euclidean Jordan
2900 algebras represented by the matrix:
2903 Domain: Euclidean Jordan algebra of dimension 2 over Algebraic
2904 Real Field (+) Euclidean Jordan algebra of dimension 3 over
2905 Algebraic Real Field
2906 Codomain: Euclidean Jordan algebra of dimension 2 over Algebraic
2908 sage: J.cartesian_projection(1)
2909 Linear operator between finite-dimensional Euclidean Jordan
2910 algebras represented by the matrix:
2914 Domain: Euclidean Jordan algebra of dimension 2 over Algebraic
2915 Real Field (+) Euclidean Jordan algebra of dimension 3 over
2916 Algebraic Real Field
2917 Codomain: Euclidean Jordan algebra of dimension 3 over Algebraic
2920 The projections work the way you'd expect on the vector
2921 representation of an element::
2923 sage: J1 = JordanSpinEJA(2)
2924 sage: J2 = ComplexHermitianEJA(2)
2925 sage: J = cartesian_product([J1,J2])
2926 sage: pi_left = J.cartesian_projection(0)
2927 sage: pi_right = J.cartesian_projection(1)
2928 sage: pi_left(J.one()).to_vector()
2930 sage: pi_right(J.one()).to_vector()
2932 sage: J.one().to_vector()
2937 The answer never changes::
2939 sage: set_random_seed()
2940 sage: J1 = random_eja()
2941 sage: J2 = random_eja()
2942 sage: J = cartesian_product([J1,J2])
2943 sage: P0 = J.cartesian_projection(0)
2944 sage: P1 = J.cartesian_projection(0)
2949 Ji
= self
.cartesian_factors()[i
]
2950 # Requires the fix on Trac 31421/31422 to work!
2951 Pi
= super().cartesian_projection(i
)
2952 return FiniteDimensionalEJAOperator(self
,Ji
,Pi
.matrix())
2955 def cartesian_embedding(self
, i
):
2959 sage: from mjo.eja.eja_algebra import (random_eja,
2960 ....: JordanSpinEJA,
2962 ....: RealSymmetricEJA)
2966 The embedding morphisms are Euclidean Jordan algebra
2969 sage: J1 = HadamardEJA(2)
2970 sage: J2 = RealSymmetricEJA(2)
2971 sage: J = cartesian_product([J1,J2])
2972 sage: J.cartesian_embedding(0)
2973 Linear operator between finite-dimensional Euclidean Jordan
2974 algebras represented by the matrix:
2980 Domain: Euclidean Jordan algebra of dimension 2 over
2981 Algebraic Real Field
2982 Codomain: Euclidean Jordan algebra of dimension 2 over
2983 Algebraic Real Field (+) Euclidean Jordan algebra of
2984 dimension 3 over Algebraic Real Field
2985 sage: J.cartesian_embedding(1)
2986 Linear operator between finite-dimensional Euclidean Jordan
2987 algebras represented by the matrix:
2993 Domain: Euclidean Jordan algebra of dimension 3 over
2994 Algebraic Real Field
2995 Codomain: Euclidean Jordan algebra of dimension 2 over
2996 Algebraic Real Field (+) Euclidean Jordan algebra of
2997 dimension 3 over Algebraic Real Field
2999 The embeddings work the way you'd expect on the vector
3000 representation of an element::
3002 sage: J1 = JordanSpinEJA(3)
3003 sage: J2 = RealSymmetricEJA(2)
3004 sage: J = cartesian_product([J1,J2])
3005 sage: iota_left = J.cartesian_embedding(0)
3006 sage: iota_right = J.cartesian_embedding(1)
3007 sage: iota_left(J1.zero()) == J.zero()
3009 sage: iota_right(J2.zero()) == J.zero()
3011 sage: J1.one().to_vector()
3013 sage: iota_left(J1.one()).to_vector()
3015 sage: J2.one().to_vector()
3017 sage: iota_right(J2.one()).to_vector()
3019 sage: J.one().to_vector()
3024 The answer never changes::
3026 sage: set_random_seed()
3027 sage: J1 = random_eja()
3028 sage: J2 = random_eja()
3029 sage: J = cartesian_product([J1,J2])
3030 sage: E0 = J.cartesian_embedding(0)
3031 sage: E1 = J.cartesian_embedding(0)
3035 Composing a projection with the corresponding inclusion should
3036 produce the identity map, and mismatching them should produce
3039 sage: set_random_seed()
3040 sage: J1 = random_eja()
3041 sage: J2 = random_eja()
3042 sage: J = cartesian_product([J1,J2])
3043 sage: iota_left = J.cartesian_embedding(0)
3044 sage: iota_right = J.cartesian_embedding(1)
3045 sage: pi_left = J.cartesian_projection(0)
3046 sage: pi_right = J.cartesian_projection(1)
3047 sage: pi_left*iota_left == J1.one().operator()
3049 sage: pi_right*iota_right == J2.one().operator()
3051 sage: (pi_left*iota_right).is_zero()
3053 sage: (pi_right*iota_left).is_zero()
3057 Ji
= self
.cartesian_factors()[i
]
3058 # Requires the fix on Trac 31421/31422 to work!
3059 Ei
= super().cartesian_embedding(i
)
3060 return FiniteDimensionalEJAOperator(Ji
,self
,Ei
.matrix())
3063 def cartesian_jordan_product(self
, x
, y
):
3065 The componentwise Jordan product.
3067 We project ``x`` and ``y`` onto our factors, and add up the
3068 Jordan products from the subalgebras. This may still be useful
3069 after (if) the default Jordan product in the Cartesian product
3070 algebra is overridden.
3074 sage: from mjo.eja.eja_algebra import (HadamardEJA,
3075 ....: JordanSpinEJA)
3079 sage: J1 = HadamardEJA(3)
3080 sage: J2 = JordanSpinEJA(3)
3081 sage: J = cartesian_product([J1,J2])
3082 sage: x1 = J1.from_vector(vector(QQ,(1,2,1)))
3083 sage: y1 = J1.from_vector(vector(QQ,(1,0,2)))
3084 sage: x2 = J2.from_vector(vector(QQ,(1,2,3)))
3085 sage: y2 = J2.from_vector(vector(QQ,(1,1,1)))
3086 sage: z1 = J.from_vector(vector(QQ,(1,2,1,1,2,3)))
3087 sage: z2 = J.from_vector(vector(QQ,(1,0,2,1,1,1)))
3088 sage: (x1*y1).to_vector()
3090 sage: (x2*y2).to_vector()
3092 sage: J.cartesian_jordan_product(z1,z2).to_vector()
3096 m
= len(self
.cartesian_factors())
3097 projections
= ( self
.cartesian_projection(i
) for i
in range(m
) )
3098 products
= ( P(x
)*P(y
) for P
in projections
)
3099 return self
._cartesian
_product
_of
_elements
(tuple(products
))
3101 def cartesian_inner_product(self
, x
, y
):
3103 The standard componentwise Cartesian inner-product.
3105 We project ``x`` and ``y`` onto our factors, and add up the
3106 inner-products from the subalgebras. This may still be useful
3107 after (if) the default inner product in the Cartesian product
3108 algebra is overridden.
3112 sage: from mjo.eja.eja_algebra import (HadamardEJA,
3113 ....: QuaternionHermitianEJA)
3117 sage: J1 = HadamardEJA(3,field=QQ)
3118 sage: J2 = QuaternionHermitianEJA(2,field=QQ,orthonormalize=False)
3119 sage: J = cartesian_product([J1,J2])
3124 sage: x1.inner_product(x2)
3126 sage: y1.inner_product(y2)
3128 sage: z1 = J._cartesian_product_of_elements((x1,y1))
3129 sage: z2 = J._cartesian_product_of_elements((x2,y2))
3130 sage: J.cartesian_inner_product(z1,z2)
3134 m
= len(self
.cartesian_factors())
3135 projections
= ( self
.cartesian_projection(i
) for i
in range(m
) )
3136 return sum( P(x
).inner_product(P(y
)) for P
in projections
)
3139 def _element_constructor_(self
, elt
):
3141 Construct an element of this algebra from an ordered tuple.
3143 We just apply the element constructor from each of our factors
3144 to the corresponding component of the tuple, and package up
3149 sage: from mjo.eja.eja_algebra import (HadamardEJA,
3150 ....: RealSymmetricEJA)
3154 sage: J1 = HadamardEJA(3)
3155 sage: J2 = RealSymmetricEJA(2)
3156 sage: J = cartesian_product([J1,J2])
3157 sage: J( (J1.matrix_basis()[1], J2.matrix_basis()[2]) )
3160 m
= len(self
.cartesian_factors())
3162 z
= tuple( self
.cartesian_factors()[i
](elt
[i
]) for i
in range(m
) )
3163 return self
._cartesian
_product
_of
_elements
(z
)
3165 raise ValueError("not an element of this algebra")
3167 Element
= CartesianProductEJAElement
3170 FiniteDimensionalEJA
.CartesianProduct
= CartesianProductEJA
3171 random_eja
= ConcreteEJA
.random_instance