2 Euclidean Jordan Algebras. These are formally-real Jordan Algebras;
3 specifically those where u^2 + v^2 = 0 implies that u = v = 0. They
4 are used in optimization, and have some additional nice methods beyond
5 what can be supported in a general Jordan Algebra.
10 sage: from mjo.eja.eja_algebra import random_eja
15 Euclidean Jordan algebra of dimension...
19 from itertools
import repeat
21 from sage
.algebras
.quatalg
.quaternion_algebra
import QuaternionAlgebra
22 from sage
.categories
.magmatic_algebras
import MagmaticAlgebras
23 from sage
.categories
.sets_cat
import cartesian_product
24 from sage
.combinat
.free_module
import (CombinatorialFreeModule
,
25 CombinatorialFreeModule_CartesianProduct
)
26 from sage
.matrix
.constructor
import matrix
27 from sage
.matrix
.matrix_space
import MatrixSpace
28 from sage
.misc
.cachefunc
import cached_method
29 from sage
.misc
.table
import table
30 from sage
.modules
.free_module
import FreeModule
, VectorSpace
31 from sage
.rings
.all
import (ZZ
, QQ
, AA
, QQbar
, RR
, RLF
, CLF
,
34 from mjo
.eja
.eja_element
import FiniteDimensionalEJAElement
35 from mjo
.eja
.eja_operator
import FiniteDimensionalEJAOperator
36 from mjo
.eja
.eja_utils
import _all2list
, _mat2vec
38 class FiniteDimensionalEJA(CombinatorialFreeModule
):
40 A finite-dimensional Euclidean Jordan algebra.
44 - basis -- a tuple of basis elements in "matrix form," which
45 must be the same form as the arguments to ``jordan_product``
46 and ``inner_product``. In reality, "matrix form" can be either
47 vectors, matrices, or a Cartesian product (ordered tuple)
48 of vectors or matrices. All of these would ideally be vector
49 spaces in sage with no special-casing needed; but in reality
50 we turn vectors into column-matrices and Cartesian products
51 `(a,b)` into column matrices `(a,b)^{T}` after converting
52 `a` and `b` themselves.
54 - jordan_product -- function of two elements (in matrix form)
55 that returns their jordan product in this algebra; this will
56 be applied to ``basis`` to compute a multiplication table for
59 - inner_product -- function of two elements (in matrix form) that
60 returns their inner product. This will be applied to ``basis`` to
61 compute an inner-product table (basically a matrix) for this algebra.
64 Element
= FiniteDimensionalEJAElement
73 cartesian_product
=False,
78 # Keep track of whether or not the matrix basis consists of
79 # tuples, since we need special cases for them damned near
80 # everywhere. This is INDEPENDENT of whether or not the
81 # algebra is a cartesian product, since a subalgebra of a
82 # cartesian product will have a basis of tuples, but will not
83 # in general itself be a cartesian product algebra.
84 self
._matrix
_basis
_is
_cartesian
= False
87 if hasattr(basis
[0], 'cartesian_factors'):
88 self
._matrix
_basis
_is
_cartesian
= True
91 if not field
.is_subring(RR
):
92 # Note: this does return true for the real algebraic
93 # field, the rationals, and any quadratic field where
94 # we've specified a real embedding.
95 raise ValueError("scalar field is not real")
97 # If the basis given to us wasn't over the field that it's
98 # supposed to be over, fix that. Or, you know, crash.
99 if not cartesian_product
:
100 # The field for a cartesian product algebra comes from one
101 # of its factors and is the same for all factors, so
102 # there's no need to "reapply" it on product algebras.
103 if self
._matrix
_basis
_is
_cartesian
:
104 # OK since if n == 0, the basis does not consist of tuples.
105 P
= basis
[0].parent()
106 basis
= tuple( P(tuple(b_i
.change_ring(field
) for b_i
in b
))
109 basis
= tuple( b
.change_ring(field
) for b
in basis
)
113 # Check commutativity of the Jordan and inner-products.
114 # This has to be done before we build the multiplication
115 # and inner-product tables/matrices, because we take
116 # advantage of symmetry in the process.
117 if not all( jordan_product(bi
,bj
) == jordan_product(bj
,bi
)
120 raise ValueError("Jordan product is not commutative")
122 if not all( inner_product(bi
,bj
) == inner_product(bj
,bi
)
125 raise ValueError("inner-product is not commutative")
128 category
= MagmaticAlgebras(field
).FiniteDimensional()
129 category
= category
.WithBasis().Unital()
131 # Element subalgebras can take advantage of this.
132 category
= category
.Associative()
133 if cartesian_product
:
134 category
= category
.CartesianProducts()
136 # Call the superclass constructor so that we can use its from_vector()
137 # method to build our multiplication table.
138 CombinatorialFreeModule
.__init
__(self
,
145 # Now comes all of the hard work. We'll be constructing an
146 # ambient vector space V that our (vectorized) basis lives in,
147 # as well as a subspace W of V spanned by those (vectorized)
148 # basis elements. The W-coordinates are the coefficients that
149 # we see in things like x = 1*e1 + 2*e2.
154 degree
= len(_all2list(basis
[0]))
156 # Build an ambient space that fits our matrix basis when
157 # written out as "long vectors."
158 V
= VectorSpace(field
, degree
)
160 # The matrix that will hole the orthonormal -> unorthonormal
161 # coordinate transformation.
162 self
._deortho
_matrix
= None
165 # Save a copy of the un-orthonormalized basis for later.
166 # Convert it to ambient V (vector) coordinates while we're
167 # at it, because we'd have to do it later anyway.
168 deortho_vector_basis
= tuple( V(_all2list(b
)) for b
in basis
)
170 from mjo
.eja
.eja_utils
import gram_schmidt
171 basis
= tuple(gram_schmidt(basis
, inner_product
))
173 # Save the (possibly orthonormalized) matrix basis for
175 self
._matrix
_basis
= basis
177 # Now create the vector space for the algebra, which will have
178 # its own set of non-ambient coordinates (in terms of the
180 vector_basis
= tuple( V(_all2list(b
)) for b
in basis
)
181 W
= V
.span_of_basis( vector_basis
, check
=check_axioms
)
184 # Now "W" is the vector space of our algebra coordinates. The
185 # variables "X1", "X2",... refer to the entries of vectors in
186 # W. Thus to convert back and forth between the orthonormal
187 # coordinates and the given ones, we need to stick the original
189 U
= V
.span_of_basis( deortho_vector_basis
, check
=check_axioms
)
190 self
._deortho
_matrix
= matrix( U
.coordinate_vector(q
)
191 for q
in vector_basis
)
194 # Now we actually compute the multiplication and inner-product
195 # tables/matrices using the possibly-orthonormalized basis.
196 self
._inner
_product
_matrix
= matrix
.identity(field
, n
)
197 self
._multiplication
_table
= [ [0 for j
in range(i
+1)]
200 # Note: the Jordan and inner-products are defined in terms
201 # of the ambient basis. It's important that their arguments
202 # are in ambient coordinates as well.
205 # ortho basis w.r.t. ambient coords
209 # The jordan product returns a matrixy answer, so we
210 # have to convert it to the algebra coordinates.
211 elt
= jordan_product(q_i
, q_j
)
212 elt
= W
.coordinate_vector(V(_all2list(elt
)))
213 self
._multiplication
_table
[i
][j
] = self
.from_vector(elt
)
215 if not orthonormalize
:
216 # If we're orthonormalizing the basis with respect
217 # to an inner-product, then the inner-product
218 # matrix with respect to the resulting basis is
219 # just going to be the identity.
220 ip
= inner_product(q_i
, q_j
)
221 self
._inner
_product
_matrix
[i
,j
] = ip
222 self
._inner
_product
_matrix
[j
,i
] = ip
224 self
._inner
_product
_matrix
._cache
= {'hermitian': True}
225 self
._inner
_product
_matrix
.set_immutable()
228 if not self
._is
_jordanian
():
229 raise ValueError("Jordan identity does not hold")
230 if not self
._inner
_product
_is
_associative
():
231 raise ValueError("inner product is not associative")
234 def _coerce_map_from_base_ring(self
):
236 Disable the map from the base ring into the algebra.
238 Performing a nonsense conversion like this automatically
239 is counterpedagogical. The fallback is to try the usual
240 element constructor, which should also fail.
244 sage: from mjo.eja.eja_algebra import random_eja
248 sage: set_random_seed()
249 sage: J = random_eja()
251 Traceback (most recent call last):
253 ValueError: not an element of this algebra
259 def product_on_basis(self
, i
, j
):
260 # We only stored the lower-triangular portion of the
261 # multiplication table.
263 return self
._multiplication
_table
[i
][j
]
265 return self
._multiplication
_table
[j
][i
]
267 def inner_product(self
, x
, y
):
269 The inner product associated with this Euclidean Jordan algebra.
271 Defaults to the trace inner product, but can be overridden by
272 subclasses if they are sure that the necessary properties are
277 sage: from mjo.eja.eja_algebra import (random_eja,
279 ....: BilinearFormEJA)
283 Our inner product is "associative," which means the following for
284 a symmetric bilinear form::
286 sage: set_random_seed()
287 sage: J = random_eja()
288 sage: x,y,z = J.random_elements(3)
289 sage: (x*y).inner_product(z) == y.inner_product(x*z)
294 Ensure that this is the usual inner product for the algebras
297 sage: set_random_seed()
298 sage: J = HadamardEJA.random_instance()
299 sage: x,y = J.random_elements(2)
300 sage: actual = x.inner_product(y)
301 sage: expected = x.to_vector().inner_product(y.to_vector())
302 sage: actual == expected
305 Ensure that this is one-half of the trace inner-product in a
306 BilinearFormEJA that isn't just the reals (when ``n`` isn't
307 one). This is in Faraut and Koranyi, and also my "On the
310 sage: set_random_seed()
311 sage: J = BilinearFormEJA.random_instance()
312 sage: n = J.dimension()
313 sage: x = J.random_element()
314 sage: y = J.random_element()
315 sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
319 B
= self
._inner
_product
_matrix
320 return (B
*x
.to_vector()).inner_product(y
.to_vector())
323 def is_associative(self
):
325 Return whether or not this algebra's Jordan product is associative.
329 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
333 sage: J = ComplexHermitianEJA(3, field=QQ, orthonormalize=False)
334 sage: J.is_associative()
336 sage: x = sum(J.gens())
337 sage: A = x.subalgebra_generated_by(orthonormalize=False)
338 sage: A.is_associative()
342 return "Associative" in self
.category().axioms()
344 def _is_jordanian(self
):
346 Whether or not this algebra's multiplication table respects the
347 Jordan identity `(x^{2})(xy) = x(x^{2}y)`.
349 We only check one arrangement of `x` and `y`, so for a
350 ``True`` result to be truly true, you should also check
351 :meth:`is_commutative`. This method should of course always
352 return ``True``, unless this algebra was constructed with
353 ``check_axioms=False`` and passed an invalid multiplication table.
355 return all( (self
.gens()[i
]**2)*(self
.gens()[i
]*self
.gens()[j
])
357 (self
.gens()[i
])*((self
.gens()[i
]**2)*self
.gens()[j
])
358 for i
in range(self
.dimension())
359 for j
in range(self
.dimension()) )
361 def _inner_product_is_associative(self
):
363 Return whether or not this algebra's inner product `B` is
364 associative; that is, whether or not `B(xy,z) = B(x,yz)`.
366 This method should of course always return ``True``, unless
367 this algebra was constructed with ``check_axioms=False`` and
368 passed an invalid Jordan or inner-product.
371 # Used to check whether or not something is zero in an inexact
372 # ring. This number is sufficient to allow the construction of
373 # QuaternionHermitianEJA(2, field=RDF) with check_axioms=True.
376 for i
in range(self
.dimension()):
377 for j
in range(self
.dimension()):
378 for k
in range(self
.dimension()):
382 diff
= (x
*y
).inner_product(z
) - x
.inner_product(y
*z
)
384 if self
.base_ring().is_exact():
388 if diff
.abs() > epsilon
:
393 def _element_constructor_(self
, elt
):
395 Construct an element of this algebra from its vector or matrix
398 This gets called only after the parent element _call_ method
399 fails to find a coercion for the argument.
403 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
405 ....: RealSymmetricEJA)
409 The identity in `S^n` is converted to the identity in the EJA::
411 sage: J = RealSymmetricEJA(3)
412 sage: I = matrix.identity(QQ,3)
413 sage: J(I) == J.one()
416 This skew-symmetric matrix can't be represented in the EJA::
418 sage: J = RealSymmetricEJA(3)
419 sage: A = matrix(QQ,3, lambda i,j: i-j)
421 Traceback (most recent call last):
423 ValueError: not an element of this algebra
425 Tuples work as well, provided that the matrix basis for the
426 algebra consists of them::
428 sage: J1 = HadamardEJA(3)
429 sage: J2 = RealSymmetricEJA(2)
430 sage: J = cartesian_product([J1,J2])
431 sage: J( (J1.matrix_basis()[1], J2.matrix_basis()[2]) )
436 Ensure that we can convert any element of the two non-matrix
437 simple algebras (whose matrix representations are columns)
438 back and forth faithfully::
440 sage: set_random_seed()
441 sage: J = HadamardEJA.random_instance()
442 sage: x = J.random_element()
443 sage: J(x.to_vector().column()) == x
445 sage: J = JordanSpinEJA.random_instance()
446 sage: x = J.random_element()
447 sage: J(x.to_vector().column()) == x
450 We cannot coerce elements between algebras just because their
451 matrix representations are compatible::
453 sage: J1 = HadamardEJA(3)
454 sage: J2 = JordanSpinEJA(3)
456 Traceback (most recent call last):
458 ValueError: not an element of this algebra
460 Traceback (most recent call last):
462 ValueError: not an element of this algebra
465 msg
= "not an element of this algebra"
466 if elt
in self
.base_ring():
467 # Ensure that no base ring -> algebra coercion is performed
468 # by this method. There's some stupidity in sage that would
469 # otherwise propagate to this method; for example, sage thinks
470 # that the integer 3 belongs to the space of 2-by-2 matrices.
471 raise ValueError(msg
)
474 # Try to convert a vector into a column-matrix...
476 except (AttributeError, TypeError):
477 # and ignore failure, because we weren't really expecting
478 # a vector as an argument anyway.
481 if elt
not in self
.matrix_space():
482 raise ValueError(msg
)
484 # Thanks for nothing! Matrix spaces aren't vector spaces in
485 # Sage, so we have to figure out its matrix-basis coordinates
486 # ourselves. We use the basis space's ring instead of the
487 # element's ring because the basis space might be an algebraic
488 # closure whereas the base ring of the 3-by-3 identity matrix
489 # could be QQ instead of QQbar.
491 # And, we also have to handle Cartesian product bases (when
492 # the matrix basis consists of tuples) here. The "good news"
493 # is that we're already converting everything to long vectors,
494 # and that strategy works for tuples as well.
496 # We pass check=False because the matrix basis is "guaranteed"
497 # to be linearly independent... right? Ha ha.
499 V
= VectorSpace(self
.base_ring(), len(elt
))
500 W
= V
.span_of_basis( (V(_all2list(s
)) for s
in self
.matrix_basis()),
504 coords
= W
.coordinate_vector(V(elt
))
505 except ArithmeticError: # vector is not in free module
506 raise ValueError(msg
)
508 return self
.from_vector(coords
)
512 Return a string representation of ``self``.
516 sage: from mjo.eja.eja_algebra import JordanSpinEJA
520 Ensure that it says what we think it says::
522 sage: JordanSpinEJA(2, field=AA)
523 Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
524 sage: JordanSpinEJA(3, field=RDF)
525 Euclidean Jordan algebra of dimension 3 over Real Double Field
528 fmt
= "Euclidean Jordan algebra of dimension {} over {}"
529 return fmt
.format(self
.dimension(), self
.base_ring())
533 def characteristic_polynomial_of(self
):
535 Return the algebra's "characteristic polynomial of" function,
536 which is itself a multivariate polynomial that, when evaluated
537 at the coordinates of some algebra element, returns that
538 element's characteristic polynomial.
540 The resulting polynomial has `n+1` variables, where `n` is the
541 dimension of this algebra. The first `n` variables correspond to
542 the coordinates of an algebra element: when evaluated at the
543 coordinates of an algebra element with respect to a certain
544 basis, the result is a univariate polynomial (in the one
545 remaining variable ``t``), namely the characteristic polynomial
550 sage: from mjo.eja.eja_algebra import JordanSpinEJA, TrivialEJA
554 The characteristic polynomial in the spin algebra is given in
555 Alizadeh, Example 11.11::
557 sage: J = JordanSpinEJA(3)
558 sage: p = J.characteristic_polynomial_of(); p
559 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
560 sage: xvec = J.one().to_vector()
564 By definition, the characteristic polynomial is a monic
565 degree-zero polynomial in a rank-zero algebra. Note that
566 Cayley-Hamilton is indeed satisfied since the polynomial
567 ``1`` evaluates to the identity element of the algebra on
570 sage: J = TrivialEJA()
571 sage: J.characteristic_polynomial_of()
578 # The list of coefficient polynomials a_0, a_1, a_2, ..., a_(r-1).
579 a
= self
._charpoly
_coefficients
()
581 # We go to a bit of trouble here to reorder the
582 # indeterminates, so that it's easier to evaluate the
583 # characteristic polynomial at x's coordinates and get back
584 # something in terms of t, which is what we want.
585 S
= PolynomialRing(self
.base_ring(),'t')
589 S
= PolynomialRing(S
, R
.variable_names())
592 return (t
**r
+ sum( a
[k
]*(t
**k
) for k
in range(r
) ))
594 def coordinate_polynomial_ring(self
):
596 The multivariate polynomial ring in which this algebra's
597 :meth:`characteristic_polynomial_of` lives.
601 sage: from mjo.eja.eja_algebra import (HadamardEJA,
602 ....: RealSymmetricEJA)
606 sage: J = HadamardEJA(2)
607 sage: J.coordinate_polynomial_ring()
608 Multivariate Polynomial Ring in X1, X2...
609 sage: J = RealSymmetricEJA(3,field=QQ,orthonormalize=False)
610 sage: J.coordinate_polynomial_ring()
611 Multivariate Polynomial Ring in X1, X2, X3, X4, X5, X6...
614 var_names
= tuple( "X%d" % z
for z
in range(1, self
.dimension()+1) )
615 return PolynomialRing(self
.base_ring(), var_names
)
617 def inner_product(self
, x
, y
):
619 The inner product associated with this Euclidean Jordan algebra.
621 Defaults to the trace inner product, but can be overridden by
622 subclasses if they are sure that the necessary properties are
627 sage: from mjo.eja.eja_algebra import (random_eja,
629 ....: BilinearFormEJA)
633 Our inner product is "associative," which means the following for
634 a symmetric bilinear form::
636 sage: set_random_seed()
637 sage: J = random_eja()
638 sage: x,y,z = J.random_elements(3)
639 sage: (x*y).inner_product(z) == y.inner_product(x*z)
644 Ensure that this is the usual inner product for the algebras
647 sage: set_random_seed()
648 sage: J = HadamardEJA.random_instance()
649 sage: x,y = J.random_elements(2)
650 sage: actual = x.inner_product(y)
651 sage: expected = x.to_vector().inner_product(y.to_vector())
652 sage: actual == expected
655 Ensure that this is one-half of the trace inner-product in a
656 BilinearFormEJA that isn't just the reals (when ``n`` isn't
657 one). This is in Faraut and Koranyi, and also my "On the
660 sage: set_random_seed()
661 sage: J = BilinearFormEJA.random_instance()
662 sage: n = J.dimension()
663 sage: x = J.random_element()
664 sage: y = J.random_element()
665 sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
668 B
= self
._inner
_product
_matrix
669 return (B
*x
.to_vector()).inner_product(y
.to_vector())
672 def is_trivial(self
):
674 Return whether or not this algebra is trivial.
676 A trivial algebra contains only the zero element.
680 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
685 sage: J = ComplexHermitianEJA(3)
691 sage: J = TrivialEJA()
696 return self
.dimension() == 0
699 def multiplication_table(self
):
701 Return a visual representation of this algebra's multiplication
702 table (on basis elements).
706 sage: from mjo.eja.eja_algebra import JordanSpinEJA
710 sage: J = JordanSpinEJA(4)
711 sage: J.multiplication_table()
712 +----++----+----+----+----+
713 | * || e0 | e1 | e2 | e3 |
714 +====++====+====+====+====+
715 | e0 || e0 | e1 | e2 | e3 |
716 +----++----+----+----+----+
717 | e1 || e1 | e0 | 0 | 0 |
718 +----++----+----+----+----+
719 | e2 || e2 | 0 | e0 | 0 |
720 +----++----+----+----+----+
721 | e3 || e3 | 0 | 0 | e0 |
722 +----++----+----+----+----+
726 # Prepend the header row.
727 M
= [["*"] + list(self
.gens())]
729 # And to each subsequent row, prepend an entry that belongs to
730 # the left-side "header column."
731 M
+= [ [self
.gens()[i
]] + [ self
.product_on_basis(i
,j
)
735 return table(M
, header_row
=True, header_column
=True, frame
=True)
738 def matrix_basis(self
):
740 Return an (often more natural) representation of this algebras
741 basis as an ordered tuple of matrices.
743 Every finite-dimensional Euclidean Jordan Algebra is a, up to
744 Jordan isomorphism, a direct sum of five simple
745 algebras---four of which comprise Hermitian matrices. And the
746 last type of algebra can of course be thought of as `n`-by-`1`
747 column matrices (ambiguusly called column vectors) to avoid
748 special cases. As a result, matrices (and column vectors) are
749 a natural representation format for Euclidean Jordan algebra
752 But, when we construct an algebra from a basis of matrices,
753 those matrix representations are lost in favor of coordinate
754 vectors *with respect to* that basis. We could eventually
755 convert back if we tried hard enough, but having the original
756 representations handy is valuable enough that we simply store
757 them and return them from this method.
759 Why implement this for non-matrix algebras? Avoiding special
760 cases for the :class:`BilinearFormEJA` pays with simplicity in
761 its own right. But mainly, we would like to be able to assume
762 that elements of a :class:`CartesianProductEJA` can be displayed
763 nicely, without having to have special classes for direct sums
764 one of whose components was a matrix algebra.
768 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
769 ....: RealSymmetricEJA)
773 sage: J = RealSymmetricEJA(2)
775 Finite family {0: e0, 1: e1, 2: e2}
776 sage: J.matrix_basis()
778 [1 0] [ 0 0.7071067811865475?] [0 0]
779 [0 0], [0.7071067811865475? 0], [0 1]
784 sage: J = JordanSpinEJA(2)
786 Finite family {0: e0, 1: e1}
787 sage: J.matrix_basis()
793 return self
._matrix
_basis
796 def matrix_space(self
):
798 Return the matrix space in which this algebra's elements live, if
799 we think of them as matrices (including column vectors of the
802 Generally this will be an `n`-by-`1` column-vector space,
803 except when the algebra is trivial. There it's `n`-by-`n`
804 (where `n` is zero), to ensure that two elements of the matrix
805 space (empty matrices) can be multiplied.
807 Matrix algebras override this with something more useful.
809 if self
.is_trivial():
810 return MatrixSpace(self
.base_ring(), 0)
812 return self
.matrix_basis()[0].parent()
818 Return the unit element of this algebra.
822 sage: from mjo.eja.eja_algebra import (HadamardEJA,
827 We can compute unit element in the Hadamard EJA::
829 sage: J = HadamardEJA(5)
831 e0 + e1 + e2 + e3 + e4
833 The unit element in the Hadamard EJA is inherited in the
834 subalgebras generated by its elements::
836 sage: J = HadamardEJA(5)
838 e0 + e1 + e2 + e3 + e4
839 sage: x = sum(J.gens())
840 sage: A = x.subalgebra_generated_by(orthonormalize=False)
843 sage: A.one().superalgebra_element()
844 e0 + e1 + e2 + e3 + e4
848 The identity element acts like the identity, regardless of
849 whether or not we orthonormalize::
851 sage: set_random_seed()
852 sage: J = random_eja()
853 sage: x = J.random_element()
854 sage: J.one()*x == x and x*J.one() == x
856 sage: A = x.subalgebra_generated_by()
857 sage: y = A.random_element()
858 sage: A.one()*y == y and y*A.one() == y
863 sage: set_random_seed()
864 sage: J = random_eja(field=QQ, orthonormalize=False)
865 sage: x = J.random_element()
866 sage: J.one()*x == x and x*J.one() == x
868 sage: A = x.subalgebra_generated_by(orthonormalize=False)
869 sage: y = A.random_element()
870 sage: A.one()*y == y and y*A.one() == y
873 The matrix of the unit element's operator is the identity,
874 regardless of the base field and whether or not we
877 sage: set_random_seed()
878 sage: J = random_eja()
879 sage: actual = J.one().operator().matrix()
880 sage: expected = matrix.identity(J.base_ring(), J.dimension())
881 sage: actual == expected
883 sage: x = J.random_element()
884 sage: A = x.subalgebra_generated_by()
885 sage: actual = A.one().operator().matrix()
886 sage: expected = matrix.identity(A.base_ring(), A.dimension())
887 sage: actual == expected
892 sage: set_random_seed()
893 sage: J = random_eja(field=QQ, orthonormalize=False)
894 sage: actual = J.one().operator().matrix()
895 sage: expected = matrix.identity(J.base_ring(), J.dimension())
896 sage: actual == expected
898 sage: x = J.random_element()
899 sage: A = x.subalgebra_generated_by(orthonormalize=False)
900 sage: actual = A.one().operator().matrix()
901 sage: expected = matrix.identity(A.base_ring(), A.dimension())
902 sage: actual == expected
905 Ensure that the cached unit element (often precomputed by
906 hand) agrees with the computed one::
908 sage: set_random_seed()
909 sage: J = random_eja()
910 sage: cached = J.one()
911 sage: J.one.clear_cache()
912 sage: J.one() == cached
917 sage: set_random_seed()
918 sage: J = random_eja(field=QQ, orthonormalize=False)
919 sage: cached = J.one()
920 sage: J.one.clear_cache()
921 sage: J.one() == cached
925 # We can brute-force compute the matrices of the operators
926 # that correspond to the basis elements of this algebra.
927 # If some linear combination of those basis elements is the
928 # algebra identity, then the same linear combination of
929 # their matrices has to be the identity matrix.
931 # Of course, matrices aren't vectors in sage, so we have to
932 # appeal to the "long vectors" isometry.
933 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
935 # Now we use basic linear algebra to find the coefficients,
936 # of the matrices-as-vectors-linear-combination, which should
937 # work for the original algebra basis too.
938 A
= matrix(self
.base_ring(), oper_vecs
)
940 # We used the isometry on the left-hand side already, but we
941 # still need to do it for the right-hand side. Recall that we
942 # wanted something that summed to the identity matrix.
943 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
945 # Now if there's an identity element in the algebra, this
946 # should work. We solve on the left to avoid having to
947 # transpose the matrix "A".
948 return self
.from_vector(A
.solve_left(b
))
951 def peirce_decomposition(self
, c
):
953 The Peirce decomposition of this algebra relative to the
956 In the future, this can be extended to a complete system of
957 orthogonal idempotents.
961 - ``c`` -- an idempotent of this algebra.
965 A triple (J0, J5, J1) containing two subalgebras and one subspace
968 - ``J0`` -- the algebra on the eigenspace of ``c.operator()``
969 corresponding to the eigenvalue zero.
971 - ``J5`` -- the eigenspace (NOT a subalgebra) of ``c.operator()``
972 corresponding to the eigenvalue one-half.
974 - ``J1`` -- the algebra on the eigenspace of ``c.operator()``
975 corresponding to the eigenvalue one.
977 These are the only possible eigenspaces for that operator, and this
978 algebra is a direct sum of them. The spaces ``J0`` and ``J1`` are
979 orthogonal, and are subalgebras of this algebra with the appropriate
984 sage: from mjo.eja.eja_algebra import random_eja, RealSymmetricEJA
988 The canonical example comes from the symmetric matrices, which
989 decompose into diagonal and off-diagonal parts::
991 sage: J = RealSymmetricEJA(3)
992 sage: C = matrix(QQ, [ [1,0,0],
996 sage: J0,J5,J1 = J.peirce_decomposition(c)
998 Euclidean Jordan algebra of dimension 1...
1000 Vector space of degree 6 and dimension 2...
1002 Euclidean Jordan algebra of dimension 3...
1003 sage: J0.one().to_matrix()
1007 sage: orig_df = AA.options.display_format
1008 sage: AA.options.display_format = 'radical'
1009 sage: J.from_vector(J5.basis()[0]).to_matrix()
1013 sage: J.from_vector(J5.basis()[1]).to_matrix()
1017 sage: AA.options.display_format = orig_df
1018 sage: J1.one().to_matrix()
1025 Every algebra decomposes trivially with respect to its identity
1028 sage: set_random_seed()
1029 sage: J = random_eja()
1030 sage: J0,J5,J1 = J.peirce_decomposition(J.one())
1031 sage: J0.dimension() == 0 and J5.dimension() == 0
1033 sage: J1.superalgebra() == J and J1.dimension() == J.dimension()
1036 The decomposition is into eigenspaces, and its components are
1037 therefore necessarily orthogonal. Moreover, the identity
1038 elements in the two subalgebras are the projections onto their
1039 respective subspaces of the superalgebra's identity element::
1041 sage: set_random_seed()
1042 sage: J = random_eja()
1043 sage: x = J.random_element()
1044 sage: if not J.is_trivial():
1045 ....: while x.is_nilpotent():
1046 ....: x = J.random_element()
1047 sage: c = x.subalgebra_idempotent()
1048 sage: J0,J5,J1 = J.peirce_decomposition(c)
1050 sage: for (w,y,z) in zip(J0.basis(), J5.basis(), J1.basis()):
1051 ....: w = w.superalgebra_element()
1052 ....: y = J.from_vector(y)
1053 ....: z = z.superalgebra_element()
1054 ....: ipsum += w.inner_product(y).abs()
1055 ....: ipsum += w.inner_product(z).abs()
1056 ....: ipsum += y.inner_product(z).abs()
1059 sage: J1(c) == J1.one()
1061 sage: J0(J.one() - c) == J0.one()
1065 if not c
.is_idempotent():
1066 raise ValueError("element is not idempotent: %s" % c
)
1068 # Default these to what they should be if they turn out to be
1069 # trivial, because eigenspaces_left() won't return eigenvalues
1070 # corresponding to trivial spaces (e.g. it returns only the
1071 # eigenspace corresponding to lambda=1 if you take the
1072 # decomposition relative to the identity element).
1073 trivial
= self
.subalgebra(())
1074 J0
= trivial
# eigenvalue zero
1075 J5
= VectorSpace(self
.base_ring(), 0) # eigenvalue one-half
1076 J1
= trivial
# eigenvalue one
1078 for (eigval
, eigspace
) in c
.operator().matrix().right_eigenspaces():
1079 if eigval
== ~
(self
.base_ring()(2)):
1082 gens
= tuple( self
.from_vector(b
) for b
in eigspace
.basis() )
1083 subalg
= self
.subalgebra(gens
, check_axioms
=False)
1089 raise ValueError("unexpected eigenvalue: %s" % eigval
)
1094 def random_element(self
, thorough
=False):
1096 Return a random element of this algebra.
1098 Our algebra superclass method only returns a linear
1099 combination of at most two basis elements. We instead
1100 want the vector space "random element" method that
1101 returns a more diverse selection.
1105 - ``thorough`` -- (boolean; default False) whether or not we
1106 should generate irrational coefficients for the random
1107 element when our base ring is irrational; this slows the
1108 algebra operations to a crawl, but any truly random method
1112 # For a general base ring... maybe we can trust this to do the
1113 # right thing? Unlikely, but.
1114 V
= self
.vector_space()
1115 v
= V
.random_element()
1117 if self
.base_ring() is AA
:
1118 # The "random element" method of the algebraic reals is
1119 # stupid at the moment, and only returns integers between
1120 # -2 and 2, inclusive:
1122 # https://trac.sagemath.org/ticket/30875
1124 # Instead, we implement our own "random vector" method,
1125 # and then coerce that into the algebra. We use the vector
1126 # space degree here instead of the dimension because a
1127 # subalgebra could (for example) be spanned by only two
1128 # vectors, each with five coordinates. We need to
1129 # generate all five coordinates.
1131 v
*= QQbar
.random_element().real()
1133 v
*= QQ
.random_element()
1135 return self
.from_vector(V
.coordinate_vector(v
))
1137 def random_elements(self
, count
, thorough
=False):
1139 Return ``count`` random elements as a tuple.
1143 - ``thorough`` -- (boolean; default False) whether or not we
1144 should generate irrational coefficients for the random
1145 elements when our base ring is irrational; this slows the
1146 algebra operations to a crawl, but any truly random method
1151 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1155 sage: J = JordanSpinEJA(3)
1156 sage: x,y,z = J.random_elements(3)
1157 sage: all( [ x in J, y in J, z in J ])
1159 sage: len( J.random_elements(10) ) == 10
1163 return tuple( self
.random_element(thorough
)
1164 for idx
in range(count
) )
1168 def _charpoly_coefficients(self
):
1170 The `r` polynomial coefficients of the "characteristic polynomial
1175 sage: from mjo.eja.eja_algebra import random_eja
1179 The theory shows that these are all homogeneous polynomials of
1182 sage: set_random_seed()
1183 sage: J = random_eja()
1184 sage: all(p.is_homogeneous() for p in J._charpoly_coefficients())
1188 n
= self
.dimension()
1189 R
= self
.coordinate_polynomial_ring()
1191 F
= R
.fraction_field()
1194 # From a result in my book, these are the entries of the
1195 # basis representation of L_x.
1196 return sum( vars[k
]*self
.gens()[k
].operator().matrix()[i
,j
]
1199 L_x
= matrix(F
, n
, n
, L_x_i_j
)
1202 if self
.rank
.is_in_cache():
1204 # There's no need to pad the system with redundant
1205 # columns if we *know* they'll be redundant.
1208 # Compute an extra power in case the rank is equal to
1209 # the dimension (otherwise, we would stop at x^(r-1)).
1210 x_powers
= [ (L_x
**k
)*self
.one().to_vector()
1211 for k
in range(n
+1) ]
1212 A
= matrix
.column(F
, x_powers
[:n
])
1213 AE
= A
.extended_echelon_form()
1220 # The theory says that only the first "r" coefficients are
1221 # nonzero, and they actually live in the original polynomial
1222 # ring and not the fraction field. We negate them because in
1223 # the actual characteristic polynomial, they get moved to the
1224 # other side where x^r lives. We don't bother to trim A_rref
1225 # down to a square matrix and solve the resulting system,
1226 # because the upper-left r-by-r portion of A_rref is
1227 # guaranteed to be the identity matrix, so e.g.
1229 # A_rref.solve_right(Y)
1231 # would just be returning Y.
1232 return (-E
*b
)[:r
].change_ring(R
)
1237 Return the rank of this EJA.
1239 This is a cached method because we know the rank a priori for
1240 all of the algebras we can construct. Thus we can avoid the
1241 expensive ``_charpoly_coefficients()`` call unless we truly
1242 need to compute the whole characteristic polynomial.
1246 sage: from mjo.eja.eja_algebra import (HadamardEJA,
1247 ....: JordanSpinEJA,
1248 ....: RealSymmetricEJA,
1249 ....: ComplexHermitianEJA,
1250 ....: QuaternionHermitianEJA,
1255 The rank of the Jordan spin algebra is always two::
1257 sage: JordanSpinEJA(2).rank()
1259 sage: JordanSpinEJA(3).rank()
1261 sage: JordanSpinEJA(4).rank()
1264 The rank of the `n`-by-`n` Hermitian real, complex, or
1265 quaternion matrices is `n`::
1267 sage: RealSymmetricEJA(4).rank()
1269 sage: ComplexHermitianEJA(3).rank()
1271 sage: QuaternionHermitianEJA(2).rank()
1276 Ensure that every EJA that we know how to construct has a
1277 positive integer rank, unless the algebra is trivial in
1278 which case its rank will be zero::
1280 sage: set_random_seed()
1281 sage: J = random_eja()
1285 sage: r > 0 or (r == 0 and J.is_trivial())
1288 Ensure that computing the rank actually works, since the ranks
1289 of all simple algebras are known and will be cached by default::
1291 sage: set_random_seed() # long time
1292 sage: J = random_eja() # long time
1293 sage: cached = J.rank() # long time
1294 sage: J.rank.clear_cache() # long time
1295 sage: J.rank() == cached # long time
1299 return len(self
._charpoly
_coefficients
())
1302 def subalgebra(self
, basis
, **kwargs
):
1304 Create a subalgebra of this algebra from the given basis.
1306 from mjo
.eja
.eja_subalgebra
import FiniteDimensionalEJASubalgebra
1307 return FiniteDimensionalEJASubalgebra(self
, basis
, **kwargs
)
1310 def vector_space(self
):
1312 Return the vector space that underlies this algebra.
1316 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1320 sage: J = RealSymmetricEJA(2)
1321 sage: J.vector_space()
1322 Vector space of dimension 3 over...
1325 return self
.zero().to_vector().parent().ambient_vector_space()
1329 class RationalBasisEJA(FiniteDimensionalEJA
):
1331 New class for algebras whose supplied basis elements have all rational entries.
1335 sage: from mjo.eja.eja_algebra import BilinearFormEJA
1339 The supplied basis is orthonormalized by default::
1341 sage: B = matrix(QQ, [[1, 0, 0], [0, 25, -32], [0, -32, 41]])
1342 sage: J = BilinearFormEJA(B)
1343 sage: J.matrix_basis()
1360 # Abuse the check_field parameter to check that the entries of
1361 # out basis (in ambient coordinates) are in the field QQ.
1362 if not all( all(b_i
in QQ
for b_i
in b
.list()) for b
in basis
):
1363 raise TypeError("basis not rational")
1365 self
._rational
_algebra
= None
1367 # There's no point in constructing the extra algebra if this
1368 # one is already rational.
1370 # Note: the same Jordan and inner-products work here,
1371 # because they are necessarily defined with respect to
1372 # ambient coordinates and not any particular basis.
1373 self
._rational
_algebra
= FiniteDimensionalEJA(
1378 orthonormalize
=False,
1382 super().__init
__(basis
,
1386 check_field
=check_field
,
1390 def _charpoly_coefficients(self
):
1394 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
1395 ....: JordanSpinEJA)
1399 The base ring of the resulting polynomial coefficients is what
1400 it should be, and not the rationals (unless the algebra was
1401 already over the rationals)::
1403 sage: J = JordanSpinEJA(3)
1404 sage: J._charpoly_coefficients()
1405 (X1^2 - X2^2 - X3^2, -2*X1)
1406 sage: a0 = J._charpoly_coefficients()[0]
1408 Algebraic Real Field
1409 sage: a0.base_ring()
1410 Algebraic Real Field
1413 if self
._rational
_algebra
is None:
1414 # There's no need to construct *another* algebra over the
1415 # rationals if this one is already over the
1416 # rationals. Likewise, if we never orthonormalized our
1417 # basis, we might as well just use the given one.
1418 return super()._charpoly
_coefficients
()
1420 # Do the computation over the rationals. The answer will be
1421 # the same, because all we've done is a change of basis.
1422 # Then, change back from QQ to our real base ring
1423 a
= ( a_i
.change_ring(self
.base_ring())
1424 for a_i
in self
._rational
_algebra
._charpoly
_coefficients
() )
1426 if self
._deortho
_matrix
is None:
1427 # This can happen if our base ring was, say, AA and we
1428 # chose not to (or didn't need to) orthonormalize. It's
1429 # still faster to do the computations over QQ even if
1430 # the numbers in the boxes stay the same.
1433 # Otherwise, convert the coordinate variables back to the
1434 # deorthonormalized ones.
1435 R
= self
.coordinate_polynomial_ring()
1436 from sage
.modules
.free_module_element
import vector
1437 X
= vector(R
, R
.gens())
1438 BX
= self
._deortho
_matrix
*X
1440 subs_dict
= { X[i]: BX[i] for i in range(len(X)) }
1441 return tuple( a_i
.subs(subs_dict
) for a_i
in a
)
1443 class ConcreteEJA(RationalBasisEJA
):
1445 A class for the Euclidean Jordan algebras that we know by name.
1447 These are the Jordan algebras whose basis, multiplication table,
1448 rank, and so on are known a priori. More to the point, they are
1449 the Euclidean Jordan algebras for which we are able to conjure up
1450 a "random instance."
1454 sage: from mjo.eja.eja_algebra import ConcreteEJA
1458 Our basis is normalized with respect to the algebra's inner
1459 product, unless we specify otherwise::
1461 sage: set_random_seed()
1462 sage: J = ConcreteEJA.random_instance()
1463 sage: all( b.norm() == 1 for b in J.gens() )
1466 Since our basis is orthonormal with respect to the algebra's inner
1467 product, and since we know that this algebra is an EJA, any
1468 left-multiplication operator's matrix will be symmetric because
1469 natural->EJA basis representation is an isometry and within the
1470 EJA the operator is self-adjoint by the Jordan axiom::
1472 sage: set_random_seed()
1473 sage: J = ConcreteEJA.random_instance()
1474 sage: x = J.random_element()
1475 sage: x.operator().is_self_adjoint()
1480 def _max_random_instance_size():
1482 Return an integer "size" that is an upper bound on the size of
1483 this algebra when it is used in a random test
1484 case. Unfortunately, the term "size" is ambiguous -- when
1485 dealing with `R^n` under either the Hadamard or Jordan spin
1486 product, the "size" refers to the dimension `n`. When dealing
1487 with a matrix algebra (real symmetric or complex/quaternion
1488 Hermitian), it refers to the size of the matrix, which is far
1489 less than the dimension of the underlying vector space.
1491 This method must be implemented in each subclass.
1493 raise NotImplementedError
1496 def random_instance(cls
, *args
, **kwargs
):
1498 Return a random instance of this type of algebra.
1500 This method should be implemented in each subclass.
1502 from sage
.misc
.prandom
import choice
1503 eja_class
= choice(cls
.__subclasses
__())
1505 # These all bubble up to the RationalBasisEJA superclass
1506 # constructor, so any (kw)args valid there are also valid
1508 return eja_class
.random_instance(*args
, **kwargs
)
1513 def dimension_over_reals():
1515 The dimension of this matrix's base ring over the reals.
1517 The reals are dimension one over themselves, obviously; that's
1518 just `\mathbb{R}^{1}`. Likewise, the complex numbers `a + bi`
1519 have dimension two. Finally, the quaternions have dimension
1520 four over the reals.
1522 This is used to determine the size of the matrix returned from
1523 :meth:`real_embed`, among other things.
1525 raise NotImplementedError
1528 def real_embed(cls
,M
):
1530 Embed the matrix ``M`` into a space of real matrices.
1532 The matrix ``M`` can have entries in any field at the moment:
1533 the real numbers, complex numbers, or quaternions. And although
1534 they are not a field, we can probably support octonions at some
1535 point, too. This function returns a real matrix that "acts like"
1536 the original with respect to matrix multiplication; i.e.
1538 real_embed(M*N) = real_embed(M)*real_embed(N)
1541 if M
.ncols() != M
.nrows():
1542 raise ValueError("the matrix 'M' must be square")
1547 def real_unembed(cls
,M
):
1549 The inverse of :meth:`real_embed`.
1551 if M
.ncols() != M
.nrows():
1552 raise ValueError("the matrix 'M' must be square")
1553 if not ZZ(M
.nrows()).mod(cls
.dimension_over_reals()).is_zero():
1554 raise ValueError("the matrix 'M' must be a real embedding")
1558 def jordan_product(X
,Y
):
1559 return (X
*Y
+ Y
*X
)/2
1562 def trace_inner_product(cls
,X
,Y
):
1564 Compute the trace inner-product of two real-embeddings.
1568 sage: from mjo.eja.eja_algebra import (RealSymmetricEJA,
1569 ....: ComplexHermitianEJA,
1570 ....: QuaternionHermitianEJA)
1574 This gives the same answer as it would if we computed the trace
1575 from the unembedded (original) matrices::
1577 sage: set_random_seed()
1578 sage: J = RealSymmetricEJA.random_instance()
1579 sage: x,y = J.random_elements(2)
1580 sage: Xe = x.to_matrix()
1581 sage: Ye = y.to_matrix()
1582 sage: X = J.real_unembed(Xe)
1583 sage: Y = J.real_unembed(Ye)
1584 sage: expected = (X*Y).trace()
1585 sage: actual = J.trace_inner_product(Xe,Ye)
1586 sage: actual == expected
1591 sage: set_random_seed()
1592 sage: J = ComplexHermitianEJA.random_instance()
1593 sage: x,y = J.random_elements(2)
1594 sage: Xe = x.to_matrix()
1595 sage: Ye = y.to_matrix()
1596 sage: X = J.real_unembed(Xe)
1597 sage: Y = J.real_unembed(Ye)
1598 sage: expected = (X*Y).trace().real()
1599 sage: actual = J.trace_inner_product(Xe,Ye)
1600 sage: actual == expected
1605 sage: set_random_seed()
1606 sage: J = QuaternionHermitianEJA.random_instance()
1607 sage: x,y = J.random_elements(2)
1608 sage: Xe = x.to_matrix()
1609 sage: Ye = y.to_matrix()
1610 sage: X = J.real_unembed(Xe)
1611 sage: Y = J.real_unembed(Ye)
1612 sage: expected = (X*Y).trace().coefficient_tuple()[0]
1613 sage: actual = J.trace_inner_product(Xe,Ye)
1614 sage: actual == expected
1618 Xu
= cls
.real_unembed(X
)
1619 Yu
= cls
.real_unembed(Y
)
1620 tr
= (Xu
*Yu
).trace()
1623 # Works in QQ, AA, RDF, et cetera.
1625 except AttributeError:
1626 # A quaternion doesn't have a real() method, but does
1627 # have coefficient_tuple() method that returns the
1628 # coefficients of 1, i, j, and k -- in that order.
1629 return tr
.coefficient_tuple()[0]
1632 class RealMatrixEJA(MatrixEJA
):
1634 def dimension_over_reals():
1638 class RealSymmetricEJA(ConcreteEJA
, RealMatrixEJA
):
1640 The rank-n simple EJA consisting of real symmetric n-by-n
1641 matrices, the usual symmetric Jordan product, and the trace inner
1642 product. It has dimension `(n^2 + n)/2` over the reals.
1646 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1650 sage: J = RealSymmetricEJA(2)
1651 sage: e0, e1, e2 = J.gens()
1659 In theory, our "field" can be any subfield of the reals::
1661 sage: RealSymmetricEJA(2, field=RDF)
1662 Euclidean Jordan algebra of dimension 3 over Real Double Field
1663 sage: RealSymmetricEJA(2, field=RR)
1664 Euclidean Jordan algebra of dimension 3 over Real Field with
1665 53 bits of precision
1669 The dimension of this algebra is `(n^2 + n) / 2`::
1671 sage: set_random_seed()
1672 sage: n_max = RealSymmetricEJA._max_random_instance_size()
1673 sage: n = ZZ.random_element(1, n_max)
1674 sage: J = RealSymmetricEJA(n)
1675 sage: J.dimension() == (n^2 + n)/2
1678 The Jordan multiplication is what we think it is::
1680 sage: set_random_seed()
1681 sage: J = RealSymmetricEJA.random_instance()
1682 sage: x,y = J.random_elements(2)
1683 sage: actual = (x*y).to_matrix()
1684 sage: X = x.to_matrix()
1685 sage: Y = y.to_matrix()
1686 sage: expected = (X*Y + Y*X)/2
1687 sage: actual == expected
1689 sage: J(expected) == x*y
1692 We can change the generator prefix::
1694 sage: RealSymmetricEJA(3, prefix='q').gens()
1695 (q0, q1, q2, q3, q4, q5)
1697 We can construct the (trivial) algebra of rank zero::
1699 sage: RealSymmetricEJA(0)
1700 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1704 def _denormalized_basis(cls
, n
):
1706 Return a basis for the space of real symmetric n-by-n matrices.
1710 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1714 sage: set_random_seed()
1715 sage: n = ZZ.random_element(1,5)
1716 sage: B = RealSymmetricEJA._denormalized_basis(n)
1717 sage: all( M.is_symmetric() for M in B)
1721 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1725 for j
in range(i
+1):
1726 Eij
= matrix(ZZ
, n
, lambda k
,l
: k
==i
and l
==j
)
1730 Sij
= Eij
+ Eij
.transpose()
1736 def _max_random_instance_size():
1737 return 4 # Dimension 10
1740 def random_instance(cls
, **kwargs
):
1742 Return a random instance of this type of algebra.
1744 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1745 return cls(n
, **kwargs
)
1747 def __init__(self
, n
, **kwargs
):
1748 # We know this is a valid EJA, but will double-check
1749 # if the user passes check_axioms=True.
1750 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
1752 super(RealSymmetricEJA
, self
).__init
__(self
._denormalized
_basis
(n
),
1753 self
.jordan_product
,
1754 self
.trace_inner_product
,
1757 # TODO: this could be factored out somehow, but is left here
1758 # because the MatrixEJA is not presently a subclass of the
1759 # FDEJA class that defines rank() and one().
1760 self
.rank
.set_cache(n
)
1761 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
1762 self
.one
.set_cache(self(idV
))
1766 class ComplexMatrixEJA(MatrixEJA
):
1767 # A manual dictionary-cache for the complex_extension() method,
1768 # since apparently @classmethods can't also be @cached_methods.
1769 _complex_extension
= {}
1772 def complex_extension(cls
,field
):
1774 The complex field that we embed/unembed, as an extension
1775 of the given ``field``.
1777 if field
in cls
._complex
_extension
:
1778 return cls
._complex
_extension
[field
]
1780 # Sage doesn't know how to adjoin the complex "i" (the root of
1781 # x^2 + 1) to a field in a general way. Here, we just enumerate
1782 # all of the cases that I have cared to support so far.
1784 # Sage doesn't know how to embed AA into QQbar, i.e. how
1785 # to adjoin sqrt(-1) to AA.
1787 elif not field
.is_exact():
1789 F
= field
.complex_field()
1791 # Works for QQ and... maybe some other fields.
1792 R
= PolynomialRing(field
, 'z')
1794 F
= field
.extension(z
**2 + 1, 'I', embedding
=CLF(-1).sqrt())
1796 cls
._complex
_extension
[field
] = F
1800 def dimension_over_reals():
1804 def real_embed(cls
,M
):
1806 Embed the n-by-n complex matrix ``M`` into the space of real
1807 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1808 bi` to the block matrix ``[[a,b],[-b,a]]``.
1812 sage: from mjo.eja.eja_algebra import ComplexMatrixEJA
1816 sage: F = QuadraticField(-1, 'I')
1817 sage: x1 = F(4 - 2*i)
1818 sage: x2 = F(1 + 2*i)
1821 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1822 sage: ComplexMatrixEJA.real_embed(M)
1831 Embedding is a homomorphism (isomorphism, in fact)::
1833 sage: set_random_seed()
1834 sage: n = ZZ.random_element(3)
1835 sage: F = QuadraticField(-1, 'I')
1836 sage: X = random_matrix(F, n)
1837 sage: Y = random_matrix(F, n)
1838 sage: Xe = ComplexMatrixEJA.real_embed(X)
1839 sage: Ye = ComplexMatrixEJA.real_embed(Y)
1840 sage: XYe = ComplexMatrixEJA.real_embed(X*Y)
1845 super(ComplexMatrixEJA
,cls
).real_embed(M
)
1848 # We don't need any adjoined elements...
1849 field
= M
.base_ring().base_ring()
1855 blocks
.append(matrix(field
, 2, [ [ a
, b
],
1858 return matrix
.block(field
, n
, blocks
)
1862 def real_unembed(cls
,M
):
1864 The inverse of _embed_complex_matrix().
1868 sage: from mjo.eja.eja_algebra import ComplexMatrixEJA
1872 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1873 ....: [-2, 1, -4, 3],
1874 ....: [ 9, 10, 11, 12],
1875 ....: [-10, 9, -12, 11] ])
1876 sage: ComplexMatrixEJA.real_unembed(A)
1878 [ 10*I + 9 12*I + 11]
1882 Unembedding is the inverse of embedding::
1884 sage: set_random_seed()
1885 sage: F = QuadraticField(-1, 'I')
1886 sage: M = random_matrix(F, 3)
1887 sage: Me = ComplexMatrixEJA.real_embed(M)
1888 sage: ComplexMatrixEJA.real_unembed(Me) == M
1892 super(ComplexMatrixEJA
,cls
).real_unembed(M
)
1894 d
= cls
.dimension_over_reals()
1895 F
= cls
.complex_extension(M
.base_ring())
1898 # Go top-left to bottom-right (reading order), converting every
1899 # 2-by-2 block we see to a single complex element.
1901 for k
in range(n
/d
):
1902 for j
in range(n
/d
):
1903 submat
= M
[d
*k
:d
*k
+d
,d
*j
:d
*j
+d
]
1904 if submat
[0,0] != submat
[1,1]:
1905 raise ValueError('bad on-diagonal submatrix')
1906 if submat
[0,1] != -submat
[1,0]:
1907 raise ValueError('bad off-diagonal submatrix')
1908 z
= submat
[0,0] + submat
[0,1]*i
1911 return matrix(F
, n
/d
, elements
)
1914 class ComplexHermitianEJA(ConcreteEJA
, ComplexMatrixEJA
):
1916 The rank-n simple EJA consisting of complex Hermitian n-by-n
1917 matrices over the real numbers, the usual symmetric Jordan product,
1918 and the real-part-of-trace inner product. It has dimension `n^2` over
1923 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1927 In theory, our "field" can be any subfield of the reals::
1929 sage: ComplexHermitianEJA(2, field=RDF)
1930 Euclidean Jordan algebra of dimension 4 over Real Double Field
1931 sage: ComplexHermitianEJA(2, field=RR)
1932 Euclidean Jordan algebra of dimension 4 over Real Field with
1933 53 bits of precision
1937 The dimension of this algebra is `n^2`::
1939 sage: set_random_seed()
1940 sage: n_max = ComplexHermitianEJA._max_random_instance_size()
1941 sage: n = ZZ.random_element(1, n_max)
1942 sage: J = ComplexHermitianEJA(n)
1943 sage: J.dimension() == n^2
1946 The Jordan multiplication is what we think it is::
1948 sage: set_random_seed()
1949 sage: J = ComplexHermitianEJA.random_instance()
1950 sage: x,y = J.random_elements(2)
1951 sage: actual = (x*y).to_matrix()
1952 sage: X = x.to_matrix()
1953 sage: Y = y.to_matrix()
1954 sage: expected = (X*Y + Y*X)/2
1955 sage: actual == expected
1957 sage: J(expected) == x*y
1960 We can change the generator prefix::
1962 sage: ComplexHermitianEJA(2, prefix='z').gens()
1965 We can construct the (trivial) algebra of rank zero::
1967 sage: ComplexHermitianEJA(0)
1968 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1973 def _denormalized_basis(cls
, n
):
1975 Returns a basis for the space of complex Hermitian n-by-n matrices.
1977 Why do we embed these? Basically, because all of numerical linear
1978 algebra assumes that you're working with vectors consisting of `n`
1979 entries from a field and scalars from the same field. There's no way
1980 to tell SageMath that (for example) the vectors contain complex
1981 numbers, while the scalar field is real.
1985 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1989 sage: set_random_seed()
1990 sage: n = ZZ.random_element(1,5)
1991 sage: B = ComplexHermitianEJA._denormalized_basis(n)
1992 sage: all( M.is_symmetric() for M in B)
1997 R
= PolynomialRing(field
, 'z')
1999 F
= field
.extension(z
**2 + 1, 'I')
2002 # This is like the symmetric case, but we need to be careful:
2004 # * We want conjugate-symmetry, not just symmetry.
2005 # * The diagonal will (as a result) be real.
2008 Eij
= matrix
.zero(F
,n
)
2010 for j
in range(i
+1):
2014 Sij
= cls
.real_embed(Eij
)
2017 # The second one has a minus because it's conjugated.
2018 Eij
[j
,i
] = 1 # Eij = Eij + Eij.transpose()
2019 Sij_real
= cls
.real_embed(Eij
)
2021 # Eij = I*Eij - I*Eij.transpose()
2024 Sij_imag
= cls
.real_embed(Eij
)
2030 # Since we embedded these, we can drop back to the "field" that we
2031 # started with instead of the complex extension "F".
2032 return tuple( s
.change_ring(field
) for s
in S
)
2035 def __init__(self
, n
, **kwargs
):
2036 # We know this is a valid EJA, but will double-check
2037 # if the user passes check_axioms=True.
2038 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2040 super(ComplexHermitianEJA
, self
).__init
__(self
._denormalized
_basis
(n
),
2041 self
.jordan_product
,
2042 self
.trace_inner_product
,
2044 # TODO: this could be factored out somehow, but is left here
2045 # because the MatrixEJA is not presently a subclass of the
2046 # FDEJA class that defines rank() and one().
2047 self
.rank
.set_cache(n
)
2048 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
2049 self
.one
.set_cache(self(idV
))
2052 def _max_random_instance_size():
2053 return 3 # Dimension 9
2056 def random_instance(cls
, **kwargs
):
2058 Return a random instance of this type of algebra.
2060 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2061 return cls(n
, **kwargs
)
2063 class QuaternionMatrixEJA(MatrixEJA
):
2065 # A manual dictionary-cache for the quaternion_extension() method,
2066 # since apparently @classmethods can't also be @cached_methods.
2067 _quaternion_extension
= {}
2070 def quaternion_extension(cls
,field
):
2072 The quaternion field that we embed/unembed, as an extension
2073 of the given ``field``.
2075 if field
in cls
._quaternion
_extension
:
2076 return cls
._quaternion
_extension
[field
]
2078 Q
= QuaternionAlgebra(field
,-1,-1)
2080 cls
._quaternion
_extension
[field
] = Q
2084 def dimension_over_reals():
2088 def real_embed(cls
,M
):
2090 Embed the n-by-n quaternion matrix ``M`` into the space of real
2091 matrices of size 4n-by-4n by first sending each quaternion entry `z
2092 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
2093 c+di],[-c + di, a-bi]]`, and then embedding those into a real
2098 sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA
2102 sage: Q = QuaternionAlgebra(QQ,-1,-1)
2103 sage: i,j,k = Q.gens()
2104 sage: x = 1 + 2*i + 3*j + 4*k
2105 sage: M = matrix(Q, 1, [[x]])
2106 sage: QuaternionMatrixEJA.real_embed(M)
2112 Embedding is a homomorphism (isomorphism, in fact)::
2114 sage: set_random_seed()
2115 sage: n = ZZ.random_element(2)
2116 sage: Q = QuaternionAlgebra(QQ,-1,-1)
2117 sage: X = random_matrix(Q, n)
2118 sage: Y = random_matrix(Q, n)
2119 sage: Xe = QuaternionMatrixEJA.real_embed(X)
2120 sage: Ye = QuaternionMatrixEJA.real_embed(Y)
2121 sage: XYe = QuaternionMatrixEJA.real_embed(X*Y)
2126 super(QuaternionMatrixEJA
,cls
).real_embed(M
)
2127 quaternions
= M
.base_ring()
2130 F
= QuadraticField(-1, 'I')
2135 t
= z
.coefficient_tuple()
2140 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
2141 [-c
+ d
*i
, a
- b
*i
]])
2142 realM
= ComplexMatrixEJA
.real_embed(cplxM
)
2143 blocks
.append(realM
)
2145 # We should have real entries by now, so use the realest field
2146 # we've got for the return value.
2147 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
2152 def real_unembed(cls
,M
):
2154 The inverse of _embed_quaternion_matrix().
2158 sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA
2162 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
2163 ....: [-2, 1, -4, 3],
2164 ....: [-3, 4, 1, -2],
2165 ....: [-4, -3, 2, 1]])
2166 sage: QuaternionMatrixEJA.real_unembed(M)
2167 [1 + 2*i + 3*j + 4*k]
2171 Unembedding is the inverse of embedding::
2173 sage: set_random_seed()
2174 sage: Q = QuaternionAlgebra(QQ, -1, -1)
2175 sage: M = random_matrix(Q, 3)
2176 sage: Me = QuaternionMatrixEJA.real_embed(M)
2177 sage: QuaternionMatrixEJA.real_unembed(Me) == M
2181 super(QuaternionMatrixEJA
,cls
).real_unembed(M
)
2183 d
= cls
.dimension_over_reals()
2185 # Use the base ring of the matrix to ensure that its entries can be
2186 # multiplied by elements of the quaternion algebra.
2187 Q
= cls
.quaternion_extension(M
.base_ring())
2190 # Go top-left to bottom-right (reading order), converting every
2191 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
2194 for l
in range(n
/d
):
2195 for m
in range(n
/d
):
2196 submat
= ComplexMatrixEJA
.real_unembed(
2197 M
[d
*l
:d
*l
+d
,d
*m
:d
*m
+d
] )
2198 if submat
[0,0] != submat
[1,1].conjugate():
2199 raise ValueError('bad on-diagonal submatrix')
2200 if submat
[0,1] != -submat
[1,0].conjugate():
2201 raise ValueError('bad off-diagonal submatrix')
2202 z
= submat
[0,0].real()
2203 z
+= submat
[0,0].imag()*i
2204 z
+= submat
[0,1].real()*j
2205 z
+= submat
[0,1].imag()*k
2208 return matrix(Q
, n
/d
, elements
)
2211 class QuaternionHermitianEJA(ConcreteEJA
, QuaternionMatrixEJA
):
2213 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
2214 matrices, the usual symmetric Jordan product, and the
2215 real-part-of-trace inner product. It has dimension `2n^2 - n` over
2220 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2224 In theory, our "field" can be any subfield of the reals::
2226 sage: QuaternionHermitianEJA(2, field=RDF)
2227 Euclidean Jordan algebra of dimension 6 over Real Double Field
2228 sage: QuaternionHermitianEJA(2, field=RR)
2229 Euclidean Jordan algebra of dimension 6 over Real Field with
2230 53 bits of precision
2234 The dimension of this algebra is `2*n^2 - n`::
2236 sage: set_random_seed()
2237 sage: n_max = QuaternionHermitianEJA._max_random_instance_size()
2238 sage: n = ZZ.random_element(1, n_max)
2239 sage: J = QuaternionHermitianEJA(n)
2240 sage: J.dimension() == 2*(n^2) - n
2243 The Jordan multiplication is what we think it is::
2245 sage: set_random_seed()
2246 sage: J = QuaternionHermitianEJA.random_instance()
2247 sage: x,y = J.random_elements(2)
2248 sage: actual = (x*y).to_matrix()
2249 sage: X = x.to_matrix()
2250 sage: Y = y.to_matrix()
2251 sage: expected = (X*Y + Y*X)/2
2252 sage: actual == expected
2254 sage: J(expected) == x*y
2257 We can change the generator prefix::
2259 sage: QuaternionHermitianEJA(2, prefix='a').gens()
2260 (a0, a1, a2, a3, a4, a5)
2262 We can construct the (trivial) algebra of rank zero::
2264 sage: QuaternionHermitianEJA(0)
2265 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2269 def _denormalized_basis(cls
, n
):
2271 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
2273 Why do we embed these? Basically, because all of numerical
2274 linear algebra assumes that you're working with vectors consisting
2275 of `n` entries from a field and scalars from the same field. There's
2276 no way to tell SageMath that (for example) the vectors contain
2277 complex numbers, while the scalar field is real.
2281 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2285 sage: set_random_seed()
2286 sage: n = ZZ.random_element(1,5)
2287 sage: B = QuaternionHermitianEJA._denormalized_basis(n)
2288 sage: all( M.is_symmetric() for M in B )
2293 Q
= QuaternionAlgebra(QQ
,-1,-1)
2296 # This is like the symmetric case, but we need to be careful:
2298 # * We want conjugate-symmetry, not just symmetry.
2299 # * The diagonal will (as a result) be real.
2302 Eij
= matrix
.zero(Q
,n
)
2304 for j
in range(i
+1):
2308 Sij
= cls
.real_embed(Eij
)
2311 # The second, third, and fourth ones have a minus
2312 # because they're conjugated.
2313 # Eij = Eij + Eij.transpose()
2315 Sij_real
= cls
.real_embed(Eij
)
2317 # Eij = I*(Eij - Eij.transpose())
2320 Sij_I
= cls
.real_embed(Eij
)
2322 # Eij = J*(Eij - Eij.transpose())
2325 Sij_J
= cls
.real_embed(Eij
)
2327 # Eij = K*(Eij - Eij.transpose())
2330 Sij_K
= cls
.real_embed(Eij
)
2336 # Since we embedded these, we can drop back to the "field" that we
2337 # started with instead of the quaternion algebra "Q".
2338 return tuple( s
.change_ring(field
) for s
in S
)
2341 def __init__(self
, n
, **kwargs
):
2342 # We know this is a valid EJA, but will double-check
2343 # if the user passes check_axioms=True.
2344 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2346 super(QuaternionHermitianEJA
, self
).__init
__(self
._denormalized
_basis
(n
),
2347 self
.jordan_product
,
2348 self
.trace_inner_product
,
2350 # TODO: this could be factored out somehow, but is left here
2351 # because the MatrixEJA is not presently a subclass of the
2352 # FDEJA class that defines rank() and one().
2353 self
.rank
.set_cache(n
)
2354 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
2355 self
.one
.set_cache(self(idV
))
2359 def _max_random_instance_size():
2361 The maximum rank of a random QuaternionHermitianEJA.
2363 return 2 # Dimension 6
2366 def random_instance(cls
, **kwargs
):
2368 Return a random instance of this type of algebra.
2370 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2371 return cls(n
, **kwargs
)
2374 class HadamardEJA(ConcreteEJA
):
2376 Return the Euclidean Jordan Algebra corresponding to the set
2377 `R^n` under the Hadamard product.
2379 Note: this is nothing more than the Cartesian product of ``n``
2380 copies of the spin algebra. Once Cartesian product algebras
2381 are implemented, this can go.
2385 sage: from mjo.eja.eja_algebra import HadamardEJA
2389 This multiplication table can be verified by hand::
2391 sage: J = HadamardEJA(3)
2392 sage: e0,e1,e2 = J.gens()
2408 We can change the generator prefix::
2410 sage: HadamardEJA(3, prefix='r').gens()
2414 def __init__(self
, n
, **kwargs
):
2416 jordan_product
= lambda x
,y
: x
2417 inner_product
= lambda x
,y
: x
2419 def jordan_product(x
,y
):
2421 return P( xi
*yi
for (xi
,yi
) in zip(x
,y
) )
2423 def inner_product(x
,y
):
2426 # New defaults for keyword arguments. Don't orthonormalize
2427 # because our basis is already orthonormal with respect to our
2428 # inner-product. Don't check the axioms, because we know this
2429 # is a valid EJA... but do double-check if the user passes
2430 # check_axioms=True. Note: we DON'T override the "check_field"
2431 # default here, because the user can pass in a field!
2432 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2433 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2435 column_basis
= tuple( b
.column() for b
in FreeModule(ZZ
, n
).basis() )
2436 super().__init
__(column_basis
,
2441 self
.rank
.set_cache(n
)
2444 self
.one
.set_cache( self
.zero() )
2446 self
.one
.set_cache( sum(self
.gens()) )
2449 def _max_random_instance_size():
2451 The maximum dimension of a random HadamardEJA.
2456 def random_instance(cls
, **kwargs
):
2458 Return a random instance of this type of algebra.
2460 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2461 return cls(n
, **kwargs
)
2464 class BilinearFormEJA(ConcreteEJA
):
2466 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2467 with the half-trace inner product and jordan product ``x*y =
2468 (<Bx,y>,y_bar>, x0*y_bar + y0*x_bar)`` where `B = 1 \times B22` is
2469 a symmetric positive-definite "bilinear form" matrix. Its
2470 dimension is the size of `B`, and it has rank two in dimensions
2471 larger than two. It reduces to the ``JordanSpinEJA`` when `B` is
2472 the identity matrix of order ``n``.
2474 We insist that the one-by-one upper-left identity block of `B` be
2475 passed in as well so that we can be passed a matrix of size zero
2476 to construct a trivial algebra.
2480 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
2481 ....: JordanSpinEJA)
2485 When no bilinear form is specified, the identity matrix is used,
2486 and the resulting algebra is the Jordan spin algebra::
2488 sage: B = matrix.identity(AA,3)
2489 sage: J0 = BilinearFormEJA(B)
2490 sage: J1 = JordanSpinEJA(3)
2491 sage: J0.multiplication_table() == J0.multiplication_table()
2494 An error is raised if the matrix `B` does not correspond to a
2495 positive-definite bilinear form::
2497 sage: B = matrix.random(QQ,2,3)
2498 sage: J = BilinearFormEJA(B)
2499 Traceback (most recent call last):
2501 ValueError: bilinear form is not positive-definite
2502 sage: B = matrix.zero(QQ,3)
2503 sage: J = BilinearFormEJA(B)
2504 Traceback (most recent call last):
2506 ValueError: bilinear form is not positive-definite
2510 We can create a zero-dimensional algebra::
2512 sage: B = matrix.identity(AA,0)
2513 sage: J = BilinearFormEJA(B)
2517 We can check the multiplication condition given in the Jordan, von
2518 Neumann, and Wigner paper (and also discussed on my "On the
2519 symmetry..." paper). Note that this relies heavily on the standard
2520 choice of basis, as does anything utilizing the bilinear form
2521 matrix. We opt not to orthonormalize the basis, because if we
2522 did, we would have to normalize the `s_{i}` in a similar manner::
2524 sage: set_random_seed()
2525 sage: n = ZZ.random_element(5)
2526 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2527 sage: B11 = matrix.identity(QQ,1)
2528 sage: B22 = M.transpose()*M
2529 sage: B = block_matrix(2,2,[ [B11,0 ],
2531 sage: J = BilinearFormEJA(B, orthonormalize=False)
2532 sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
2533 sage: V = J.vector_space()
2534 sage: sis = [ J( V([0] + (M.inverse()*ei).list()).column() )
2535 ....: for ei in eis ]
2536 sage: actual = [ sis[i]*sis[j]
2537 ....: for i in range(n-1)
2538 ....: for j in range(n-1) ]
2539 sage: expected = [ J.one() if i == j else J.zero()
2540 ....: for i in range(n-1)
2541 ....: for j in range(n-1) ]
2542 sage: actual == expected
2546 def __init__(self
, B
, **kwargs
):
2547 # The matrix "B" is supplied by the user in most cases,
2548 # so it makes sense to check whether or not its positive-
2549 # definite unless we are specifically asked not to...
2550 if ("check_axioms" not in kwargs
) or kwargs
["check_axioms"]:
2551 if not B
.is_positive_definite():
2552 raise ValueError("bilinear form is not positive-definite")
2554 # However, all of the other data for this EJA is computed
2555 # by us in manner that guarantees the axioms are
2556 # satisfied. So, again, unless we are specifically asked to
2557 # verify things, we'll skip the rest of the checks.
2558 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2560 def inner_product(x
,y
):
2561 return (y
.T
*B
*x
)[0,0]
2563 def jordan_product(x
,y
):
2569 z0
= inner_product(y
,x
)
2570 zbar
= y0
*xbar
+ x0
*ybar
2571 return P([z0
] + zbar
.list())
2574 column_basis
= tuple( b
.column() for b
in FreeModule(ZZ
, n
).basis() )
2575 super(BilinearFormEJA
, self
).__init
__(column_basis
,
2580 # The rank of this algebra is two, unless we're in a
2581 # one-dimensional ambient space (because the rank is bounded
2582 # by the ambient dimension).
2583 self
.rank
.set_cache(min(n
,2))
2586 self
.one
.set_cache( self
.zero() )
2588 self
.one
.set_cache( self
.monomial(0) )
2591 def _max_random_instance_size():
2593 The maximum dimension of a random BilinearFormEJA.
2598 def random_instance(cls
, **kwargs
):
2600 Return a random instance of this algebra.
2602 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2604 B
= matrix
.identity(ZZ
, n
)
2605 return cls(B
, **kwargs
)
2607 B11
= matrix
.identity(ZZ
, 1)
2608 M
= matrix
.random(ZZ
, n
-1)
2609 I
= matrix
.identity(ZZ
, n
-1)
2611 while alpha
.is_zero():
2612 alpha
= ZZ
.random_element().abs()
2613 B22
= M
.transpose()*M
+ alpha
*I
2615 from sage
.matrix
.special
import block_matrix
2616 B
= block_matrix(2,2, [ [B11
, ZZ(0) ],
2619 return cls(B
, **kwargs
)
2622 class JordanSpinEJA(BilinearFormEJA
):
2624 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2625 with the usual inner product and jordan product ``x*y =
2626 (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2631 sage: from mjo.eja.eja_algebra import JordanSpinEJA
2635 This multiplication table can be verified by hand::
2637 sage: J = JordanSpinEJA(4)
2638 sage: e0,e1,e2,e3 = J.gens()
2654 We can change the generator prefix::
2656 sage: JordanSpinEJA(2, prefix='B').gens()
2661 Ensure that we have the usual inner product on `R^n`::
2663 sage: set_random_seed()
2664 sage: J = JordanSpinEJA.random_instance()
2665 sage: x,y = J.random_elements(2)
2666 sage: actual = x.inner_product(y)
2667 sage: expected = x.to_vector().inner_product(y.to_vector())
2668 sage: actual == expected
2672 def __init__(self
, n
, **kwargs
):
2673 # This is a special case of the BilinearFormEJA with the
2674 # identity matrix as its bilinear form.
2675 B
= matrix
.identity(ZZ
, n
)
2677 # Don't orthonormalize because our basis is already
2678 # orthonormal with respect to our inner-product.
2679 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2681 # But also don't pass check_field=False here, because the user
2682 # can pass in a field!
2683 super(JordanSpinEJA
, self
).__init
__(B
, **kwargs
)
2686 def _max_random_instance_size():
2688 The maximum dimension of a random JordanSpinEJA.
2693 def random_instance(cls
, **kwargs
):
2695 Return a random instance of this type of algebra.
2697 Needed here to override the implementation for ``BilinearFormEJA``.
2699 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2700 return cls(n
, **kwargs
)
2703 class TrivialEJA(ConcreteEJA
):
2705 The trivial Euclidean Jordan algebra consisting of only a zero element.
2709 sage: from mjo.eja.eja_algebra import TrivialEJA
2713 sage: J = TrivialEJA()
2720 sage: 7*J.one()*12*J.one()
2722 sage: J.one().inner_product(J.one())
2724 sage: J.one().norm()
2726 sage: J.one().subalgebra_generated_by()
2727 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2732 def __init__(self
, **kwargs
):
2733 jordan_product
= lambda x
,y
: x
2734 inner_product
= lambda x
,y
: 0
2737 # New defaults for keyword arguments
2738 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2739 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2741 super(TrivialEJA
, self
).__init
__(basis
,
2745 # The rank is zero using my definition, namely the dimension of the
2746 # largest subalgebra generated by any element.
2747 self
.rank
.set_cache(0)
2748 self
.one
.set_cache( self
.zero() )
2751 def random_instance(cls
, **kwargs
):
2752 # We don't take a "size" argument so the superclass method is
2753 # inappropriate for us.
2754 return cls(**kwargs
)
2757 class CartesianProductEJA(CombinatorialFreeModule_CartesianProduct
,
2758 FiniteDimensionalEJA
):
2760 The external (orthogonal) direct sum of two or more Euclidean
2761 Jordan algebras. Every Euclidean Jordan algebra decomposes into an
2762 orthogonal direct sum of simple Euclidean Jordan algebras which is
2763 then isometric to a Cartesian product, so no generality is lost by
2764 providing only this construction.
2768 sage: from mjo.eja.eja_algebra import (random_eja,
2769 ....: CartesianProductEJA,
2771 ....: JordanSpinEJA,
2772 ....: RealSymmetricEJA)
2776 The Jordan product is inherited from our factors and implemented by
2777 our CombinatorialFreeModule Cartesian product superclass::
2779 sage: set_random_seed()
2780 sage: J1 = HadamardEJA(2)
2781 sage: J2 = RealSymmetricEJA(2)
2782 sage: J = cartesian_product([J1,J2])
2783 sage: x,y = J.random_elements(2)
2787 The ability to retrieve the original factors is implemented by our
2788 CombinatorialFreeModule Cartesian product superclass::
2790 sage: J1 = HadamardEJA(2, field=QQ)
2791 sage: J2 = JordanSpinEJA(3, field=QQ)
2792 sage: J = cartesian_product([J1,J2])
2793 sage: J.cartesian_factors()
2794 (Euclidean Jordan algebra of dimension 2 over Rational Field,
2795 Euclidean Jordan algebra of dimension 3 over Rational Field)
2797 You can provide more than two factors::
2799 sage: J1 = HadamardEJA(2)
2800 sage: J2 = JordanSpinEJA(3)
2801 sage: J3 = RealSymmetricEJA(3)
2802 sage: cartesian_product([J1,J2,J3])
2803 Euclidean Jordan algebra of dimension 2 over Algebraic Real
2804 Field (+) Euclidean Jordan algebra of dimension 3 over Algebraic
2805 Real Field (+) Euclidean Jordan algebra of dimension 6 over
2806 Algebraic Real Field
2808 Rank is additive on a Cartesian product::
2810 sage: J1 = HadamardEJA(1)
2811 sage: J2 = RealSymmetricEJA(2)
2812 sage: J = cartesian_product([J1,J2])
2813 sage: J1.rank.clear_cache()
2814 sage: J2.rank.clear_cache()
2815 sage: J.rank.clear_cache()
2818 sage: J.rank() == J1.rank() + J2.rank()
2821 The same rank computation works over the rationals, with whatever
2824 sage: J1 = HadamardEJA(1, field=QQ, orthonormalize=False)
2825 sage: J2 = RealSymmetricEJA(2, field=QQ, orthonormalize=False)
2826 sage: J = cartesian_product([J1,J2])
2827 sage: J1.rank.clear_cache()
2828 sage: J2.rank.clear_cache()
2829 sage: J.rank.clear_cache()
2832 sage: J.rank() == J1.rank() + J2.rank()
2835 The product algebra will be associative if and only if all of its
2836 components are associative::
2838 sage: J1 = HadamardEJA(2)
2839 sage: J1.is_associative()
2841 sage: J2 = HadamardEJA(3)
2842 sage: J2.is_associative()
2844 sage: J3 = RealSymmetricEJA(3)
2845 sage: J3.is_associative()
2847 sage: CP1 = cartesian_product([J1,J2])
2848 sage: CP1.is_associative()
2850 sage: CP2 = cartesian_product([J1,J3])
2851 sage: CP2.is_associative()
2856 All factors must share the same base field::
2858 sage: J1 = HadamardEJA(2, field=QQ)
2859 sage: J2 = RealSymmetricEJA(2)
2860 sage: CartesianProductEJA((J1,J2))
2861 Traceback (most recent call last):
2863 ValueError: all factors must share the same base field
2865 The cached unit element is the same one that would be computed::
2867 sage: set_random_seed() # long time
2868 sage: J1 = random_eja() # long time
2869 sage: J2 = random_eja() # long time
2870 sage: J = cartesian_product([J1,J2]) # long time
2871 sage: actual = J.one() # long time
2872 sage: J.one.clear_cache() # long time
2873 sage: expected = J.one() # long time
2874 sage: actual == expected # long time
2878 Element
= FiniteDimensionalEJAElement
2881 def __init__(self
, algebras
, **kwargs
):
2882 CombinatorialFreeModule_CartesianProduct
.__init
__(self
,
2885 field
= algebras
[0].base_ring()
2886 if not all( J
.base_ring() == field
for J
in algebras
):
2887 raise ValueError("all factors must share the same base field")
2889 associative
= all( m
.is_associative() for m
in algebras
)
2891 # The definition of matrix_space() and self.basis() relies
2892 # only on the stuff in the CFM_CartesianProduct class, which
2893 # we've already initialized.
2894 Js
= self
.cartesian_factors()
2896 MS
= self
.matrix_space()
2898 MS(tuple( self
.cartesian_projection(i
)(b
).to_matrix()
2899 for i
in range(m
) ))
2900 for b
in self
.basis()
2903 # Define jordan/inner products that operate on that matrix_basis.
2904 def jordan_product(x
,y
):
2906 (Js
[i
](x
[i
])*Js
[i
](y
[i
])).to_matrix() for i
in range(m
)
2909 def inner_product(x
, y
):
2911 Js
[i
](x
[i
]).inner_product(Js
[i
](y
[i
])) for i
in range(m
)
2914 # There's no need to check the field since it already came
2915 # from an EJA. Likewise the axioms are guaranteed to be
2916 # satisfied, unless the guy writing this class sucks.
2918 # If you want the basis to be orthonormalized, orthonormalize
2920 FiniteDimensionalEJA
.__init
__(self
,
2925 orthonormalize
=False,
2926 associative
=associative
,
2927 cartesian_product
=True,
2931 ones
= tuple(J
.one() for J
in algebras
)
2932 self
.one
.set_cache(self
._cartesian
_product
_of
_elements
(ones
))
2933 self
.rank
.set_cache(sum(J
.rank() for J
in algebras
))
2935 def matrix_space(self
):
2937 Return the space that our matrix basis lives in as a Cartesian
2942 sage: from mjo.eja.eja_algebra import (HadamardEJA,
2943 ....: RealSymmetricEJA)
2947 sage: J1 = HadamardEJA(1)
2948 sage: J2 = RealSymmetricEJA(2)
2949 sage: J = cartesian_product([J1,J2])
2950 sage: J.matrix_space()
2951 The Cartesian product of (Full MatrixSpace of 1 by 1 dense
2952 matrices over Algebraic Real Field, Full MatrixSpace of 2
2953 by 2 dense matrices over Algebraic Real Field)
2956 from sage
.categories
.cartesian_product
import cartesian_product
2957 return cartesian_product( [J
.matrix_space()
2958 for J
in self
.cartesian_factors()] )
2961 def cartesian_projection(self
, i
):
2965 sage: from mjo.eja.eja_algebra import (random_eja,
2966 ....: JordanSpinEJA,
2968 ....: RealSymmetricEJA,
2969 ....: ComplexHermitianEJA)
2973 The projection morphisms are Euclidean Jordan algebra
2976 sage: J1 = HadamardEJA(2)
2977 sage: J2 = RealSymmetricEJA(2)
2978 sage: J = cartesian_product([J1,J2])
2979 sage: J.cartesian_projection(0)
2980 Linear operator between finite-dimensional Euclidean Jordan
2981 algebras represented by the matrix:
2984 Domain: Euclidean Jordan algebra of dimension 2 over Algebraic
2985 Real Field (+) Euclidean Jordan algebra of dimension 3 over
2986 Algebraic Real Field
2987 Codomain: Euclidean Jordan algebra of dimension 2 over Algebraic
2989 sage: J.cartesian_projection(1)
2990 Linear operator between finite-dimensional Euclidean Jordan
2991 algebras represented by the matrix:
2995 Domain: Euclidean Jordan algebra of dimension 2 over Algebraic
2996 Real Field (+) Euclidean Jordan algebra of dimension 3 over
2997 Algebraic Real Field
2998 Codomain: Euclidean Jordan algebra of dimension 3 over Algebraic
3001 The projections work the way you'd expect on the vector
3002 representation of an element::
3004 sage: J1 = JordanSpinEJA(2)
3005 sage: J2 = ComplexHermitianEJA(2)
3006 sage: J = cartesian_product([J1,J2])
3007 sage: pi_left = J.cartesian_projection(0)
3008 sage: pi_right = J.cartesian_projection(1)
3009 sage: pi_left(J.one()).to_vector()
3011 sage: pi_right(J.one()).to_vector()
3013 sage: J.one().to_vector()
3018 The answer never changes::
3020 sage: set_random_seed()
3021 sage: J1 = random_eja()
3022 sage: J2 = random_eja()
3023 sage: J = cartesian_product([J1,J2])
3024 sage: P0 = J.cartesian_projection(0)
3025 sage: P1 = J.cartesian_projection(0)
3030 Ji
= self
.cartesian_factors()[i
]
3031 # Requires the fix on Trac 31421/31422 to work!
3032 Pi
= super().cartesian_projection(i
)
3033 return FiniteDimensionalEJAOperator(self
,Ji
,Pi
.matrix())
3036 def cartesian_embedding(self
, i
):
3040 sage: from mjo.eja.eja_algebra import (random_eja,
3041 ....: JordanSpinEJA,
3043 ....: RealSymmetricEJA)
3047 The embedding morphisms are Euclidean Jordan algebra
3050 sage: J1 = HadamardEJA(2)
3051 sage: J2 = RealSymmetricEJA(2)
3052 sage: J = cartesian_product([J1,J2])
3053 sage: J.cartesian_embedding(0)
3054 Linear operator between finite-dimensional Euclidean Jordan
3055 algebras represented by the matrix:
3061 Domain: Euclidean Jordan algebra of dimension 2 over
3062 Algebraic Real Field
3063 Codomain: Euclidean Jordan algebra of dimension 2 over
3064 Algebraic Real Field (+) Euclidean Jordan algebra of
3065 dimension 3 over Algebraic Real Field
3066 sage: J.cartesian_embedding(1)
3067 Linear operator between finite-dimensional Euclidean Jordan
3068 algebras represented by the matrix:
3074 Domain: Euclidean Jordan algebra of dimension 3 over
3075 Algebraic Real Field
3076 Codomain: Euclidean Jordan algebra of dimension 2 over
3077 Algebraic Real Field (+) Euclidean Jordan algebra of
3078 dimension 3 over Algebraic Real Field
3080 The embeddings work the way you'd expect on the vector
3081 representation of an element::
3083 sage: J1 = JordanSpinEJA(3)
3084 sage: J2 = RealSymmetricEJA(2)
3085 sage: J = cartesian_product([J1,J2])
3086 sage: iota_left = J.cartesian_embedding(0)
3087 sage: iota_right = J.cartesian_embedding(1)
3088 sage: iota_left(J1.zero()) == J.zero()
3090 sage: iota_right(J2.zero()) == J.zero()
3092 sage: J1.one().to_vector()
3094 sage: iota_left(J1.one()).to_vector()
3096 sage: J2.one().to_vector()
3098 sage: iota_right(J2.one()).to_vector()
3100 sage: J.one().to_vector()
3105 The answer never changes::
3107 sage: set_random_seed()
3108 sage: J1 = random_eja()
3109 sage: J2 = random_eja()
3110 sage: J = cartesian_product([J1,J2])
3111 sage: E0 = J.cartesian_embedding(0)
3112 sage: E1 = J.cartesian_embedding(0)
3116 Composing a projection with the corresponding inclusion should
3117 produce the identity map, and mismatching them should produce
3120 sage: set_random_seed()
3121 sage: J1 = random_eja()
3122 sage: J2 = random_eja()
3123 sage: J = cartesian_product([J1,J2])
3124 sage: iota_left = J.cartesian_embedding(0)
3125 sage: iota_right = J.cartesian_embedding(1)
3126 sage: pi_left = J.cartesian_projection(0)
3127 sage: pi_right = J.cartesian_projection(1)
3128 sage: pi_left*iota_left == J1.one().operator()
3130 sage: pi_right*iota_right == J2.one().operator()
3132 sage: (pi_left*iota_right).is_zero()
3134 sage: (pi_right*iota_left).is_zero()
3138 Ji
= self
.cartesian_factors()[i
]
3139 # Requires the fix on Trac 31421/31422 to work!
3140 Ei
= super().cartesian_embedding(i
)
3141 return FiniteDimensionalEJAOperator(Ji
,self
,Ei
.matrix())
3145 FiniteDimensionalEJA
.CartesianProduct
= CartesianProductEJA
3147 random_eja
= ConcreteEJA
.random_instance