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 "matrix form," which
46 must be the same form as the arguments to ``jordan_product``
47 and ``inner_product``. In reality, "matrix form" can be either
48 vectors, matrices, or a Cartesian product (ordered tuple)
49 of vectors or matrices. All of these would ideally be vector
50 spaces in sage with no special-casing needed; but in reality
51 we turn vectors into column-matrices and Cartesian products
52 `(a,b)` into column matrices `(a,b)^{T}` after converting
53 `a` and `b` themselves.
55 - jordan_product -- function of two elements (in matrix form)
56 that returns their jordan product in this algebra; this will
57 be applied to ``basis`` to compute a multiplication table for
60 - inner_product -- function of two elements (in matrix form) that
61 returns their inner product. This will be applied to ``basis`` to
62 compute an inner-product table (basically a matrix) for this algebra.
65 Element
= FiniteDimensionalEJAElement
74 cartesian_product
=False,
80 if not field
.is_subring(RR
):
81 # Note: this does return true for the real algebraic
82 # field, the rationals, and any quadratic field where
83 # we've specified a real embedding.
84 raise ValueError("scalar field is not real")
86 # If the basis given to us wasn't over the field that it's
87 # supposed to be over, fix that. Or, you know, crash.
88 if not cartesian_product
:
89 # The field for a cartesian product algebra comes from one
90 # of its factors and is the same for all factors, so
91 # there's no need to "reapply" it on product algebras.
92 basis
= tuple( b
.change_ring(field
) for b
in basis
)
96 # Check commutativity of the Jordan and inner-products.
97 # This has to be done before we build the multiplication
98 # and inner-product tables/matrices, because we take
99 # advantage of symmetry in the process.
100 if not all( jordan_product(bi
,bj
) == jordan_product(bj
,bi
)
103 raise ValueError("Jordan product is not commutative")
105 if not all( inner_product(bi
,bj
) == inner_product(bj
,bi
)
108 raise ValueError("inner-product is not commutative")
111 category
= MagmaticAlgebras(field
).FiniteDimensional()
112 category
= category
.WithBasis().Unital()
114 # Element subalgebras can take advantage of this.
115 category
= category
.Associative()
116 if cartesian_product
:
117 category
= category
.CartesianProducts()
119 # Call the superclass constructor so that we can use its from_vector()
120 # method to build our multiplication table.
122 super().__init
__(field
,
128 # Now comes all of the hard work. We'll be constructing an
129 # ambient vector space V that our (vectorized) basis lives in,
130 # as well as a subspace W of V spanned by those (vectorized)
131 # basis elements. The W-coordinates are the coefficients that
132 # we see in things like x = 1*e1 + 2*e2.
136 # flatten a vector, matrix, or cartesian product of those
137 # things into a long list.
138 if cartesian_product
:
139 return sum(( b_i
.list() for b_i
in b
), [])
145 degree
= len(flatten(basis
[0]))
147 # Build an ambient space that fits our matrix basis when
148 # written out as "long vectors."
149 V
= VectorSpace(field
, degree
)
151 # The matrix that will hole the orthonormal -> unorthonormal
152 # coordinate transformation.
153 self
._deortho
_matrix
= None
156 # Save a copy of the un-orthonormalized basis for later.
157 # Convert it to ambient V (vector) coordinates while we're
158 # at it, because we'd have to do it later anyway.
159 deortho_vector_basis
= tuple( V(flatten(b
)) for b
in basis
)
161 from mjo
.eja
.eja_utils
import gram_schmidt
162 basis
= tuple(gram_schmidt(basis
, inner_product
))
164 # Save the (possibly orthonormalized) matrix basis for
166 self
._matrix
_basis
= basis
168 # Now create the vector space for the algebra, which will have
169 # its own set of non-ambient coordinates (in terms of the
171 vector_basis
= tuple( V(flatten(b
)) for b
in basis
)
172 W
= V
.span_of_basis( vector_basis
, check
=check_axioms
)
175 # Now "W" is the vector space of our algebra coordinates. The
176 # variables "X1", "X2",... refer to the entries of vectors in
177 # W. Thus to convert back and forth between the orthonormal
178 # coordinates and the given ones, we need to stick the original
180 U
= V
.span_of_basis( deortho_vector_basis
, check
=check_axioms
)
181 self
._deortho
_matrix
= matrix( U
.coordinate_vector(q
)
182 for q
in vector_basis
)
185 # Now we actually compute the multiplication and inner-product
186 # tables/matrices using the possibly-orthonormalized basis.
187 self
._inner
_product
_matrix
= matrix
.identity(field
, n
)
188 self
._multiplication
_table
= [ [0 for j
in range(i
+1)]
191 # Note: the Jordan and inner-products are defined in terms
192 # of the ambient basis. It's important that their arguments
193 # are in ambient coordinates as well.
196 # ortho basis w.r.t. ambient coords
200 # The jordan product returns a matrixy answer, so we
201 # have to convert it to the algebra coordinates.
202 elt
= jordan_product(q_i
, q_j
)
203 elt
= W
.coordinate_vector(V(flatten(elt
)))
204 self
._multiplication
_table
[i
][j
] = self
.from_vector(elt
)
206 if not orthonormalize
:
207 # If we're orthonormalizing the basis with respect
208 # to an inner-product, then the inner-product
209 # matrix with respect to the resulting basis is
210 # just going to be the identity.
211 ip
= inner_product(q_i
, q_j
)
212 self
._inner
_product
_matrix
[i
,j
] = ip
213 self
._inner
_product
_matrix
[j
,i
] = ip
215 self
._inner
_product
_matrix
._cache
= {'hermitian': True}
216 self
._inner
_product
_matrix
.set_immutable()
219 if not self
._is
_jordanian
():
220 raise ValueError("Jordan identity does not hold")
221 if not self
._inner
_product
_is
_associative
():
222 raise ValueError("inner product is not associative")
225 def _coerce_map_from_base_ring(self
):
227 Disable the map from the base ring into the algebra.
229 Performing a nonsense conversion like this automatically
230 is counterpedagogical. The fallback is to try the usual
231 element constructor, which should also fail.
235 sage: from mjo.eja.eja_algebra import random_eja
239 sage: set_random_seed()
240 sage: J = random_eja()
242 Traceback (most recent call last):
244 ValueError: not an element of this algebra
250 def product_on_basis(self
, i
, j
):
251 # We only stored the lower-triangular portion of the
252 # multiplication table.
254 return self
._multiplication
_table
[i
][j
]
256 return self
._multiplication
_table
[j
][i
]
258 def inner_product(self
, x
, y
):
260 The inner product associated with this Euclidean Jordan algebra.
262 Defaults to the trace inner product, but can be overridden by
263 subclasses if they are sure that the necessary properties are
268 sage: from mjo.eja.eja_algebra import (random_eja,
270 ....: BilinearFormEJA)
274 Our inner product is "associative," which means the following for
275 a symmetric bilinear form::
277 sage: set_random_seed()
278 sage: J = random_eja()
279 sage: x,y,z = J.random_elements(3)
280 sage: (x*y).inner_product(z) == y.inner_product(x*z)
285 Ensure that this is the usual inner product for the algebras
288 sage: set_random_seed()
289 sage: J = HadamardEJA.random_instance()
290 sage: x,y = J.random_elements(2)
291 sage: actual = x.inner_product(y)
292 sage: expected = x.to_vector().inner_product(y.to_vector())
293 sage: actual == expected
296 Ensure that this is one-half of the trace inner-product in a
297 BilinearFormEJA that isn't just the reals (when ``n`` isn't
298 one). This is in Faraut and Koranyi, and also my "On the
301 sage: set_random_seed()
302 sage: J = BilinearFormEJA.random_instance()
303 sage: n = J.dimension()
304 sage: x = J.random_element()
305 sage: y = J.random_element()
306 sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
310 B
= self
._inner
_product
_matrix
311 return (B
*x
.to_vector()).inner_product(y
.to_vector())
314 def is_associative(self
):
316 Return whether or not this algebra's Jordan product is associative.
320 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
324 sage: J = ComplexHermitianEJA(3, field=QQ, orthonormalize=False)
325 sage: J.is_associative()
327 sage: x = sum(J.gens())
328 sage: A = x.subalgebra_generated_by(orthonormalize=False)
329 sage: A.is_associative()
333 return "Associative" in self
.category().axioms()
335 def _is_jordanian(self
):
337 Whether or not this algebra's multiplication table respects the
338 Jordan identity `(x^{2})(xy) = x(x^{2}y)`.
340 We only check one arrangement of `x` and `y`, so for a
341 ``True`` result to be truly true, you should also check
342 :meth:`is_commutative`. This method should of course always
343 return ``True``, unless this algebra was constructed with
344 ``check_axioms=False`` and passed an invalid multiplication table.
346 return all( (self
.gens()[i
]**2)*(self
.gens()[i
]*self
.gens()[j
])
348 (self
.gens()[i
])*((self
.gens()[i
]**2)*self
.gens()[j
])
349 for i
in range(self
.dimension())
350 for j
in range(self
.dimension()) )
352 def _inner_product_is_associative(self
):
354 Return whether or not this algebra's inner product `B` is
355 associative; that is, whether or not `B(xy,z) = B(x,yz)`.
357 This method should of course always return ``True``, unless
358 this algebra was constructed with ``check_axioms=False`` and
359 passed an invalid Jordan or inner-product.
362 # Used to check whether or not something is zero in an inexact
363 # ring. This number is sufficient to allow the construction of
364 # QuaternionHermitianEJA(2, field=RDF) with check_axioms=True.
367 for i
in range(self
.dimension()):
368 for j
in range(self
.dimension()):
369 for k
in range(self
.dimension()):
373 diff
= (x
*y
).inner_product(z
) - x
.inner_product(y
*z
)
375 if self
.base_ring().is_exact():
379 if diff
.abs() > epsilon
:
384 def _element_constructor_(self
, elt
):
386 Construct an element of this algebra from its vector or matrix
389 This gets called only after the parent element _call_ method
390 fails to find a coercion for the argument.
394 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
396 ....: RealSymmetricEJA)
400 The identity in `S^n` is converted to the identity in the EJA::
402 sage: J = RealSymmetricEJA(3)
403 sage: I = matrix.identity(QQ,3)
404 sage: J(I) == J.one()
407 This skew-symmetric matrix can't be represented in the EJA::
409 sage: J = RealSymmetricEJA(3)
410 sage: A = matrix(QQ,3, lambda i,j: i-j)
412 Traceback (most recent call last):
414 ValueError: not an element of this algebra
418 Ensure that we can convert any element of the two non-matrix
419 simple algebras (whose matrix representations are columns)
420 back and forth faithfully::
422 sage: set_random_seed()
423 sage: J = HadamardEJA.random_instance()
424 sage: x = J.random_element()
425 sage: J(x.to_vector().column()) == x
427 sage: J = JordanSpinEJA.random_instance()
428 sage: x = J.random_element()
429 sage: J(x.to_vector().column()) == x
433 msg
= "not an element of this algebra"
435 # The superclass implementation of random_element()
436 # needs to be able to coerce "0" into the algebra.
438 elif elt
in self
.base_ring():
439 # Ensure that no base ring -> algebra coercion is performed
440 # by this method. There's some stupidity in sage that would
441 # otherwise propagate to this method; for example, sage thinks
442 # that the integer 3 belongs to the space of 2-by-2 matrices.
443 raise ValueError(msg
)
447 except (AttributeError, TypeError):
448 # Try to convert a vector into a column-matrix
451 if elt
not in self
.matrix_space():
452 raise ValueError(msg
)
454 # Thanks for nothing! Matrix spaces aren't vector spaces in
455 # Sage, so we have to figure out its matrix-basis coordinates
456 # ourselves. We use the basis space's ring instead of the
457 # element's ring because the basis space might be an algebraic
458 # closure whereas the base ring of the 3-by-3 identity matrix
459 # could be QQ instead of QQbar.
461 # We pass check=False because the matrix basis is "guaranteed"
462 # to be linearly independent... right? Ha ha.
463 V
= VectorSpace(self
.base_ring(), elt
.nrows()*elt
.ncols())
464 W
= V
.span_of_basis( (_mat2vec(s
) for s
in self
.matrix_basis()),
468 coords
= W
.coordinate_vector(_mat2vec(elt
))
469 except ArithmeticError: # vector is not in free module
470 raise ValueError(msg
)
472 return self
.from_vector(coords
)
476 Return a string representation of ``self``.
480 sage: from mjo.eja.eja_algebra import JordanSpinEJA
484 Ensure that it says what we think it says::
486 sage: JordanSpinEJA(2, field=AA)
487 Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
488 sage: JordanSpinEJA(3, field=RDF)
489 Euclidean Jordan algebra of dimension 3 over Real Double Field
492 fmt
= "Euclidean Jordan algebra of dimension {} over {}"
493 return fmt
.format(self
.dimension(), self
.base_ring())
497 def characteristic_polynomial_of(self
):
499 Return the algebra's "characteristic polynomial of" function,
500 which is itself a multivariate polynomial that, when evaluated
501 at the coordinates of some algebra element, returns that
502 element's characteristic polynomial.
504 The resulting polynomial has `n+1` variables, where `n` is the
505 dimension of this algebra. The first `n` variables correspond to
506 the coordinates of an algebra element: when evaluated at the
507 coordinates of an algebra element with respect to a certain
508 basis, the result is a univariate polynomial (in the one
509 remaining variable ``t``), namely the characteristic polynomial
514 sage: from mjo.eja.eja_algebra import JordanSpinEJA, TrivialEJA
518 The characteristic polynomial in the spin algebra is given in
519 Alizadeh, Example 11.11::
521 sage: J = JordanSpinEJA(3)
522 sage: p = J.characteristic_polynomial_of(); p
523 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
524 sage: xvec = J.one().to_vector()
528 By definition, the characteristic polynomial is a monic
529 degree-zero polynomial in a rank-zero algebra. Note that
530 Cayley-Hamilton is indeed satisfied since the polynomial
531 ``1`` evaluates to the identity element of the algebra on
534 sage: J = TrivialEJA()
535 sage: J.characteristic_polynomial_of()
542 # The list of coefficient polynomials a_0, a_1, a_2, ..., a_(r-1).
543 a
= self
._charpoly
_coefficients
()
545 # We go to a bit of trouble here to reorder the
546 # indeterminates, so that it's easier to evaluate the
547 # characteristic polynomial at x's coordinates and get back
548 # something in terms of t, which is what we want.
549 S
= PolynomialRing(self
.base_ring(),'t')
553 S
= PolynomialRing(S
, R
.variable_names())
556 return (t
**r
+ sum( a
[k
]*(t
**k
) for k
in range(r
) ))
558 def coordinate_polynomial_ring(self
):
560 The multivariate polynomial ring in which this algebra's
561 :meth:`characteristic_polynomial_of` lives.
565 sage: from mjo.eja.eja_algebra import (HadamardEJA,
566 ....: RealSymmetricEJA)
570 sage: J = HadamardEJA(2)
571 sage: J.coordinate_polynomial_ring()
572 Multivariate Polynomial Ring in X1, X2...
573 sage: J = RealSymmetricEJA(3,field=QQ,orthonormalize=False)
574 sage: J.coordinate_polynomial_ring()
575 Multivariate Polynomial Ring in X1, X2, X3, X4, X5, X6...
578 var_names
= tuple( "X%d" % z
for z
in range(1, self
.dimension()+1) )
579 return PolynomialRing(self
.base_ring(), var_names
)
581 def inner_product(self
, x
, y
):
583 The inner product associated with this Euclidean Jordan algebra.
585 Defaults to the trace inner product, but can be overridden by
586 subclasses if they are sure that the necessary properties are
591 sage: from mjo.eja.eja_algebra import (random_eja,
593 ....: BilinearFormEJA)
597 Our inner product is "associative," which means the following for
598 a symmetric bilinear form::
600 sage: set_random_seed()
601 sage: J = random_eja()
602 sage: x,y,z = J.random_elements(3)
603 sage: (x*y).inner_product(z) == y.inner_product(x*z)
608 Ensure that this is the usual inner product for the algebras
611 sage: set_random_seed()
612 sage: J = HadamardEJA.random_instance()
613 sage: x,y = J.random_elements(2)
614 sage: actual = x.inner_product(y)
615 sage: expected = x.to_vector().inner_product(y.to_vector())
616 sage: actual == expected
619 Ensure that this is one-half of the trace inner-product in a
620 BilinearFormEJA that isn't just the reals (when ``n`` isn't
621 one). This is in Faraut and Koranyi, and also my "On the
624 sage: set_random_seed()
625 sage: J = BilinearFormEJA.random_instance()
626 sage: n = J.dimension()
627 sage: x = J.random_element()
628 sage: y = J.random_element()
629 sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
632 B
= self
._inner
_product
_matrix
633 return (B
*x
.to_vector()).inner_product(y
.to_vector())
636 def is_trivial(self
):
638 Return whether or not this algebra is trivial.
640 A trivial algebra contains only the zero element.
644 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
649 sage: J = ComplexHermitianEJA(3)
655 sage: J = TrivialEJA()
660 return self
.dimension() == 0
663 def multiplication_table(self
):
665 Return a visual representation of this algebra's multiplication
666 table (on basis elements).
670 sage: from mjo.eja.eja_algebra import JordanSpinEJA
674 sage: J = JordanSpinEJA(4)
675 sage: J.multiplication_table()
676 +----++----+----+----+----+
677 | * || e0 | e1 | e2 | e3 |
678 +====++====+====+====+====+
679 | e0 || e0 | e1 | e2 | e3 |
680 +----++----+----+----+----+
681 | e1 || e1 | e0 | 0 | 0 |
682 +----++----+----+----+----+
683 | e2 || e2 | 0 | e0 | 0 |
684 +----++----+----+----+----+
685 | e3 || e3 | 0 | 0 | e0 |
686 +----++----+----+----+----+
690 # Prepend the header row.
691 M
= [["*"] + list(self
.gens())]
693 # And to each subsequent row, prepend an entry that belongs to
694 # the left-side "header column."
695 M
+= [ [self
.gens()[i
]] + [ self
.product_on_basis(i
,j
)
699 return table(M
, header_row
=True, header_column
=True, frame
=True)
702 def matrix_basis(self
):
704 Return an (often more natural) representation of this algebras
705 basis as an ordered tuple of matrices.
707 Every finite-dimensional Euclidean Jordan Algebra is a, up to
708 Jordan isomorphism, a direct sum of five simple
709 algebras---four of which comprise Hermitian matrices. And the
710 last type of algebra can of course be thought of as `n`-by-`1`
711 column matrices (ambiguusly called column vectors) to avoid
712 special cases. As a result, matrices (and column vectors) are
713 a natural representation format for Euclidean Jordan algebra
716 But, when we construct an algebra from a basis of matrices,
717 those matrix representations are lost in favor of coordinate
718 vectors *with respect to* that basis. We could eventually
719 convert back if we tried hard enough, but having the original
720 representations handy is valuable enough that we simply store
721 them and return them from this method.
723 Why implement this for non-matrix algebras? Avoiding special
724 cases for the :class:`BilinearFormEJA` pays with simplicity in
725 its own right. But mainly, we would like to be able to assume
726 that elements of a :class:`CartesianProductEJA` can be displayed
727 nicely, without having to have special classes for direct sums
728 one of whose components was a matrix algebra.
732 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
733 ....: RealSymmetricEJA)
737 sage: J = RealSymmetricEJA(2)
739 Finite family {0: e0, 1: e1, 2: e2}
740 sage: J.matrix_basis()
742 [1 0] [ 0 0.7071067811865475?] [0 0]
743 [0 0], [0.7071067811865475? 0], [0 1]
748 sage: J = JordanSpinEJA(2)
750 Finite family {0: e0, 1: e1}
751 sage: J.matrix_basis()
757 return self
._matrix
_basis
760 def matrix_space(self
):
762 Return the matrix space in which this algebra's elements live, if
763 we think of them as matrices (including column vectors of the
766 Generally this will be an `n`-by-`1` column-vector space,
767 except when the algebra is trivial. There it's `n`-by-`n`
768 (where `n` is zero), to ensure that two elements of the matrix
769 space (empty matrices) can be multiplied.
771 Matrix algebras override this with something more useful.
773 if self
.is_trivial():
774 return MatrixSpace(self
.base_ring(), 0)
776 return self
.matrix_basis()[0].parent()
782 Return the unit element of this algebra.
786 sage: from mjo.eja.eja_algebra import (HadamardEJA,
791 We can compute unit element in the Hadamard EJA::
793 sage: J = HadamardEJA(5)
795 e0 + e1 + e2 + e3 + e4
797 The unit element in the Hadamard EJA is inherited in the
798 subalgebras generated by its elements::
800 sage: J = HadamardEJA(5)
802 e0 + e1 + e2 + e3 + e4
803 sage: x = sum(J.gens())
804 sage: A = x.subalgebra_generated_by(orthonormalize=False)
807 sage: A.one().superalgebra_element()
808 e0 + e1 + e2 + e3 + e4
812 The identity element acts like the identity, regardless of
813 whether or not we orthonormalize::
815 sage: set_random_seed()
816 sage: J = random_eja()
817 sage: x = J.random_element()
818 sage: J.one()*x == x and x*J.one() == x
820 sage: A = x.subalgebra_generated_by()
821 sage: y = A.random_element()
822 sage: A.one()*y == y and y*A.one() == y
827 sage: set_random_seed()
828 sage: J = random_eja(field=QQ, orthonormalize=False)
829 sage: x = J.random_element()
830 sage: J.one()*x == x and x*J.one() == x
832 sage: A = x.subalgebra_generated_by(orthonormalize=False)
833 sage: y = A.random_element()
834 sage: A.one()*y == y and y*A.one() == y
837 The matrix of the unit element's operator is the identity,
838 regardless of the base field and whether or not we
841 sage: set_random_seed()
842 sage: J = random_eja()
843 sage: actual = J.one().operator().matrix()
844 sage: expected = matrix.identity(J.base_ring(), J.dimension())
845 sage: actual == expected
847 sage: x = J.random_element()
848 sage: A = x.subalgebra_generated_by()
849 sage: actual = A.one().operator().matrix()
850 sage: expected = matrix.identity(A.base_ring(), A.dimension())
851 sage: actual == expected
856 sage: set_random_seed()
857 sage: J = random_eja(field=QQ, orthonormalize=False)
858 sage: actual = J.one().operator().matrix()
859 sage: expected = matrix.identity(J.base_ring(), J.dimension())
860 sage: actual == expected
862 sage: x = J.random_element()
863 sage: A = x.subalgebra_generated_by(orthonormalize=False)
864 sage: actual = A.one().operator().matrix()
865 sage: expected = matrix.identity(A.base_ring(), A.dimension())
866 sage: actual == expected
869 Ensure that the cached unit element (often precomputed by
870 hand) agrees with the computed one::
872 sage: set_random_seed()
873 sage: J = random_eja()
874 sage: cached = J.one()
875 sage: J.one.clear_cache()
876 sage: J.one() == cached
881 sage: set_random_seed()
882 sage: J = random_eja(field=QQ, orthonormalize=False)
883 sage: cached = J.one()
884 sage: J.one.clear_cache()
885 sage: J.one() == cached
889 # We can brute-force compute the matrices of the operators
890 # that correspond to the basis elements of this algebra.
891 # If some linear combination of those basis elements is the
892 # algebra identity, then the same linear combination of
893 # their matrices has to be the identity matrix.
895 # Of course, matrices aren't vectors in sage, so we have to
896 # appeal to the "long vectors" isometry.
897 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
899 # Now we use basic linear algebra to find the coefficients,
900 # of the matrices-as-vectors-linear-combination, which should
901 # work for the original algebra basis too.
902 A
= matrix(self
.base_ring(), oper_vecs
)
904 # We used the isometry on the left-hand side already, but we
905 # still need to do it for the right-hand side. Recall that we
906 # wanted something that summed to the identity matrix.
907 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
909 # Now if there's an identity element in the algebra, this
910 # should work. We solve on the left to avoid having to
911 # transpose the matrix "A".
912 return self
.from_vector(A
.solve_left(b
))
915 def peirce_decomposition(self
, c
):
917 The Peirce decomposition of this algebra relative to the
920 In the future, this can be extended to a complete system of
921 orthogonal idempotents.
925 - ``c`` -- an idempotent of this algebra.
929 A triple (J0, J5, J1) containing two subalgebras and one subspace
932 - ``J0`` -- the algebra on the eigenspace of ``c.operator()``
933 corresponding to the eigenvalue zero.
935 - ``J5`` -- the eigenspace (NOT a subalgebra) of ``c.operator()``
936 corresponding to the eigenvalue one-half.
938 - ``J1`` -- the algebra on the eigenspace of ``c.operator()``
939 corresponding to the eigenvalue one.
941 These are the only possible eigenspaces for that operator, and this
942 algebra is a direct sum of them. The spaces ``J0`` and ``J1`` are
943 orthogonal, and are subalgebras of this algebra with the appropriate
948 sage: from mjo.eja.eja_algebra import random_eja, RealSymmetricEJA
952 The canonical example comes from the symmetric matrices, which
953 decompose into diagonal and off-diagonal parts::
955 sage: J = RealSymmetricEJA(3)
956 sage: C = matrix(QQ, [ [1,0,0],
960 sage: J0,J5,J1 = J.peirce_decomposition(c)
962 Euclidean Jordan algebra of dimension 1...
964 Vector space of degree 6 and dimension 2...
966 Euclidean Jordan algebra of dimension 3...
967 sage: J0.one().to_matrix()
971 sage: orig_df = AA.options.display_format
972 sage: AA.options.display_format = 'radical'
973 sage: J.from_vector(J5.basis()[0]).to_matrix()
977 sage: J.from_vector(J5.basis()[1]).to_matrix()
981 sage: AA.options.display_format = orig_df
982 sage: J1.one().to_matrix()
989 Every algebra decomposes trivially with respect to its identity
992 sage: set_random_seed()
993 sage: J = random_eja()
994 sage: J0,J5,J1 = J.peirce_decomposition(J.one())
995 sage: J0.dimension() == 0 and J5.dimension() == 0
997 sage: J1.superalgebra() == J and J1.dimension() == J.dimension()
1000 The decomposition is into eigenspaces, and its components are
1001 therefore necessarily orthogonal. Moreover, the identity
1002 elements in the two subalgebras are the projections onto their
1003 respective subspaces of the superalgebra's identity element::
1005 sage: set_random_seed()
1006 sage: J = random_eja()
1007 sage: x = J.random_element()
1008 sage: if not J.is_trivial():
1009 ....: while x.is_nilpotent():
1010 ....: x = J.random_element()
1011 sage: c = x.subalgebra_idempotent()
1012 sage: J0,J5,J1 = J.peirce_decomposition(c)
1014 sage: for (w,y,z) in zip(J0.basis(), J5.basis(), J1.basis()):
1015 ....: w = w.superalgebra_element()
1016 ....: y = J.from_vector(y)
1017 ....: z = z.superalgebra_element()
1018 ....: ipsum += w.inner_product(y).abs()
1019 ....: ipsum += w.inner_product(z).abs()
1020 ....: ipsum += y.inner_product(z).abs()
1023 sage: J1(c) == J1.one()
1025 sage: J0(J.one() - c) == J0.one()
1029 if not c
.is_idempotent():
1030 raise ValueError("element is not idempotent: %s" % c
)
1032 from mjo
.eja
.eja_subalgebra
import FiniteDimensionalEJASubalgebra
1034 # Default these to what they should be if they turn out to be
1035 # trivial, because eigenspaces_left() won't return eigenvalues
1036 # corresponding to trivial spaces (e.g. it returns only the
1037 # eigenspace corresponding to lambda=1 if you take the
1038 # decomposition relative to the identity element).
1039 trivial
= FiniteDimensionalEJASubalgebra(self
, ())
1040 J0
= trivial
# eigenvalue zero
1041 J5
= VectorSpace(self
.base_ring(), 0) # eigenvalue one-half
1042 J1
= trivial
# eigenvalue one
1044 for (eigval
, eigspace
) in c
.operator().matrix().right_eigenspaces():
1045 if eigval
== ~
(self
.base_ring()(2)):
1048 gens
= tuple( self
.from_vector(b
) for b
in eigspace
.basis() )
1049 subalg
= FiniteDimensionalEJASubalgebra(self
,
1057 raise ValueError("unexpected eigenvalue: %s" % eigval
)
1062 def random_element(self
, thorough
=False):
1064 Return a random element of this algebra.
1066 Our algebra superclass method only returns a linear
1067 combination of at most two basis elements. We instead
1068 want the vector space "random element" method that
1069 returns a more diverse selection.
1073 - ``thorough`` -- (boolean; default False) whether or not we
1074 should generate irrational coefficients for the random
1075 element when our base ring is irrational; this slows the
1076 algebra operations to a crawl, but any truly random method
1080 # For a general base ring... maybe we can trust this to do the
1081 # right thing? Unlikely, but.
1082 V
= self
.vector_space()
1083 v
= V
.random_element()
1085 if self
.base_ring() is AA
:
1086 # The "random element" method of the algebraic reals is
1087 # stupid at the moment, and only returns integers between
1088 # -2 and 2, inclusive:
1090 # https://trac.sagemath.org/ticket/30875
1092 # Instead, we implement our own "random vector" method,
1093 # and then coerce that into the algebra. We use the vector
1094 # space degree here instead of the dimension because a
1095 # subalgebra could (for example) be spanned by only two
1096 # vectors, each with five coordinates. We need to
1097 # generate all five coordinates.
1099 v
*= QQbar
.random_element().real()
1101 v
*= QQ
.random_element()
1103 return self
.from_vector(V
.coordinate_vector(v
))
1105 def random_elements(self
, count
, thorough
=False):
1107 Return ``count`` random elements as a tuple.
1111 - ``thorough`` -- (boolean; default False) whether or not we
1112 should generate irrational coefficients for the random
1113 elements when our base ring is irrational; this slows the
1114 algebra operations to a crawl, but any truly random method
1119 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1123 sage: J = JordanSpinEJA(3)
1124 sage: x,y,z = J.random_elements(3)
1125 sage: all( [ x in J, y in J, z in J ])
1127 sage: len( J.random_elements(10) ) == 10
1131 return tuple( self
.random_element(thorough
)
1132 for idx
in range(count
) )
1136 def _charpoly_coefficients(self
):
1138 The `r` polynomial coefficients of the "characteristic polynomial
1143 sage: from mjo.eja.eja_algebra import random_eja
1147 The theory shows that these are all homogeneous polynomials of
1150 sage: set_random_seed()
1151 sage: J = random_eja()
1152 sage: all(p.is_homogeneous() for p in J._charpoly_coefficients())
1156 n
= self
.dimension()
1157 R
= self
.coordinate_polynomial_ring()
1159 F
= R
.fraction_field()
1162 # From a result in my book, these are the entries of the
1163 # basis representation of L_x.
1164 return sum( vars[k
]*self
.gens()[k
].operator().matrix()[i
,j
]
1167 L_x
= matrix(F
, n
, n
, L_x_i_j
)
1170 if self
.rank
.is_in_cache():
1172 # There's no need to pad the system with redundant
1173 # columns if we *know* they'll be redundant.
1176 # Compute an extra power in case the rank is equal to
1177 # the dimension (otherwise, we would stop at x^(r-1)).
1178 x_powers
= [ (L_x
**k
)*self
.one().to_vector()
1179 for k
in range(n
+1) ]
1180 A
= matrix
.column(F
, x_powers
[:n
])
1181 AE
= A
.extended_echelon_form()
1188 # The theory says that only the first "r" coefficients are
1189 # nonzero, and they actually live in the original polynomial
1190 # ring and not the fraction field. We negate them because in
1191 # the actual characteristic polynomial, they get moved to the
1192 # other side where x^r lives. We don't bother to trim A_rref
1193 # down to a square matrix and solve the resulting system,
1194 # because the upper-left r-by-r portion of A_rref is
1195 # guaranteed to be the identity matrix, so e.g.
1197 # A_rref.solve_right(Y)
1199 # would just be returning Y.
1200 return (-E
*b
)[:r
].change_ring(R
)
1205 Return the rank of this EJA.
1207 This is a cached method because we know the rank a priori for
1208 all of the algebras we can construct. Thus we can avoid the
1209 expensive ``_charpoly_coefficients()`` call unless we truly
1210 need to compute the whole characteristic polynomial.
1214 sage: from mjo.eja.eja_algebra import (HadamardEJA,
1215 ....: JordanSpinEJA,
1216 ....: RealSymmetricEJA,
1217 ....: ComplexHermitianEJA,
1218 ....: QuaternionHermitianEJA,
1223 The rank of the Jordan spin algebra is always two::
1225 sage: JordanSpinEJA(2).rank()
1227 sage: JordanSpinEJA(3).rank()
1229 sage: JordanSpinEJA(4).rank()
1232 The rank of the `n`-by-`n` Hermitian real, complex, or
1233 quaternion matrices is `n`::
1235 sage: RealSymmetricEJA(4).rank()
1237 sage: ComplexHermitianEJA(3).rank()
1239 sage: QuaternionHermitianEJA(2).rank()
1244 Ensure that every EJA that we know how to construct has a
1245 positive integer rank, unless the algebra is trivial in
1246 which case its rank will be zero::
1248 sage: set_random_seed()
1249 sage: J = random_eja()
1253 sage: r > 0 or (r == 0 and J.is_trivial())
1256 Ensure that computing the rank actually works, since the ranks
1257 of all simple algebras are known and will be cached by default::
1259 sage: set_random_seed() # long time
1260 sage: J = random_eja() # long time
1261 sage: cached = J.rank() # long time
1262 sage: J.rank.clear_cache() # long time
1263 sage: J.rank() == cached # long time
1267 return len(self
._charpoly
_coefficients
())
1270 def vector_space(self
):
1272 Return the vector space that underlies this algebra.
1276 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1280 sage: J = RealSymmetricEJA(2)
1281 sage: J.vector_space()
1282 Vector space of dimension 3 over...
1285 return self
.zero().to_vector().parent().ambient_vector_space()
1288 Element
= FiniteDimensionalEJAElement
1290 class RationalBasisEJA(FiniteDimensionalEJA
):
1292 New class for algebras whose supplied basis elements have all rational entries.
1296 sage: from mjo.eja.eja_algebra import BilinearFormEJA
1300 The supplied basis is orthonormalized by default::
1302 sage: B = matrix(QQ, [[1, 0, 0], [0, 25, -32], [0, -32, 41]])
1303 sage: J = BilinearFormEJA(B)
1304 sage: J.matrix_basis()
1321 # Abuse the check_field parameter to check that the entries of
1322 # out basis (in ambient coordinates) are in the field QQ.
1323 if not all( all(b_i
in QQ
for b_i
in b
.list()) for b
in basis
):
1324 raise TypeError("basis not rational")
1326 self
._rational
_algebra
= None
1328 # There's no point in constructing the extra algebra if this
1329 # one is already rational.
1331 # Note: the same Jordan and inner-products work here,
1332 # because they are necessarily defined with respect to
1333 # ambient coordinates and not any particular basis.
1334 self
._rational
_algebra
= FiniteDimensionalEJA(
1339 orthonormalize
=False,
1343 super().__init
__(basis
,
1347 check_field
=check_field
,
1351 def _charpoly_coefficients(self
):
1355 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
1356 ....: JordanSpinEJA)
1360 The base ring of the resulting polynomial coefficients is what
1361 it should be, and not the rationals (unless the algebra was
1362 already over the rationals)::
1364 sage: J = JordanSpinEJA(3)
1365 sage: J._charpoly_coefficients()
1366 (X1^2 - X2^2 - X3^2, -2*X1)
1367 sage: a0 = J._charpoly_coefficients()[0]
1369 Algebraic Real Field
1370 sage: a0.base_ring()
1371 Algebraic Real Field
1374 if self
._rational
_algebra
is None:
1375 # There's no need to construct *another* algebra over the
1376 # rationals if this one is already over the
1377 # rationals. Likewise, if we never orthonormalized our
1378 # basis, we might as well just use the given one.
1379 return super()._charpoly
_coefficients
()
1381 # Do the computation over the rationals. The answer will be
1382 # the same, because all we've done is a change of basis.
1383 # Then, change back from QQ to our real base ring
1384 a
= ( a_i
.change_ring(self
.base_ring())
1385 for a_i
in self
._rational
_algebra
._charpoly
_coefficients
() )
1387 if self
._deortho
_matrix
is None:
1388 # This can happen if our base ring was, say, AA and we
1389 # chose not to (or didn't need to) orthonormalize. It's
1390 # still faster to do the computations over QQ even if
1391 # the numbers in the boxes stay the same.
1394 # Otherwise, convert the coordinate variables back to the
1395 # deorthonormalized ones.
1396 R
= self
.coordinate_polynomial_ring()
1397 from sage
.modules
.free_module_element
import vector
1398 X
= vector(R
, R
.gens())
1399 BX
= self
._deortho
_matrix
*X
1401 subs_dict
= { X[i]: BX[i] for i in range(len(X)) }
1402 return tuple( a_i
.subs(subs_dict
) for a_i
in a
)
1404 class ConcreteEJA(RationalBasisEJA
):
1406 A class for the Euclidean Jordan algebras that we know by name.
1408 These are the Jordan algebras whose basis, multiplication table,
1409 rank, and so on are known a priori. More to the point, they are
1410 the Euclidean Jordan algebras for which we are able to conjure up
1411 a "random instance."
1415 sage: from mjo.eja.eja_algebra import ConcreteEJA
1419 Our basis is normalized with respect to the algebra's inner
1420 product, unless we specify otherwise::
1422 sage: set_random_seed()
1423 sage: J = ConcreteEJA.random_instance()
1424 sage: all( b.norm() == 1 for b in J.gens() )
1427 Since our basis is orthonormal with respect to the algebra's inner
1428 product, and since we know that this algebra is an EJA, any
1429 left-multiplication operator's matrix will be symmetric because
1430 natural->EJA basis representation is an isometry and within the
1431 EJA the operator is self-adjoint by the Jordan axiom::
1433 sage: set_random_seed()
1434 sage: J = ConcreteEJA.random_instance()
1435 sage: x = J.random_element()
1436 sage: x.operator().is_self_adjoint()
1441 def _max_random_instance_size():
1443 Return an integer "size" that is an upper bound on the size of
1444 this algebra when it is used in a random test
1445 case. Unfortunately, the term "size" is ambiguous -- when
1446 dealing with `R^n` under either the Hadamard or Jordan spin
1447 product, the "size" refers to the dimension `n`. When dealing
1448 with a matrix algebra (real symmetric or complex/quaternion
1449 Hermitian), it refers to the size of the matrix, which is far
1450 less than the dimension of the underlying vector space.
1452 This method must be implemented in each subclass.
1454 raise NotImplementedError
1457 def random_instance(cls
, *args
, **kwargs
):
1459 Return a random instance of this type of algebra.
1461 This method should be implemented in each subclass.
1463 from sage
.misc
.prandom
import choice
1464 eja_class
= choice(cls
.__subclasses
__())
1466 # These all bubble up to the RationalBasisEJA superclass
1467 # constructor, so any (kw)args valid there are also valid
1469 return eja_class
.random_instance(*args
, **kwargs
)
1474 def dimension_over_reals():
1476 The dimension of this matrix's base ring over the reals.
1478 The reals are dimension one over themselves, obviously; that's
1479 just `\mathbb{R}^{1}`. Likewise, the complex numbers `a + bi`
1480 have dimension two. Finally, the quaternions have dimension
1481 four over the reals.
1483 This is used to determine the size of the matrix returned from
1484 :meth:`real_embed`, among other things.
1486 raise NotImplementedError
1489 def real_embed(cls
,M
):
1491 Embed the matrix ``M`` into a space of real matrices.
1493 The matrix ``M`` can have entries in any field at the moment:
1494 the real numbers, complex numbers, or quaternions. And although
1495 they are not a field, we can probably support octonions at some
1496 point, too. This function returns a real matrix that "acts like"
1497 the original with respect to matrix multiplication; i.e.
1499 real_embed(M*N) = real_embed(M)*real_embed(N)
1502 if M
.ncols() != M
.nrows():
1503 raise ValueError("the matrix 'M' must be square")
1508 def real_unembed(cls
,M
):
1510 The inverse of :meth:`real_embed`.
1512 if M
.ncols() != M
.nrows():
1513 raise ValueError("the matrix 'M' must be square")
1514 if not ZZ(M
.nrows()).mod(cls
.dimension_over_reals()).is_zero():
1515 raise ValueError("the matrix 'M' must be a real embedding")
1519 def jordan_product(X
,Y
):
1520 return (X
*Y
+ Y
*X
)/2
1523 def trace_inner_product(cls
,X
,Y
):
1525 Compute the trace inner-product of two real-embeddings.
1529 sage: from mjo.eja.eja_algebra import (RealSymmetricEJA,
1530 ....: ComplexHermitianEJA,
1531 ....: QuaternionHermitianEJA)
1535 This gives the same answer as it would if we computed the trace
1536 from the unembedded (original) matrices::
1538 sage: set_random_seed()
1539 sage: J = RealSymmetricEJA.random_instance()
1540 sage: x,y = J.random_elements(2)
1541 sage: Xe = x.to_matrix()
1542 sage: Ye = y.to_matrix()
1543 sage: X = J.real_unembed(Xe)
1544 sage: Y = J.real_unembed(Ye)
1545 sage: expected = (X*Y).trace()
1546 sage: actual = J.trace_inner_product(Xe,Ye)
1547 sage: actual == expected
1552 sage: set_random_seed()
1553 sage: J = ComplexHermitianEJA.random_instance()
1554 sage: x,y = J.random_elements(2)
1555 sage: Xe = x.to_matrix()
1556 sage: Ye = y.to_matrix()
1557 sage: X = J.real_unembed(Xe)
1558 sage: Y = J.real_unembed(Ye)
1559 sage: expected = (X*Y).trace().real()
1560 sage: actual = J.trace_inner_product(Xe,Ye)
1561 sage: actual == expected
1566 sage: set_random_seed()
1567 sage: J = QuaternionHermitianEJA.random_instance()
1568 sage: x,y = J.random_elements(2)
1569 sage: Xe = x.to_matrix()
1570 sage: Ye = y.to_matrix()
1571 sage: X = J.real_unembed(Xe)
1572 sage: Y = J.real_unembed(Ye)
1573 sage: expected = (X*Y).trace().coefficient_tuple()[0]
1574 sage: actual = J.trace_inner_product(Xe,Ye)
1575 sage: actual == expected
1579 Xu
= cls
.real_unembed(X
)
1580 Yu
= cls
.real_unembed(Y
)
1581 tr
= (Xu
*Yu
).trace()
1584 # Works in QQ, AA, RDF, et cetera.
1586 except AttributeError:
1587 # A quaternion doesn't have a real() method, but does
1588 # have coefficient_tuple() method that returns the
1589 # coefficients of 1, i, j, and k -- in that order.
1590 return tr
.coefficient_tuple()[0]
1593 class RealMatrixEJA(MatrixEJA
):
1595 def dimension_over_reals():
1599 class RealSymmetricEJA(ConcreteEJA
, RealMatrixEJA
):
1601 The rank-n simple EJA consisting of real symmetric n-by-n
1602 matrices, the usual symmetric Jordan product, and the trace inner
1603 product. It has dimension `(n^2 + n)/2` over the reals.
1607 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1611 sage: J = RealSymmetricEJA(2)
1612 sage: e0, e1, e2 = J.gens()
1620 In theory, our "field" can be any subfield of the reals::
1622 sage: RealSymmetricEJA(2, field=RDF)
1623 Euclidean Jordan algebra of dimension 3 over Real Double Field
1624 sage: RealSymmetricEJA(2, field=RR)
1625 Euclidean Jordan algebra of dimension 3 over Real Field with
1626 53 bits of precision
1630 The dimension of this algebra is `(n^2 + n) / 2`::
1632 sage: set_random_seed()
1633 sage: n_max = RealSymmetricEJA._max_random_instance_size()
1634 sage: n = ZZ.random_element(1, n_max)
1635 sage: J = RealSymmetricEJA(n)
1636 sage: J.dimension() == (n^2 + n)/2
1639 The Jordan multiplication is what we think it is::
1641 sage: set_random_seed()
1642 sage: J = RealSymmetricEJA.random_instance()
1643 sage: x,y = J.random_elements(2)
1644 sage: actual = (x*y).to_matrix()
1645 sage: X = x.to_matrix()
1646 sage: Y = y.to_matrix()
1647 sage: expected = (X*Y + Y*X)/2
1648 sage: actual == expected
1650 sage: J(expected) == x*y
1653 We can change the generator prefix::
1655 sage: RealSymmetricEJA(3, prefix='q').gens()
1656 (q0, q1, q2, q3, q4, q5)
1658 We can construct the (trivial) algebra of rank zero::
1660 sage: RealSymmetricEJA(0)
1661 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1665 def _denormalized_basis(cls
, n
):
1667 Return a basis for the space of real symmetric n-by-n matrices.
1671 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1675 sage: set_random_seed()
1676 sage: n = ZZ.random_element(1,5)
1677 sage: B = RealSymmetricEJA._denormalized_basis(n)
1678 sage: all( M.is_symmetric() for M in B)
1682 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1686 for j
in range(i
+1):
1687 Eij
= matrix(ZZ
, n
, lambda k
,l
: k
==i
and l
==j
)
1691 Sij
= Eij
+ Eij
.transpose()
1697 def _max_random_instance_size():
1698 return 4 # Dimension 10
1701 def random_instance(cls
, **kwargs
):
1703 Return a random instance of this type of algebra.
1705 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1706 return cls(n
, **kwargs
)
1708 def __init__(self
, n
, **kwargs
):
1709 # We know this is a valid EJA, but will double-check
1710 # if the user passes check_axioms=True.
1711 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
1713 super(RealSymmetricEJA
, self
).__init
__(self
._denormalized
_basis
(n
),
1714 self
.jordan_product
,
1715 self
.trace_inner_product
,
1718 # TODO: this could be factored out somehow, but is left here
1719 # because the MatrixEJA is not presently a subclass of the
1720 # FDEJA class that defines rank() and one().
1721 self
.rank
.set_cache(n
)
1722 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
1723 self
.one
.set_cache(self(idV
))
1727 class ComplexMatrixEJA(MatrixEJA
):
1728 # A manual dictionary-cache for the complex_extension() method,
1729 # since apparently @classmethods can't also be @cached_methods.
1730 _complex_extension
= {}
1733 def complex_extension(cls
,field
):
1735 The complex field that we embed/unembed, as an extension
1736 of the given ``field``.
1738 if field
in cls
._complex
_extension
:
1739 return cls
._complex
_extension
[field
]
1741 # Sage doesn't know how to adjoin the complex "i" (the root of
1742 # x^2 + 1) to a field in a general way. Here, we just enumerate
1743 # all of the cases that I have cared to support so far.
1745 # Sage doesn't know how to embed AA into QQbar, i.e. how
1746 # to adjoin sqrt(-1) to AA.
1748 elif not field
.is_exact():
1750 F
= field
.complex_field()
1752 # Works for QQ and... maybe some other fields.
1753 R
= PolynomialRing(field
, 'z')
1755 F
= field
.extension(z
**2 + 1, 'I', embedding
=CLF(-1).sqrt())
1757 cls
._complex
_extension
[field
] = F
1761 def dimension_over_reals():
1765 def real_embed(cls
,M
):
1767 Embed the n-by-n complex matrix ``M`` into the space of real
1768 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1769 bi` to the block matrix ``[[a,b],[-b,a]]``.
1773 sage: from mjo.eja.eja_algebra import ComplexMatrixEJA
1777 sage: F = QuadraticField(-1, 'I')
1778 sage: x1 = F(4 - 2*i)
1779 sage: x2 = F(1 + 2*i)
1782 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1783 sage: ComplexMatrixEJA.real_embed(M)
1792 Embedding is a homomorphism (isomorphism, in fact)::
1794 sage: set_random_seed()
1795 sage: n = ZZ.random_element(3)
1796 sage: F = QuadraticField(-1, 'I')
1797 sage: X = random_matrix(F, n)
1798 sage: Y = random_matrix(F, n)
1799 sage: Xe = ComplexMatrixEJA.real_embed(X)
1800 sage: Ye = ComplexMatrixEJA.real_embed(Y)
1801 sage: XYe = ComplexMatrixEJA.real_embed(X*Y)
1806 super(ComplexMatrixEJA
,cls
).real_embed(M
)
1809 # We don't need any adjoined elements...
1810 field
= M
.base_ring().base_ring()
1816 blocks
.append(matrix(field
, 2, [ [ a
, b
],
1819 return matrix
.block(field
, n
, blocks
)
1823 def real_unembed(cls
,M
):
1825 The inverse of _embed_complex_matrix().
1829 sage: from mjo.eja.eja_algebra import ComplexMatrixEJA
1833 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1834 ....: [-2, 1, -4, 3],
1835 ....: [ 9, 10, 11, 12],
1836 ....: [-10, 9, -12, 11] ])
1837 sage: ComplexMatrixEJA.real_unembed(A)
1839 [ 10*I + 9 12*I + 11]
1843 Unembedding is the inverse of embedding::
1845 sage: set_random_seed()
1846 sage: F = QuadraticField(-1, 'I')
1847 sage: M = random_matrix(F, 3)
1848 sage: Me = ComplexMatrixEJA.real_embed(M)
1849 sage: ComplexMatrixEJA.real_unembed(Me) == M
1853 super(ComplexMatrixEJA
,cls
).real_unembed(M
)
1855 d
= cls
.dimension_over_reals()
1856 F
= cls
.complex_extension(M
.base_ring())
1859 # Go top-left to bottom-right (reading order), converting every
1860 # 2-by-2 block we see to a single complex element.
1862 for k
in range(n
/d
):
1863 for j
in range(n
/d
):
1864 submat
= M
[d
*k
:d
*k
+d
,d
*j
:d
*j
+d
]
1865 if submat
[0,0] != submat
[1,1]:
1866 raise ValueError('bad on-diagonal submatrix')
1867 if submat
[0,1] != -submat
[1,0]:
1868 raise ValueError('bad off-diagonal submatrix')
1869 z
= submat
[0,0] + submat
[0,1]*i
1872 return matrix(F
, n
/d
, elements
)
1875 class ComplexHermitianEJA(ConcreteEJA
, ComplexMatrixEJA
):
1877 The rank-n simple EJA consisting of complex Hermitian n-by-n
1878 matrices over the real numbers, the usual symmetric Jordan product,
1879 and the real-part-of-trace inner product. It has dimension `n^2` over
1884 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1888 In theory, our "field" can be any subfield of the reals::
1890 sage: ComplexHermitianEJA(2, field=RDF)
1891 Euclidean Jordan algebra of dimension 4 over Real Double Field
1892 sage: ComplexHermitianEJA(2, field=RR)
1893 Euclidean Jordan algebra of dimension 4 over Real Field with
1894 53 bits of precision
1898 The dimension of this algebra is `n^2`::
1900 sage: set_random_seed()
1901 sage: n_max = ComplexHermitianEJA._max_random_instance_size()
1902 sage: n = ZZ.random_element(1, n_max)
1903 sage: J = ComplexHermitianEJA(n)
1904 sage: J.dimension() == n^2
1907 The Jordan multiplication is what we think it is::
1909 sage: set_random_seed()
1910 sage: J = ComplexHermitianEJA.random_instance()
1911 sage: x,y = J.random_elements(2)
1912 sage: actual = (x*y).to_matrix()
1913 sage: X = x.to_matrix()
1914 sage: Y = y.to_matrix()
1915 sage: expected = (X*Y + Y*X)/2
1916 sage: actual == expected
1918 sage: J(expected) == x*y
1921 We can change the generator prefix::
1923 sage: ComplexHermitianEJA(2, prefix='z').gens()
1926 We can construct the (trivial) algebra of rank zero::
1928 sage: ComplexHermitianEJA(0)
1929 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1934 def _denormalized_basis(cls
, n
):
1936 Returns a basis for the space of complex Hermitian n-by-n matrices.
1938 Why do we embed these? Basically, because all of numerical linear
1939 algebra assumes that you're working with vectors consisting of `n`
1940 entries from a field and scalars from the same field. There's no way
1941 to tell SageMath that (for example) the vectors contain complex
1942 numbers, while the scalar field is real.
1946 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1950 sage: set_random_seed()
1951 sage: n = ZZ.random_element(1,5)
1952 sage: B = ComplexHermitianEJA._denormalized_basis(n)
1953 sage: all( M.is_symmetric() for M in B)
1958 R
= PolynomialRing(field
, 'z')
1960 F
= field
.extension(z
**2 + 1, 'I')
1963 # This is like the symmetric case, but we need to be careful:
1965 # * We want conjugate-symmetry, not just symmetry.
1966 # * The diagonal will (as a result) be real.
1969 Eij
= matrix
.zero(F
,n
)
1971 for j
in range(i
+1):
1975 Sij
= cls
.real_embed(Eij
)
1978 # The second one has a minus because it's conjugated.
1979 Eij
[j
,i
] = 1 # Eij = Eij + Eij.transpose()
1980 Sij_real
= cls
.real_embed(Eij
)
1982 # Eij = I*Eij - I*Eij.transpose()
1985 Sij_imag
= cls
.real_embed(Eij
)
1991 # Since we embedded these, we can drop back to the "field" that we
1992 # started with instead of the complex extension "F".
1993 return tuple( s
.change_ring(field
) for s
in S
)
1996 def __init__(self
, n
, **kwargs
):
1997 # We know this is a valid EJA, but will double-check
1998 # if the user passes check_axioms=True.
1999 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2001 super(ComplexHermitianEJA
, self
).__init
__(self
._denormalized
_basis
(n
),
2002 self
.jordan_product
,
2003 self
.trace_inner_product
,
2005 # TODO: this could be factored out somehow, but is left here
2006 # because the MatrixEJA is not presently a subclass of the
2007 # FDEJA class that defines rank() and one().
2008 self
.rank
.set_cache(n
)
2009 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
2010 self
.one
.set_cache(self(idV
))
2013 def _max_random_instance_size():
2014 return 3 # Dimension 9
2017 def random_instance(cls
, **kwargs
):
2019 Return a random instance of this type of algebra.
2021 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2022 return cls(n
, **kwargs
)
2024 class QuaternionMatrixEJA(MatrixEJA
):
2026 # A manual dictionary-cache for the quaternion_extension() method,
2027 # since apparently @classmethods can't also be @cached_methods.
2028 _quaternion_extension
= {}
2031 def quaternion_extension(cls
,field
):
2033 The quaternion field that we embed/unembed, as an extension
2034 of the given ``field``.
2036 if field
in cls
._quaternion
_extension
:
2037 return cls
._quaternion
_extension
[field
]
2039 Q
= QuaternionAlgebra(field
,-1,-1)
2041 cls
._quaternion
_extension
[field
] = Q
2045 def dimension_over_reals():
2049 def real_embed(cls
,M
):
2051 Embed the n-by-n quaternion matrix ``M`` into the space of real
2052 matrices of size 4n-by-4n by first sending each quaternion entry `z
2053 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
2054 c+di],[-c + di, a-bi]]`, and then embedding those into a real
2059 sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA
2063 sage: Q = QuaternionAlgebra(QQ,-1,-1)
2064 sage: i,j,k = Q.gens()
2065 sage: x = 1 + 2*i + 3*j + 4*k
2066 sage: M = matrix(Q, 1, [[x]])
2067 sage: QuaternionMatrixEJA.real_embed(M)
2073 Embedding is a homomorphism (isomorphism, in fact)::
2075 sage: set_random_seed()
2076 sage: n = ZZ.random_element(2)
2077 sage: Q = QuaternionAlgebra(QQ,-1,-1)
2078 sage: X = random_matrix(Q, n)
2079 sage: Y = random_matrix(Q, n)
2080 sage: Xe = QuaternionMatrixEJA.real_embed(X)
2081 sage: Ye = QuaternionMatrixEJA.real_embed(Y)
2082 sage: XYe = QuaternionMatrixEJA.real_embed(X*Y)
2087 super(QuaternionMatrixEJA
,cls
).real_embed(M
)
2088 quaternions
= M
.base_ring()
2091 F
= QuadraticField(-1, 'I')
2096 t
= z
.coefficient_tuple()
2101 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
2102 [-c
+ d
*i
, a
- b
*i
]])
2103 realM
= ComplexMatrixEJA
.real_embed(cplxM
)
2104 blocks
.append(realM
)
2106 # We should have real entries by now, so use the realest field
2107 # we've got for the return value.
2108 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
2113 def real_unembed(cls
,M
):
2115 The inverse of _embed_quaternion_matrix().
2119 sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA
2123 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
2124 ....: [-2, 1, -4, 3],
2125 ....: [-3, 4, 1, -2],
2126 ....: [-4, -3, 2, 1]])
2127 sage: QuaternionMatrixEJA.real_unembed(M)
2128 [1 + 2*i + 3*j + 4*k]
2132 Unembedding is the inverse of embedding::
2134 sage: set_random_seed()
2135 sage: Q = QuaternionAlgebra(QQ, -1, -1)
2136 sage: M = random_matrix(Q, 3)
2137 sage: Me = QuaternionMatrixEJA.real_embed(M)
2138 sage: QuaternionMatrixEJA.real_unembed(Me) == M
2142 super(QuaternionMatrixEJA
,cls
).real_unembed(M
)
2144 d
= cls
.dimension_over_reals()
2146 # Use the base ring of the matrix to ensure that its entries can be
2147 # multiplied by elements of the quaternion algebra.
2148 Q
= cls
.quaternion_extension(M
.base_ring())
2151 # Go top-left to bottom-right (reading order), converting every
2152 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
2155 for l
in range(n
/d
):
2156 for m
in range(n
/d
):
2157 submat
= ComplexMatrixEJA
.real_unembed(
2158 M
[d
*l
:d
*l
+d
,d
*m
:d
*m
+d
] )
2159 if submat
[0,0] != submat
[1,1].conjugate():
2160 raise ValueError('bad on-diagonal submatrix')
2161 if submat
[0,1] != -submat
[1,0].conjugate():
2162 raise ValueError('bad off-diagonal submatrix')
2163 z
= submat
[0,0].real()
2164 z
+= submat
[0,0].imag()*i
2165 z
+= submat
[0,1].real()*j
2166 z
+= submat
[0,1].imag()*k
2169 return matrix(Q
, n
/d
, elements
)
2172 class QuaternionHermitianEJA(ConcreteEJA
, QuaternionMatrixEJA
):
2174 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
2175 matrices, the usual symmetric Jordan product, and the
2176 real-part-of-trace inner product. It has dimension `2n^2 - n` over
2181 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2185 In theory, our "field" can be any subfield of the reals::
2187 sage: QuaternionHermitianEJA(2, field=RDF)
2188 Euclidean Jordan algebra of dimension 6 over Real Double Field
2189 sage: QuaternionHermitianEJA(2, field=RR)
2190 Euclidean Jordan algebra of dimension 6 over Real Field with
2191 53 bits of precision
2195 The dimension of this algebra is `2*n^2 - n`::
2197 sage: set_random_seed()
2198 sage: n_max = QuaternionHermitianEJA._max_random_instance_size()
2199 sage: n = ZZ.random_element(1, n_max)
2200 sage: J = QuaternionHermitianEJA(n)
2201 sage: J.dimension() == 2*(n^2) - n
2204 The Jordan multiplication is what we think it is::
2206 sage: set_random_seed()
2207 sage: J = QuaternionHermitianEJA.random_instance()
2208 sage: x,y = J.random_elements(2)
2209 sage: actual = (x*y).to_matrix()
2210 sage: X = x.to_matrix()
2211 sage: Y = y.to_matrix()
2212 sage: expected = (X*Y + Y*X)/2
2213 sage: actual == expected
2215 sage: J(expected) == x*y
2218 We can change the generator prefix::
2220 sage: QuaternionHermitianEJA(2, prefix='a').gens()
2221 (a0, a1, a2, a3, a4, a5)
2223 We can construct the (trivial) algebra of rank zero::
2225 sage: QuaternionHermitianEJA(0)
2226 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2230 def _denormalized_basis(cls
, n
):
2232 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
2234 Why do we embed these? Basically, because all of numerical
2235 linear algebra assumes that you're working with vectors consisting
2236 of `n` entries from a field and scalars from the same field. There's
2237 no way to tell SageMath that (for example) the vectors contain
2238 complex numbers, while the scalar field is real.
2242 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2246 sage: set_random_seed()
2247 sage: n = ZZ.random_element(1,5)
2248 sage: B = QuaternionHermitianEJA._denormalized_basis(n)
2249 sage: all( M.is_symmetric() for M in B )
2254 Q
= QuaternionAlgebra(QQ
,-1,-1)
2257 # This is like the symmetric case, but we need to be careful:
2259 # * We want conjugate-symmetry, not just symmetry.
2260 # * The diagonal will (as a result) be real.
2263 Eij
= matrix
.zero(Q
,n
)
2265 for j
in range(i
+1):
2269 Sij
= cls
.real_embed(Eij
)
2272 # The second, third, and fourth ones have a minus
2273 # because they're conjugated.
2274 # Eij = Eij + Eij.transpose()
2276 Sij_real
= cls
.real_embed(Eij
)
2278 # Eij = I*(Eij - Eij.transpose())
2281 Sij_I
= cls
.real_embed(Eij
)
2283 # Eij = J*(Eij - Eij.transpose())
2286 Sij_J
= cls
.real_embed(Eij
)
2288 # Eij = K*(Eij - Eij.transpose())
2291 Sij_K
= cls
.real_embed(Eij
)
2297 # Since we embedded these, we can drop back to the "field" that we
2298 # started with instead of the quaternion algebra "Q".
2299 return tuple( s
.change_ring(field
) for s
in S
)
2302 def __init__(self
, n
, **kwargs
):
2303 # We know this is a valid EJA, but will double-check
2304 # if the user passes check_axioms=True.
2305 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2307 super(QuaternionHermitianEJA
, self
).__init
__(self
._denormalized
_basis
(n
),
2308 self
.jordan_product
,
2309 self
.trace_inner_product
,
2311 # TODO: this could be factored out somehow, but is left here
2312 # because the MatrixEJA is not presently a subclass of the
2313 # FDEJA class that defines rank() and one().
2314 self
.rank
.set_cache(n
)
2315 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
2316 self
.one
.set_cache(self(idV
))
2320 def _max_random_instance_size():
2322 The maximum rank of a random QuaternionHermitianEJA.
2324 return 2 # Dimension 6
2327 def random_instance(cls
, **kwargs
):
2329 Return a random instance of this type of algebra.
2331 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2332 return cls(n
, **kwargs
)
2335 class HadamardEJA(ConcreteEJA
):
2337 Return the Euclidean Jordan Algebra corresponding to the set
2338 `R^n` under the Hadamard product.
2340 Note: this is nothing more than the Cartesian product of ``n``
2341 copies of the spin algebra. Once Cartesian product algebras
2342 are implemented, this can go.
2346 sage: from mjo.eja.eja_algebra import HadamardEJA
2350 This multiplication table can be verified by hand::
2352 sage: J = HadamardEJA(3)
2353 sage: e0,e1,e2 = J.gens()
2369 We can change the generator prefix::
2371 sage: HadamardEJA(3, prefix='r').gens()
2375 def __init__(self
, n
, **kwargs
):
2377 jordan_product
= lambda x
,y
: x
2378 inner_product
= lambda x
,y
: x
2380 def jordan_product(x
,y
):
2382 return P( xi
*yi
for (xi
,yi
) in zip(x
,y
) )
2384 def inner_product(x
,y
):
2387 # New defaults for keyword arguments. Don't orthonormalize
2388 # because our basis is already orthonormal with respect to our
2389 # inner-product. Don't check the axioms, because we know this
2390 # is a valid EJA... but do double-check if the user passes
2391 # check_axioms=True. Note: we DON'T override the "check_field"
2392 # default here, because the user can pass in a field!
2393 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2394 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2396 column_basis
= tuple( b
.column() for b
in FreeModule(ZZ
, n
).basis() )
2397 super().__init
__(column_basis
,
2402 self
.rank
.set_cache(n
)
2405 self
.one
.set_cache( self
.zero() )
2407 self
.one
.set_cache( sum(self
.gens()) )
2410 def _max_random_instance_size():
2412 The maximum dimension of a random HadamardEJA.
2417 def random_instance(cls
, **kwargs
):
2419 Return a random instance of this type of algebra.
2421 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2422 return cls(n
, **kwargs
)
2425 class BilinearFormEJA(ConcreteEJA
):
2427 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2428 with the half-trace inner product and jordan product ``x*y =
2429 (<Bx,y>,y_bar>, x0*y_bar + y0*x_bar)`` where `B = 1 \times B22` is
2430 a symmetric positive-definite "bilinear form" matrix. Its
2431 dimension is the size of `B`, and it has rank two in dimensions
2432 larger than two. It reduces to the ``JordanSpinEJA`` when `B` is
2433 the identity matrix of order ``n``.
2435 We insist that the one-by-one upper-left identity block of `B` be
2436 passed in as well so that we can be passed a matrix of size zero
2437 to construct a trivial algebra.
2441 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
2442 ....: JordanSpinEJA)
2446 When no bilinear form is specified, the identity matrix is used,
2447 and the resulting algebra is the Jordan spin algebra::
2449 sage: B = matrix.identity(AA,3)
2450 sage: J0 = BilinearFormEJA(B)
2451 sage: J1 = JordanSpinEJA(3)
2452 sage: J0.multiplication_table() == J0.multiplication_table()
2455 An error is raised if the matrix `B` does not correspond to a
2456 positive-definite bilinear form::
2458 sage: B = matrix.random(QQ,2,3)
2459 sage: J = BilinearFormEJA(B)
2460 Traceback (most recent call last):
2462 ValueError: bilinear form is not positive-definite
2463 sage: B = matrix.zero(QQ,3)
2464 sage: J = BilinearFormEJA(B)
2465 Traceback (most recent call last):
2467 ValueError: bilinear form is not positive-definite
2471 We can create a zero-dimensional algebra::
2473 sage: B = matrix.identity(AA,0)
2474 sage: J = BilinearFormEJA(B)
2478 We can check the multiplication condition given in the Jordan, von
2479 Neumann, and Wigner paper (and also discussed on my "On the
2480 symmetry..." paper). Note that this relies heavily on the standard
2481 choice of basis, as does anything utilizing the bilinear form
2482 matrix. We opt not to orthonormalize the basis, because if we
2483 did, we would have to normalize the `s_{i}` in a similar manner::
2485 sage: set_random_seed()
2486 sage: n = ZZ.random_element(5)
2487 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2488 sage: B11 = matrix.identity(QQ,1)
2489 sage: B22 = M.transpose()*M
2490 sage: B = block_matrix(2,2,[ [B11,0 ],
2492 sage: J = BilinearFormEJA(B, orthonormalize=False)
2493 sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
2494 sage: V = J.vector_space()
2495 sage: sis = [ J( V([0] + (M.inverse()*ei).list()).column() )
2496 ....: for ei in eis ]
2497 sage: actual = [ sis[i]*sis[j]
2498 ....: for i in range(n-1)
2499 ....: for j in range(n-1) ]
2500 sage: expected = [ J.one() if i == j else J.zero()
2501 ....: for i in range(n-1)
2502 ....: for j in range(n-1) ]
2503 sage: actual == expected
2507 def __init__(self
, B
, **kwargs
):
2508 # The matrix "B" is supplied by the user in most cases,
2509 # so it makes sense to check whether or not its positive-
2510 # definite unless we are specifically asked not to...
2511 if ("check_axioms" not in kwargs
) or kwargs
["check_axioms"]:
2512 if not B
.is_positive_definite():
2513 raise ValueError("bilinear form is not positive-definite")
2515 # However, all of the other data for this EJA is computed
2516 # by us in manner that guarantees the axioms are
2517 # satisfied. So, again, unless we are specifically asked to
2518 # verify things, we'll skip the rest of the checks.
2519 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2521 def inner_product(x
,y
):
2522 return (y
.T
*B
*x
)[0,0]
2524 def jordan_product(x
,y
):
2530 z0
= inner_product(y
,x
)
2531 zbar
= y0
*xbar
+ x0
*ybar
2532 return P([z0
] + zbar
.list())
2535 column_basis
= tuple( b
.column() for b
in FreeModule(ZZ
, n
).basis() )
2536 super(BilinearFormEJA
, self
).__init
__(column_basis
,
2541 # The rank of this algebra is two, unless we're in a
2542 # one-dimensional ambient space (because the rank is bounded
2543 # by the ambient dimension).
2544 self
.rank
.set_cache(min(n
,2))
2547 self
.one
.set_cache( self
.zero() )
2549 self
.one
.set_cache( self
.monomial(0) )
2552 def _max_random_instance_size():
2554 The maximum dimension of a random BilinearFormEJA.
2559 def random_instance(cls
, **kwargs
):
2561 Return a random instance of this algebra.
2563 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2565 B
= matrix
.identity(ZZ
, n
)
2566 return cls(B
, **kwargs
)
2568 B11
= matrix
.identity(ZZ
, 1)
2569 M
= matrix
.random(ZZ
, n
-1)
2570 I
= matrix
.identity(ZZ
, n
-1)
2572 while alpha
.is_zero():
2573 alpha
= ZZ
.random_element().abs()
2574 B22
= M
.transpose()*M
+ alpha
*I
2576 from sage
.matrix
.special
import block_matrix
2577 B
= block_matrix(2,2, [ [B11
, ZZ(0) ],
2580 return cls(B
, **kwargs
)
2583 class JordanSpinEJA(BilinearFormEJA
):
2585 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2586 with the usual inner product and jordan product ``x*y =
2587 (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2592 sage: from mjo.eja.eja_algebra import JordanSpinEJA
2596 This multiplication table can be verified by hand::
2598 sage: J = JordanSpinEJA(4)
2599 sage: e0,e1,e2,e3 = J.gens()
2615 We can change the generator prefix::
2617 sage: JordanSpinEJA(2, prefix='B').gens()
2622 Ensure that we have the usual inner product on `R^n`::
2624 sage: set_random_seed()
2625 sage: J = JordanSpinEJA.random_instance()
2626 sage: x,y = J.random_elements(2)
2627 sage: actual = x.inner_product(y)
2628 sage: expected = x.to_vector().inner_product(y.to_vector())
2629 sage: actual == expected
2633 def __init__(self
, n
, **kwargs
):
2634 # This is a special case of the BilinearFormEJA with the
2635 # identity matrix as its bilinear form.
2636 B
= matrix
.identity(ZZ
, n
)
2638 # Don't orthonormalize because our basis is already
2639 # orthonormal with respect to our inner-product.
2640 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2642 # But also don't pass check_field=False here, because the user
2643 # can pass in a field!
2644 super(JordanSpinEJA
, self
).__init
__(B
, **kwargs
)
2647 def _max_random_instance_size():
2649 The maximum dimension of a random JordanSpinEJA.
2654 def random_instance(cls
, **kwargs
):
2656 Return a random instance of this type of algebra.
2658 Needed here to override the implementation for ``BilinearFormEJA``.
2660 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2661 return cls(n
, **kwargs
)
2664 class TrivialEJA(ConcreteEJA
):
2666 The trivial Euclidean Jordan algebra consisting of only a zero element.
2670 sage: from mjo.eja.eja_algebra import TrivialEJA
2674 sage: J = TrivialEJA()
2681 sage: 7*J.one()*12*J.one()
2683 sage: J.one().inner_product(J.one())
2685 sage: J.one().norm()
2687 sage: J.one().subalgebra_generated_by()
2688 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2693 def __init__(self
, **kwargs
):
2694 jordan_product
= lambda x
,y
: x
2695 inner_product
= lambda x
,y
: 0
2698 # New defaults for keyword arguments
2699 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2700 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2702 super(TrivialEJA
, self
).__init
__(basis
,
2706 # The rank is zero using my definition, namely the dimension of the
2707 # largest subalgebra generated by any element.
2708 self
.rank
.set_cache(0)
2709 self
.one
.set_cache( self
.zero() )
2712 def random_instance(cls
, **kwargs
):
2713 # We don't take a "size" argument so the superclass method is
2714 # inappropriate for us.
2715 return cls(**kwargs
)
2718 class CartesianProductEJA(CombinatorialFreeModule_CartesianProduct
,
2719 FiniteDimensionalEJA
):
2721 The external (orthogonal) direct sum of two or more Euclidean
2722 Jordan algebras. Every Euclidean Jordan algebra decomposes into an
2723 orthogonal direct sum of simple Euclidean Jordan algebras which is
2724 then isometric to a Cartesian product, so no generality is lost by
2725 providing only this construction.
2729 sage: from mjo.eja.eja_algebra import (random_eja,
2730 ....: CartesianProductEJA,
2732 ....: JordanSpinEJA,
2733 ....: RealSymmetricEJA)
2737 The Jordan product is inherited from our factors and implemented by
2738 our CombinatorialFreeModule Cartesian product superclass::
2740 sage: set_random_seed()
2741 sage: J1 = HadamardEJA(2)
2742 sage: J2 = RealSymmetricEJA(2)
2743 sage: J = cartesian_product([J1,J2])
2744 sage: x,y = J.random_elements(2)
2748 The ability to retrieve the original factors is implemented by our
2749 CombinatorialFreeModule Cartesian product superclass::
2751 sage: J1 = HadamardEJA(2, field=QQ)
2752 sage: J2 = JordanSpinEJA(3, field=QQ)
2753 sage: J = cartesian_product([J1,J2])
2754 sage: J.cartesian_factors()
2755 (Euclidean Jordan algebra of dimension 2 over Rational Field,
2756 Euclidean Jordan algebra of dimension 3 over Rational Field)
2758 You can provide more than two factors::
2760 sage: J1 = HadamardEJA(2)
2761 sage: J2 = JordanSpinEJA(3)
2762 sage: J3 = RealSymmetricEJA(3)
2763 sage: cartesian_product([J1,J2,J3])
2764 Euclidean Jordan algebra of dimension 2 over Algebraic Real
2765 Field (+) Euclidean Jordan algebra of dimension 3 over Algebraic
2766 Real Field (+) Euclidean Jordan algebra of dimension 6 over
2767 Algebraic Real Field
2769 Rank is additive on a Cartesian product::
2771 sage: J1 = HadamardEJA(1)
2772 sage: J2 = RealSymmetricEJA(2)
2773 sage: J = cartesian_product([J1,J2])
2774 sage: J1.rank.clear_cache()
2775 sage: J2.rank.clear_cache()
2776 sage: J.rank.clear_cache()
2779 sage: J.rank() == J1.rank() + J2.rank()
2782 The same rank computation works over the rationals, with whatever
2785 sage: J1 = HadamardEJA(1, field=QQ, orthonormalize=False)
2786 sage: J2 = RealSymmetricEJA(2, field=QQ, orthonormalize=False)
2787 sage: J = cartesian_product([J1,J2])
2788 sage: J1.rank.clear_cache()
2789 sage: J2.rank.clear_cache()
2790 sage: J.rank.clear_cache()
2793 sage: J.rank() == J1.rank() + J2.rank()
2796 The product algebra will be associative if and only if all of its
2797 components are associative::
2799 sage: J1 = HadamardEJA(2)
2800 sage: J1.is_associative()
2802 sage: J2 = HadamardEJA(3)
2803 sage: J2.is_associative()
2805 sage: J3 = RealSymmetricEJA(3)
2806 sage: J3.is_associative()
2808 sage: CP1 = cartesian_product([J1,J2])
2809 sage: CP1.is_associative()
2811 sage: CP2 = cartesian_product([J1,J3])
2812 sage: CP2.is_associative()
2817 All factors must share the same base field::
2819 sage: J1 = HadamardEJA(2, field=QQ)
2820 sage: J2 = RealSymmetricEJA(2)
2821 sage: CartesianProductEJA((J1,J2))
2822 Traceback (most recent call last):
2824 ValueError: all factors must share the same base field
2826 The "cached" Jordan and inner products are the componentwise
2829 sage: set_random_seed()
2830 sage: J1 = random_eja()
2831 sage: J2 = random_eja()
2832 sage: J = cartesian_product([J1,J2])
2833 sage: x,y = J.random_elements(2)
2834 sage: x*y == J.cartesian_jordan_product(x,y)
2836 sage: x.inner_product(y) == J.cartesian_inner_product(x,y)
2839 The cached unit element is the same one that would be computed::
2841 sage: set_random_seed() # long time
2842 sage: J1 = random_eja() # long time
2843 sage: J2 = random_eja() # long time
2844 sage: J = cartesian_product([J1,J2]) # long time
2845 sage: actual = J.one() # long time
2846 sage: J.one.clear_cache() # long time
2847 sage: expected = J.one() # long time
2848 sage: actual == expected # long time
2852 def __init__(self
, algebras
, **kwargs
):
2853 CombinatorialFreeModule_CartesianProduct
.__init
__(self
,
2856 field
= algebras
[0].base_ring()
2857 if not all( J
.base_ring() == field
for J
in algebras
):
2858 raise ValueError("all factors must share the same base field")
2860 associative
= all( m
.is_associative() for m
in algebras
)
2862 # The definition of matrix_space() and self.basis() relies
2863 # only on the stuff in the CFM_CartesianProduct class, which
2864 # we've already initialized.
2865 Js
= self
.cartesian_factors()
2867 MS
= self
.matrix_space()
2869 MS(tuple( self
.cartesian_projection(i
)(b
).to_matrix()
2870 for i
in range(m
) ))
2871 for b
in self
.basis()
2874 # Define jordan/inner products that operate on that matrix_basis.
2875 def jordan_product(x
,y
):
2877 (Js
[i
](x
[i
])*Js
[i
](y
[i
])).to_matrix() for i
in range(m
)
2880 def inner_product(x
, y
):
2882 Js
[i
](x
[i
]).inner_product(Js
[i
](y
[i
])) for i
in range(m
)
2885 # There's no need to check the field since it already came
2886 # from an EJA. Likewise the axioms are guaranteed to be
2887 # satisfied, unless the guy writing this class sucks.
2889 # If you want the basis to be orthonormalized, orthonormalize
2891 FiniteDimensionalEJA
.__init
__(self
,
2896 orthonormalize
=False,
2897 associative
=associative
,
2898 cartesian_product
=True,
2902 ones
= tuple(J
.one() for J
in algebras
)
2903 self
.one
.set_cache(self
._cartesian
_product
_of
_elements
(ones
))
2904 self
.rank
.set_cache(sum(J
.rank() for J
in algebras
))
2906 def matrix_space(self
):
2908 Return the space that our matrix basis lives in as a Cartesian
2913 sage: from mjo.eja.eja_algebra import (HadamardEJA,
2914 ....: RealSymmetricEJA)
2918 sage: J1 = HadamardEJA(1)
2919 sage: J2 = RealSymmetricEJA(2)
2920 sage: J = cartesian_product([J1,J2])
2921 sage: J.matrix_space()
2922 The Cartesian product of (Full MatrixSpace of 1 by 1 dense
2923 matrices over Algebraic Real Field, Full MatrixSpace of 2
2924 by 2 dense matrices over Algebraic Real Field)
2927 from sage
.categories
.cartesian_product
import cartesian_product
2928 return cartesian_product( [J
.matrix_space()
2929 for J
in self
.cartesian_factors()] )
2932 def cartesian_projection(self
, i
):
2936 sage: from mjo.eja.eja_algebra import (random_eja,
2937 ....: JordanSpinEJA,
2939 ....: RealSymmetricEJA,
2940 ....: ComplexHermitianEJA)
2944 The projection morphisms are Euclidean Jordan algebra
2947 sage: J1 = HadamardEJA(2)
2948 sage: J2 = RealSymmetricEJA(2)
2949 sage: J = cartesian_product([J1,J2])
2950 sage: J.cartesian_projection(0)
2951 Linear operator between finite-dimensional Euclidean Jordan
2952 algebras represented by the matrix:
2955 Domain: Euclidean Jordan algebra of dimension 2 over Algebraic
2956 Real Field (+) Euclidean Jordan algebra of dimension 3 over
2957 Algebraic Real Field
2958 Codomain: Euclidean Jordan algebra of dimension 2 over Algebraic
2960 sage: J.cartesian_projection(1)
2961 Linear operator between finite-dimensional Euclidean Jordan
2962 algebras represented by the matrix:
2966 Domain: Euclidean Jordan algebra of dimension 2 over Algebraic
2967 Real Field (+) Euclidean Jordan algebra of dimension 3 over
2968 Algebraic Real Field
2969 Codomain: Euclidean Jordan algebra of dimension 3 over Algebraic
2972 The projections work the way you'd expect on the vector
2973 representation of an element::
2975 sage: J1 = JordanSpinEJA(2)
2976 sage: J2 = ComplexHermitianEJA(2)
2977 sage: J = cartesian_product([J1,J2])
2978 sage: pi_left = J.cartesian_projection(0)
2979 sage: pi_right = J.cartesian_projection(1)
2980 sage: pi_left(J.one()).to_vector()
2982 sage: pi_right(J.one()).to_vector()
2984 sage: J.one().to_vector()
2989 The answer never changes::
2991 sage: set_random_seed()
2992 sage: J1 = random_eja()
2993 sage: J2 = random_eja()
2994 sage: J = cartesian_product([J1,J2])
2995 sage: P0 = J.cartesian_projection(0)
2996 sage: P1 = J.cartesian_projection(0)
3001 Ji
= self
.cartesian_factors()[i
]
3002 # Requires the fix on Trac 31421/31422 to work!
3003 Pi
= super().cartesian_projection(i
)
3004 return FiniteDimensionalEJAOperator(self
,Ji
,Pi
.matrix())
3007 def cartesian_embedding(self
, i
):
3011 sage: from mjo.eja.eja_algebra import (random_eja,
3012 ....: JordanSpinEJA,
3014 ....: RealSymmetricEJA)
3018 The embedding morphisms are Euclidean Jordan algebra
3021 sage: J1 = HadamardEJA(2)
3022 sage: J2 = RealSymmetricEJA(2)
3023 sage: J = cartesian_product([J1,J2])
3024 sage: J.cartesian_embedding(0)
3025 Linear operator between finite-dimensional Euclidean Jordan
3026 algebras represented by the matrix:
3032 Domain: Euclidean Jordan algebra of dimension 2 over
3033 Algebraic Real Field
3034 Codomain: Euclidean Jordan algebra of dimension 2 over
3035 Algebraic Real Field (+) Euclidean Jordan algebra of
3036 dimension 3 over Algebraic Real Field
3037 sage: J.cartesian_embedding(1)
3038 Linear operator between finite-dimensional Euclidean Jordan
3039 algebras represented by the matrix:
3045 Domain: Euclidean Jordan algebra of dimension 3 over
3046 Algebraic Real Field
3047 Codomain: Euclidean Jordan algebra of dimension 2 over
3048 Algebraic Real Field (+) Euclidean Jordan algebra of
3049 dimension 3 over Algebraic Real Field
3051 The embeddings work the way you'd expect on the vector
3052 representation of an element::
3054 sage: J1 = JordanSpinEJA(3)
3055 sage: J2 = RealSymmetricEJA(2)
3056 sage: J = cartesian_product([J1,J2])
3057 sage: iota_left = J.cartesian_embedding(0)
3058 sage: iota_right = J.cartesian_embedding(1)
3059 sage: iota_left(J1.zero()) == J.zero()
3061 sage: iota_right(J2.zero()) == J.zero()
3063 sage: J1.one().to_vector()
3065 sage: iota_left(J1.one()).to_vector()
3067 sage: J2.one().to_vector()
3069 sage: iota_right(J2.one()).to_vector()
3071 sage: J.one().to_vector()
3076 The answer never changes::
3078 sage: set_random_seed()
3079 sage: J1 = random_eja()
3080 sage: J2 = random_eja()
3081 sage: J = cartesian_product([J1,J2])
3082 sage: E0 = J.cartesian_embedding(0)
3083 sage: E1 = J.cartesian_embedding(0)
3087 Composing a projection with the corresponding inclusion should
3088 produce the identity map, and mismatching them should produce
3091 sage: set_random_seed()
3092 sage: J1 = random_eja()
3093 sage: J2 = random_eja()
3094 sage: J = cartesian_product([J1,J2])
3095 sage: iota_left = J.cartesian_embedding(0)
3096 sage: iota_right = J.cartesian_embedding(1)
3097 sage: pi_left = J.cartesian_projection(0)
3098 sage: pi_right = J.cartesian_projection(1)
3099 sage: pi_left*iota_left == J1.one().operator()
3101 sage: pi_right*iota_right == J2.one().operator()
3103 sage: (pi_left*iota_right).is_zero()
3105 sage: (pi_right*iota_left).is_zero()
3109 Ji
= self
.cartesian_factors()[i
]
3110 # Requires the fix on Trac 31421/31422 to work!
3111 Ei
= super().cartesian_embedding(i
)
3112 return FiniteDimensionalEJAOperator(Ji
,self
,Ei
.matrix())
3115 def cartesian_jordan_product(self
, x
, y
):
3117 The componentwise Jordan product.
3119 We project ``x`` and ``y`` onto our factors, and add up the
3120 Jordan products from the subalgebras. This may still be useful
3121 after (if) the default Jordan product in the Cartesian product
3122 algebra is overridden.
3126 sage: from mjo.eja.eja_algebra import (HadamardEJA,
3127 ....: JordanSpinEJA)
3131 sage: J1 = HadamardEJA(3)
3132 sage: J2 = JordanSpinEJA(3)
3133 sage: J = cartesian_product([J1,J2])
3134 sage: x1 = J1.from_vector(vector(QQ,(1,2,1)))
3135 sage: y1 = J1.from_vector(vector(QQ,(1,0,2)))
3136 sage: x2 = J2.from_vector(vector(QQ,(1,2,3)))
3137 sage: y2 = J2.from_vector(vector(QQ,(1,1,1)))
3138 sage: z1 = J.from_vector(vector(QQ,(1,2,1,1,2,3)))
3139 sage: z2 = J.from_vector(vector(QQ,(1,0,2,1,1,1)))
3140 sage: (x1*y1).to_vector()
3142 sage: (x2*y2).to_vector()
3144 sage: J.cartesian_jordan_product(z1,z2).to_vector()
3148 m
= len(self
.cartesian_factors())
3149 projections
= ( self
.cartesian_projection(i
) for i
in range(m
) )
3150 products
= ( P(x
)*P(y
) for P
in projections
)
3151 return self
._cartesian
_product
_of
_elements
(tuple(products
))
3153 def cartesian_inner_product(self
, x
, y
):
3155 The standard componentwise Cartesian inner-product.
3157 We project ``x`` and ``y`` onto our factors, and add up the
3158 inner-products from the subalgebras. This may still be useful
3159 after (if) the default inner product in the Cartesian product
3160 algebra is overridden.
3164 sage: from mjo.eja.eja_algebra import (HadamardEJA,
3165 ....: QuaternionHermitianEJA)
3169 sage: J1 = HadamardEJA(3,field=QQ)
3170 sage: J2 = QuaternionHermitianEJA(2,field=QQ,orthonormalize=False)
3171 sage: J = cartesian_product([J1,J2])
3176 sage: x1.inner_product(x2)
3178 sage: y1.inner_product(y2)
3180 sage: z1 = J._cartesian_product_of_elements((x1,y1))
3181 sage: z2 = J._cartesian_product_of_elements((x2,y2))
3182 sage: J.cartesian_inner_product(z1,z2)
3186 m
= len(self
.cartesian_factors())
3187 projections
= ( self
.cartesian_projection(i
) for i
in range(m
) )
3188 return sum( P(x
).inner_product(P(y
)) for P
in projections
)
3191 def _element_constructor_(self
, elt
):
3193 Construct an element of this algebra from an ordered tuple.
3195 We just apply the element constructor from each of our factors
3196 to the corresponding component of the tuple, and package up
3201 sage: from mjo.eja.eja_algebra import (HadamardEJA,
3202 ....: RealSymmetricEJA)
3206 sage: J1 = HadamardEJA(3)
3207 sage: J2 = RealSymmetricEJA(2)
3208 sage: J = cartesian_product([J1,J2])
3209 sage: J( (J1.matrix_basis()[1], J2.matrix_basis()[2]) )
3212 m
= len(self
.cartesian_factors())
3214 z
= tuple( self
.cartesian_factors()[i
](elt
[i
]) for i
in range(m
) )
3215 return self
._cartesian
_product
_of
_elements
(z
)
3217 raise ValueError("not an element of this algebra")
3219 Element
= CartesianProductEJAElement
3222 FiniteDimensionalEJA
.CartesianProduct
= CartesianProductEJA
3224 random_eja
= ConcreteEJA
.random_instance
3225 #def random_eja(*args, **kwargs):
3226 # from sage.categories.cartesian_product import cartesian_product
3227 # J1 = HadamardEJA(1, **kwargs)
3228 # J2 = RealSymmetricEJA(2, **kwargs)
3229 # J = cartesian_product([J1,J2])