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 # Default these to what they should be if they turn out to be
1033 # trivial, because eigenspaces_left() won't return eigenvalues
1034 # corresponding to trivial spaces (e.g. it returns only the
1035 # eigenspace corresponding to lambda=1 if you take the
1036 # decomposition relative to the identity element).
1037 trivial
= self
.subalgebra(())
1038 J0
= trivial
# eigenvalue zero
1039 J5
= VectorSpace(self
.base_ring(), 0) # eigenvalue one-half
1040 J1
= trivial
# eigenvalue one
1042 for (eigval
, eigspace
) in c
.operator().matrix().right_eigenspaces():
1043 if eigval
== ~
(self
.base_ring()(2)):
1046 gens
= tuple( self
.from_vector(b
) for b
in eigspace
.basis() )
1047 subalg
= self
.subalgebra(gens
, check_axioms
=False)
1053 raise ValueError("unexpected eigenvalue: %s" % eigval
)
1058 def random_element(self
, thorough
=False):
1060 Return a random element of this algebra.
1062 Our algebra superclass method only returns a linear
1063 combination of at most two basis elements. We instead
1064 want the vector space "random element" method that
1065 returns a more diverse selection.
1069 - ``thorough`` -- (boolean; default False) whether or not we
1070 should generate irrational coefficients for the random
1071 element when our base ring is irrational; this slows the
1072 algebra operations to a crawl, but any truly random method
1076 # For a general base ring... maybe we can trust this to do the
1077 # right thing? Unlikely, but.
1078 V
= self
.vector_space()
1079 v
= V
.random_element()
1081 if self
.base_ring() is AA
:
1082 # The "random element" method of the algebraic reals is
1083 # stupid at the moment, and only returns integers between
1084 # -2 and 2, inclusive:
1086 # https://trac.sagemath.org/ticket/30875
1088 # Instead, we implement our own "random vector" method,
1089 # and then coerce that into the algebra. We use the vector
1090 # space degree here instead of the dimension because a
1091 # subalgebra could (for example) be spanned by only two
1092 # vectors, each with five coordinates. We need to
1093 # generate all five coordinates.
1095 v
*= QQbar
.random_element().real()
1097 v
*= QQ
.random_element()
1099 return self
.from_vector(V
.coordinate_vector(v
))
1101 def random_elements(self
, count
, thorough
=False):
1103 Return ``count`` random elements as a tuple.
1107 - ``thorough`` -- (boolean; default False) whether or not we
1108 should generate irrational coefficients for the random
1109 elements when our base ring is irrational; this slows the
1110 algebra operations to a crawl, but any truly random method
1115 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1119 sage: J = JordanSpinEJA(3)
1120 sage: x,y,z = J.random_elements(3)
1121 sage: all( [ x in J, y in J, z in J ])
1123 sage: len( J.random_elements(10) ) == 10
1127 return tuple( self
.random_element(thorough
)
1128 for idx
in range(count
) )
1132 def _charpoly_coefficients(self
):
1134 The `r` polynomial coefficients of the "characteristic polynomial
1139 sage: from mjo.eja.eja_algebra import random_eja
1143 The theory shows that these are all homogeneous polynomials of
1146 sage: set_random_seed()
1147 sage: J = random_eja()
1148 sage: all(p.is_homogeneous() for p in J._charpoly_coefficients())
1152 n
= self
.dimension()
1153 R
= self
.coordinate_polynomial_ring()
1155 F
= R
.fraction_field()
1158 # From a result in my book, these are the entries of the
1159 # basis representation of L_x.
1160 return sum( vars[k
]*self
.gens()[k
].operator().matrix()[i
,j
]
1163 L_x
= matrix(F
, n
, n
, L_x_i_j
)
1166 if self
.rank
.is_in_cache():
1168 # There's no need to pad the system with redundant
1169 # columns if we *know* they'll be redundant.
1172 # Compute an extra power in case the rank is equal to
1173 # the dimension (otherwise, we would stop at x^(r-1)).
1174 x_powers
= [ (L_x
**k
)*self
.one().to_vector()
1175 for k
in range(n
+1) ]
1176 A
= matrix
.column(F
, x_powers
[:n
])
1177 AE
= A
.extended_echelon_form()
1184 # The theory says that only the first "r" coefficients are
1185 # nonzero, and they actually live in the original polynomial
1186 # ring and not the fraction field. We negate them because in
1187 # the actual characteristic polynomial, they get moved to the
1188 # other side where x^r lives. We don't bother to trim A_rref
1189 # down to a square matrix and solve the resulting system,
1190 # because the upper-left r-by-r portion of A_rref is
1191 # guaranteed to be the identity matrix, so e.g.
1193 # A_rref.solve_right(Y)
1195 # would just be returning Y.
1196 return (-E
*b
)[:r
].change_ring(R
)
1201 Return the rank of this EJA.
1203 This is a cached method because we know the rank a priori for
1204 all of the algebras we can construct. Thus we can avoid the
1205 expensive ``_charpoly_coefficients()`` call unless we truly
1206 need to compute the whole characteristic polynomial.
1210 sage: from mjo.eja.eja_algebra import (HadamardEJA,
1211 ....: JordanSpinEJA,
1212 ....: RealSymmetricEJA,
1213 ....: ComplexHermitianEJA,
1214 ....: QuaternionHermitianEJA,
1219 The rank of the Jordan spin algebra is always two::
1221 sage: JordanSpinEJA(2).rank()
1223 sage: JordanSpinEJA(3).rank()
1225 sage: JordanSpinEJA(4).rank()
1228 The rank of the `n`-by-`n` Hermitian real, complex, or
1229 quaternion matrices is `n`::
1231 sage: RealSymmetricEJA(4).rank()
1233 sage: ComplexHermitianEJA(3).rank()
1235 sage: QuaternionHermitianEJA(2).rank()
1240 Ensure that every EJA that we know how to construct has a
1241 positive integer rank, unless the algebra is trivial in
1242 which case its rank will be zero::
1244 sage: set_random_seed()
1245 sage: J = random_eja()
1249 sage: r > 0 or (r == 0 and J.is_trivial())
1252 Ensure that computing the rank actually works, since the ranks
1253 of all simple algebras are known and will be cached by default::
1255 sage: set_random_seed() # long time
1256 sage: J = random_eja() # long time
1257 sage: cached = J.rank() # long time
1258 sage: J.rank.clear_cache() # long time
1259 sage: J.rank() == cached # long time
1263 return len(self
._charpoly
_coefficients
())
1266 def subalgebra(self
, basis
, **kwargs
):
1268 Create a subalgebra of this algebra from the given basis.
1270 This is a simple wrapper around a subalgebra class constructor
1271 that can be overridden in subclasses.
1273 from mjo
.eja
.eja_subalgebra
import FiniteDimensionalEJASubalgebra
1274 return FiniteDimensionalEJASubalgebra(self
, basis
, **kwargs
)
1277 def vector_space(self
):
1279 Return the vector space that underlies this algebra.
1283 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1287 sage: J = RealSymmetricEJA(2)
1288 sage: J.vector_space()
1289 Vector space of dimension 3 over...
1292 return self
.zero().to_vector().parent().ambient_vector_space()
1295 Element
= FiniteDimensionalEJAElement
1297 class RationalBasisEJA(FiniteDimensionalEJA
):
1299 New class for algebras whose supplied basis elements have all rational entries.
1303 sage: from mjo.eja.eja_algebra import BilinearFormEJA
1307 The supplied basis is orthonormalized by default::
1309 sage: B = matrix(QQ, [[1, 0, 0], [0, 25, -32], [0, -32, 41]])
1310 sage: J = BilinearFormEJA(B)
1311 sage: J.matrix_basis()
1328 # Abuse the check_field parameter to check that the entries of
1329 # out basis (in ambient coordinates) are in the field QQ.
1330 if not all( all(b_i
in QQ
for b_i
in b
.list()) for b
in basis
):
1331 raise TypeError("basis not rational")
1333 self
._rational
_algebra
= None
1335 # There's no point in constructing the extra algebra if this
1336 # one is already rational.
1338 # Note: the same Jordan and inner-products work here,
1339 # because they are necessarily defined with respect to
1340 # ambient coordinates and not any particular basis.
1341 self
._rational
_algebra
= FiniteDimensionalEJA(
1346 orthonormalize
=False,
1350 super().__init
__(basis
,
1354 check_field
=check_field
,
1358 def _charpoly_coefficients(self
):
1362 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
1363 ....: JordanSpinEJA)
1367 The base ring of the resulting polynomial coefficients is what
1368 it should be, and not the rationals (unless the algebra was
1369 already over the rationals)::
1371 sage: J = JordanSpinEJA(3)
1372 sage: J._charpoly_coefficients()
1373 (X1^2 - X2^2 - X3^2, -2*X1)
1374 sage: a0 = J._charpoly_coefficients()[0]
1376 Algebraic Real Field
1377 sage: a0.base_ring()
1378 Algebraic Real Field
1381 if self
._rational
_algebra
is None:
1382 # There's no need to construct *another* algebra over the
1383 # rationals if this one is already over the
1384 # rationals. Likewise, if we never orthonormalized our
1385 # basis, we might as well just use the given one.
1386 return super()._charpoly
_coefficients
()
1388 # Do the computation over the rationals. The answer will be
1389 # the same, because all we've done is a change of basis.
1390 # Then, change back from QQ to our real base ring
1391 a
= ( a_i
.change_ring(self
.base_ring())
1392 for a_i
in self
._rational
_algebra
._charpoly
_coefficients
() )
1394 if self
._deortho
_matrix
is None:
1395 # This can happen if our base ring was, say, AA and we
1396 # chose not to (or didn't need to) orthonormalize. It's
1397 # still faster to do the computations over QQ even if
1398 # the numbers in the boxes stay the same.
1401 # Otherwise, convert the coordinate variables back to the
1402 # deorthonormalized ones.
1403 R
= self
.coordinate_polynomial_ring()
1404 from sage
.modules
.free_module_element
import vector
1405 X
= vector(R
, R
.gens())
1406 BX
= self
._deortho
_matrix
*X
1408 subs_dict
= { X[i]: BX[i] for i in range(len(X)) }
1409 return tuple( a_i
.subs(subs_dict
) for a_i
in a
)
1411 class ConcreteEJA(RationalBasisEJA
):
1413 A class for the Euclidean Jordan algebras that we know by name.
1415 These are the Jordan algebras whose basis, multiplication table,
1416 rank, and so on are known a priori. More to the point, they are
1417 the Euclidean Jordan algebras for which we are able to conjure up
1418 a "random instance."
1422 sage: from mjo.eja.eja_algebra import ConcreteEJA
1426 Our basis is normalized with respect to the algebra's inner
1427 product, unless we specify otherwise::
1429 sage: set_random_seed()
1430 sage: J = ConcreteEJA.random_instance()
1431 sage: all( b.norm() == 1 for b in J.gens() )
1434 Since our basis is orthonormal with respect to the algebra's inner
1435 product, and since we know that this algebra is an EJA, any
1436 left-multiplication operator's matrix will be symmetric because
1437 natural->EJA basis representation is an isometry and within the
1438 EJA the operator is self-adjoint by the Jordan axiom::
1440 sage: set_random_seed()
1441 sage: J = ConcreteEJA.random_instance()
1442 sage: x = J.random_element()
1443 sage: x.operator().is_self_adjoint()
1448 def _max_random_instance_size():
1450 Return an integer "size" that is an upper bound on the size of
1451 this algebra when it is used in a random test
1452 case. Unfortunately, the term "size" is ambiguous -- when
1453 dealing with `R^n` under either the Hadamard or Jordan spin
1454 product, the "size" refers to the dimension `n`. When dealing
1455 with a matrix algebra (real symmetric or complex/quaternion
1456 Hermitian), it refers to the size of the matrix, which is far
1457 less than the dimension of the underlying vector space.
1459 This method must be implemented in each subclass.
1461 raise NotImplementedError
1464 def random_instance(cls
, *args
, **kwargs
):
1466 Return a random instance of this type of algebra.
1468 This method should be implemented in each subclass.
1470 from sage
.misc
.prandom
import choice
1471 eja_class
= choice(cls
.__subclasses
__())
1473 # These all bubble up to the RationalBasisEJA superclass
1474 # constructor, so any (kw)args valid there are also valid
1476 return eja_class
.random_instance(*args
, **kwargs
)
1481 def dimension_over_reals():
1483 The dimension of this matrix's base ring over the reals.
1485 The reals are dimension one over themselves, obviously; that's
1486 just `\mathbb{R}^{1}`. Likewise, the complex numbers `a + bi`
1487 have dimension two. Finally, the quaternions have dimension
1488 four over the reals.
1490 This is used to determine the size of the matrix returned from
1491 :meth:`real_embed`, among other things.
1493 raise NotImplementedError
1496 def real_embed(cls
,M
):
1498 Embed the matrix ``M`` into a space of real matrices.
1500 The matrix ``M`` can have entries in any field at the moment:
1501 the real numbers, complex numbers, or quaternions. And although
1502 they are not a field, we can probably support octonions at some
1503 point, too. This function returns a real matrix that "acts like"
1504 the original with respect to matrix multiplication; i.e.
1506 real_embed(M*N) = real_embed(M)*real_embed(N)
1509 if M
.ncols() != M
.nrows():
1510 raise ValueError("the matrix 'M' must be square")
1515 def real_unembed(cls
,M
):
1517 The inverse of :meth:`real_embed`.
1519 if M
.ncols() != M
.nrows():
1520 raise ValueError("the matrix 'M' must be square")
1521 if not ZZ(M
.nrows()).mod(cls
.dimension_over_reals()).is_zero():
1522 raise ValueError("the matrix 'M' must be a real embedding")
1526 def jordan_product(X
,Y
):
1527 return (X
*Y
+ Y
*X
)/2
1530 def trace_inner_product(cls
,X
,Y
):
1532 Compute the trace inner-product of two real-embeddings.
1536 sage: from mjo.eja.eja_algebra import (RealSymmetricEJA,
1537 ....: ComplexHermitianEJA,
1538 ....: QuaternionHermitianEJA)
1542 This gives the same answer as it would if we computed the trace
1543 from the unembedded (original) matrices::
1545 sage: set_random_seed()
1546 sage: J = RealSymmetricEJA.random_instance()
1547 sage: x,y = J.random_elements(2)
1548 sage: Xe = x.to_matrix()
1549 sage: Ye = y.to_matrix()
1550 sage: X = J.real_unembed(Xe)
1551 sage: Y = J.real_unembed(Ye)
1552 sage: expected = (X*Y).trace()
1553 sage: actual = J.trace_inner_product(Xe,Ye)
1554 sage: actual == expected
1559 sage: set_random_seed()
1560 sage: J = ComplexHermitianEJA.random_instance()
1561 sage: x,y = J.random_elements(2)
1562 sage: Xe = x.to_matrix()
1563 sage: Ye = y.to_matrix()
1564 sage: X = J.real_unembed(Xe)
1565 sage: Y = J.real_unembed(Ye)
1566 sage: expected = (X*Y).trace().real()
1567 sage: actual = J.trace_inner_product(Xe,Ye)
1568 sage: actual == expected
1573 sage: set_random_seed()
1574 sage: J = QuaternionHermitianEJA.random_instance()
1575 sage: x,y = J.random_elements(2)
1576 sage: Xe = x.to_matrix()
1577 sage: Ye = y.to_matrix()
1578 sage: X = J.real_unembed(Xe)
1579 sage: Y = J.real_unembed(Ye)
1580 sage: expected = (X*Y).trace().coefficient_tuple()[0]
1581 sage: actual = J.trace_inner_product(Xe,Ye)
1582 sage: actual == expected
1586 Xu
= cls
.real_unembed(X
)
1587 Yu
= cls
.real_unembed(Y
)
1588 tr
= (Xu
*Yu
).trace()
1591 # Works in QQ, AA, RDF, et cetera.
1593 except AttributeError:
1594 # A quaternion doesn't have a real() method, but does
1595 # have coefficient_tuple() method that returns the
1596 # coefficients of 1, i, j, and k -- in that order.
1597 return tr
.coefficient_tuple()[0]
1600 class RealMatrixEJA(MatrixEJA
):
1602 def dimension_over_reals():
1606 class RealSymmetricEJA(ConcreteEJA
, RealMatrixEJA
):
1608 The rank-n simple EJA consisting of real symmetric n-by-n
1609 matrices, the usual symmetric Jordan product, and the trace inner
1610 product. It has dimension `(n^2 + n)/2` over the reals.
1614 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1618 sage: J = RealSymmetricEJA(2)
1619 sage: e0, e1, e2 = J.gens()
1627 In theory, our "field" can be any subfield of the reals::
1629 sage: RealSymmetricEJA(2, field=RDF)
1630 Euclidean Jordan algebra of dimension 3 over Real Double Field
1631 sage: RealSymmetricEJA(2, field=RR)
1632 Euclidean Jordan algebra of dimension 3 over Real Field with
1633 53 bits of precision
1637 The dimension of this algebra is `(n^2 + n) / 2`::
1639 sage: set_random_seed()
1640 sage: n_max = RealSymmetricEJA._max_random_instance_size()
1641 sage: n = ZZ.random_element(1, n_max)
1642 sage: J = RealSymmetricEJA(n)
1643 sage: J.dimension() == (n^2 + n)/2
1646 The Jordan multiplication is what we think it is::
1648 sage: set_random_seed()
1649 sage: J = RealSymmetricEJA.random_instance()
1650 sage: x,y = J.random_elements(2)
1651 sage: actual = (x*y).to_matrix()
1652 sage: X = x.to_matrix()
1653 sage: Y = y.to_matrix()
1654 sage: expected = (X*Y + Y*X)/2
1655 sage: actual == expected
1657 sage: J(expected) == x*y
1660 We can change the generator prefix::
1662 sage: RealSymmetricEJA(3, prefix='q').gens()
1663 (q0, q1, q2, q3, q4, q5)
1665 We can construct the (trivial) algebra of rank zero::
1667 sage: RealSymmetricEJA(0)
1668 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1672 def _denormalized_basis(cls
, n
):
1674 Return a basis for the space of real symmetric n-by-n matrices.
1678 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1682 sage: set_random_seed()
1683 sage: n = ZZ.random_element(1,5)
1684 sage: B = RealSymmetricEJA._denormalized_basis(n)
1685 sage: all( M.is_symmetric() for M in B)
1689 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1693 for j
in range(i
+1):
1694 Eij
= matrix(ZZ
, n
, lambda k
,l
: k
==i
and l
==j
)
1698 Sij
= Eij
+ Eij
.transpose()
1704 def _max_random_instance_size():
1705 return 4 # Dimension 10
1708 def random_instance(cls
, **kwargs
):
1710 Return a random instance of this type of algebra.
1712 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1713 return cls(n
, **kwargs
)
1715 def __init__(self
, n
, **kwargs
):
1716 # We know this is a valid EJA, but will double-check
1717 # if the user passes check_axioms=True.
1718 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
1720 super(RealSymmetricEJA
, self
).__init
__(self
._denormalized
_basis
(n
),
1721 self
.jordan_product
,
1722 self
.trace_inner_product
,
1725 # TODO: this could be factored out somehow, but is left here
1726 # because the MatrixEJA is not presently a subclass of the
1727 # FDEJA class that defines rank() and one().
1728 self
.rank
.set_cache(n
)
1729 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
1730 self
.one
.set_cache(self(idV
))
1734 class ComplexMatrixEJA(MatrixEJA
):
1735 # A manual dictionary-cache for the complex_extension() method,
1736 # since apparently @classmethods can't also be @cached_methods.
1737 _complex_extension
= {}
1740 def complex_extension(cls
,field
):
1742 The complex field that we embed/unembed, as an extension
1743 of the given ``field``.
1745 if field
in cls
._complex
_extension
:
1746 return cls
._complex
_extension
[field
]
1748 # Sage doesn't know how to adjoin the complex "i" (the root of
1749 # x^2 + 1) to a field in a general way. Here, we just enumerate
1750 # all of the cases that I have cared to support so far.
1752 # Sage doesn't know how to embed AA into QQbar, i.e. how
1753 # to adjoin sqrt(-1) to AA.
1755 elif not field
.is_exact():
1757 F
= field
.complex_field()
1759 # Works for QQ and... maybe some other fields.
1760 R
= PolynomialRing(field
, 'z')
1762 F
= field
.extension(z
**2 + 1, 'I', embedding
=CLF(-1).sqrt())
1764 cls
._complex
_extension
[field
] = F
1768 def dimension_over_reals():
1772 def real_embed(cls
,M
):
1774 Embed the n-by-n complex matrix ``M`` into the space of real
1775 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1776 bi` to the block matrix ``[[a,b],[-b,a]]``.
1780 sage: from mjo.eja.eja_algebra import ComplexMatrixEJA
1784 sage: F = QuadraticField(-1, 'I')
1785 sage: x1 = F(4 - 2*i)
1786 sage: x2 = F(1 + 2*i)
1789 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1790 sage: ComplexMatrixEJA.real_embed(M)
1799 Embedding is a homomorphism (isomorphism, in fact)::
1801 sage: set_random_seed()
1802 sage: n = ZZ.random_element(3)
1803 sage: F = QuadraticField(-1, 'I')
1804 sage: X = random_matrix(F, n)
1805 sage: Y = random_matrix(F, n)
1806 sage: Xe = ComplexMatrixEJA.real_embed(X)
1807 sage: Ye = ComplexMatrixEJA.real_embed(Y)
1808 sage: XYe = ComplexMatrixEJA.real_embed(X*Y)
1813 super(ComplexMatrixEJA
,cls
).real_embed(M
)
1816 # We don't need any adjoined elements...
1817 field
= M
.base_ring().base_ring()
1823 blocks
.append(matrix(field
, 2, [ [ a
, b
],
1826 return matrix
.block(field
, n
, blocks
)
1830 def real_unembed(cls
,M
):
1832 The inverse of _embed_complex_matrix().
1836 sage: from mjo.eja.eja_algebra import ComplexMatrixEJA
1840 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1841 ....: [-2, 1, -4, 3],
1842 ....: [ 9, 10, 11, 12],
1843 ....: [-10, 9, -12, 11] ])
1844 sage: ComplexMatrixEJA.real_unembed(A)
1846 [ 10*I + 9 12*I + 11]
1850 Unembedding is the inverse of embedding::
1852 sage: set_random_seed()
1853 sage: F = QuadraticField(-1, 'I')
1854 sage: M = random_matrix(F, 3)
1855 sage: Me = ComplexMatrixEJA.real_embed(M)
1856 sage: ComplexMatrixEJA.real_unembed(Me) == M
1860 super(ComplexMatrixEJA
,cls
).real_unembed(M
)
1862 d
= cls
.dimension_over_reals()
1863 F
= cls
.complex_extension(M
.base_ring())
1866 # Go top-left to bottom-right (reading order), converting every
1867 # 2-by-2 block we see to a single complex element.
1869 for k
in range(n
/d
):
1870 for j
in range(n
/d
):
1871 submat
= M
[d
*k
:d
*k
+d
,d
*j
:d
*j
+d
]
1872 if submat
[0,0] != submat
[1,1]:
1873 raise ValueError('bad on-diagonal submatrix')
1874 if submat
[0,1] != -submat
[1,0]:
1875 raise ValueError('bad off-diagonal submatrix')
1876 z
= submat
[0,0] + submat
[0,1]*i
1879 return matrix(F
, n
/d
, elements
)
1882 class ComplexHermitianEJA(ConcreteEJA
, ComplexMatrixEJA
):
1884 The rank-n simple EJA consisting of complex Hermitian n-by-n
1885 matrices over the real numbers, the usual symmetric Jordan product,
1886 and the real-part-of-trace inner product. It has dimension `n^2` over
1891 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1895 In theory, our "field" can be any subfield of the reals::
1897 sage: ComplexHermitianEJA(2, field=RDF)
1898 Euclidean Jordan algebra of dimension 4 over Real Double Field
1899 sage: ComplexHermitianEJA(2, field=RR)
1900 Euclidean Jordan algebra of dimension 4 over Real Field with
1901 53 bits of precision
1905 The dimension of this algebra is `n^2`::
1907 sage: set_random_seed()
1908 sage: n_max = ComplexHermitianEJA._max_random_instance_size()
1909 sage: n = ZZ.random_element(1, n_max)
1910 sage: J = ComplexHermitianEJA(n)
1911 sage: J.dimension() == n^2
1914 The Jordan multiplication is what we think it is::
1916 sage: set_random_seed()
1917 sage: J = ComplexHermitianEJA.random_instance()
1918 sage: x,y = J.random_elements(2)
1919 sage: actual = (x*y).to_matrix()
1920 sage: X = x.to_matrix()
1921 sage: Y = y.to_matrix()
1922 sage: expected = (X*Y + Y*X)/2
1923 sage: actual == expected
1925 sage: J(expected) == x*y
1928 We can change the generator prefix::
1930 sage: ComplexHermitianEJA(2, prefix='z').gens()
1933 We can construct the (trivial) algebra of rank zero::
1935 sage: ComplexHermitianEJA(0)
1936 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1941 def _denormalized_basis(cls
, n
):
1943 Returns a basis for the space of complex Hermitian n-by-n matrices.
1945 Why do we embed these? Basically, because all of numerical linear
1946 algebra assumes that you're working with vectors consisting of `n`
1947 entries from a field and scalars from the same field. There's no way
1948 to tell SageMath that (for example) the vectors contain complex
1949 numbers, while the scalar field is real.
1953 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1957 sage: set_random_seed()
1958 sage: n = ZZ.random_element(1,5)
1959 sage: B = ComplexHermitianEJA._denormalized_basis(n)
1960 sage: all( M.is_symmetric() for M in B)
1965 R
= PolynomialRing(field
, 'z')
1967 F
= field
.extension(z
**2 + 1, 'I')
1970 # This is like the symmetric case, but we need to be careful:
1972 # * We want conjugate-symmetry, not just symmetry.
1973 # * The diagonal will (as a result) be real.
1976 Eij
= matrix
.zero(F
,n
)
1978 for j
in range(i
+1):
1982 Sij
= cls
.real_embed(Eij
)
1985 # The second one has a minus because it's conjugated.
1986 Eij
[j
,i
] = 1 # Eij = Eij + Eij.transpose()
1987 Sij_real
= cls
.real_embed(Eij
)
1989 # Eij = I*Eij - I*Eij.transpose()
1992 Sij_imag
= cls
.real_embed(Eij
)
1998 # Since we embedded these, we can drop back to the "field" that we
1999 # started with instead of the complex extension "F".
2000 return tuple( s
.change_ring(field
) for s
in S
)
2003 def __init__(self
, n
, **kwargs
):
2004 # We know this is a valid EJA, but will double-check
2005 # if the user passes check_axioms=True.
2006 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2008 super(ComplexHermitianEJA
, self
).__init
__(self
._denormalized
_basis
(n
),
2009 self
.jordan_product
,
2010 self
.trace_inner_product
,
2012 # TODO: this could be factored out somehow, but is left here
2013 # because the MatrixEJA is not presently a subclass of the
2014 # FDEJA class that defines rank() and one().
2015 self
.rank
.set_cache(n
)
2016 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
2017 self
.one
.set_cache(self(idV
))
2020 def _max_random_instance_size():
2021 return 3 # Dimension 9
2024 def random_instance(cls
, **kwargs
):
2026 Return a random instance of this type of algebra.
2028 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2029 return cls(n
, **kwargs
)
2031 class QuaternionMatrixEJA(MatrixEJA
):
2033 # A manual dictionary-cache for the quaternion_extension() method,
2034 # since apparently @classmethods can't also be @cached_methods.
2035 _quaternion_extension
= {}
2038 def quaternion_extension(cls
,field
):
2040 The quaternion field that we embed/unembed, as an extension
2041 of the given ``field``.
2043 if field
in cls
._quaternion
_extension
:
2044 return cls
._quaternion
_extension
[field
]
2046 Q
= QuaternionAlgebra(field
,-1,-1)
2048 cls
._quaternion
_extension
[field
] = Q
2052 def dimension_over_reals():
2056 def real_embed(cls
,M
):
2058 Embed the n-by-n quaternion matrix ``M`` into the space of real
2059 matrices of size 4n-by-4n by first sending each quaternion entry `z
2060 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
2061 c+di],[-c + di, a-bi]]`, and then embedding those into a real
2066 sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA
2070 sage: Q = QuaternionAlgebra(QQ,-1,-1)
2071 sage: i,j,k = Q.gens()
2072 sage: x = 1 + 2*i + 3*j + 4*k
2073 sage: M = matrix(Q, 1, [[x]])
2074 sage: QuaternionMatrixEJA.real_embed(M)
2080 Embedding is a homomorphism (isomorphism, in fact)::
2082 sage: set_random_seed()
2083 sage: n = ZZ.random_element(2)
2084 sage: Q = QuaternionAlgebra(QQ,-1,-1)
2085 sage: X = random_matrix(Q, n)
2086 sage: Y = random_matrix(Q, n)
2087 sage: Xe = QuaternionMatrixEJA.real_embed(X)
2088 sage: Ye = QuaternionMatrixEJA.real_embed(Y)
2089 sage: XYe = QuaternionMatrixEJA.real_embed(X*Y)
2094 super(QuaternionMatrixEJA
,cls
).real_embed(M
)
2095 quaternions
= M
.base_ring()
2098 F
= QuadraticField(-1, 'I')
2103 t
= z
.coefficient_tuple()
2108 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
2109 [-c
+ d
*i
, a
- b
*i
]])
2110 realM
= ComplexMatrixEJA
.real_embed(cplxM
)
2111 blocks
.append(realM
)
2113 # We should have real entries by now, so use the realest field
2114 # we've got for the return value.
2115 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
2120 def real_unembed(cls
,M
):
2122 The inverse of _embed_quaternion_matrix().
2126 sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA
2130 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
2131 ....: [-2, 1, -4, 3],
2132 ....: [-3, 4, 1, -2],
2133 ....: [-4, -3, 2, 1]])
2134 sage: QuaternionMatrixEJA.real_unembed(M)
2135 [1 + 2*i + 3*j + 4*k]
2139 Unembedding is the inverse of embedding::
2141 sage: set_random_seed()
2142 sage: Q = QuaternionAlgebra(QQ, -1, -1)
2143 sage: M = random_matrix(Q, 3)
2144 sage: Me = QuaternionMatrixEJA.real_embed(M)
2145 sage: QuaternionMatrixEJA.real_unembed(Me) == M
2149 super(QuaternionMatrixEJA
,cls
).real_unembed(M
)
2151 d
= cls
.dimension_over_reals()
2153 # Use the base ring of the matrix to ensure that its entries can be
2154 # multiplied by elements of the quaternion algebra.
2155 Q
= cls
.quaternion_extension(M
.base_ring())
2158 # Go top-left to bottom-right (reading order), converting every
2159 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
2162 for l
in range(n
/d
):
2163 for m
in range(n
/d
):
2164 submat
= ComplexMatrixEJA
.real_unembed(
2165 M
[d
*l
:d
*l
+d
,d
*m
:d
*m
+d
] )
2166 if submat
[0,0] != submat
[1,1].conjugate():
2167 raise ValueError('bad on-diagonal submatrix')
2168 if submat
[0,1] != -submat
[1,0].conjugate():
2169 raise ValueError('bad off-diagonal submatrix')
2170 z
= submat
[0,0].real()
2171 z
+= submat
[0,0].imag()*i
2172 z
+= submat
[0,1].real()*j
2173 z
+= submat
[0,1].imag()*k
2176 return matrix(Q
, n
/d
, elements
)
2179 class QuaternionHermitianEJA(ConcreteEJA
, QuaternionMatrixEJA
):
2181 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
2182 matrices, the usual symmetric Jordan product, and the
2183 real-part-of-trace inner product. It has dimension `2n^2 - n` over
2188 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2192 In theory, our "field" can be any subfield of the reals::
2194 sage: QuaternionHermitianEJA(2, field=RDF)
2195 Euclidean Jordan algebra of dimension 6 over Real Double Field
2196 sage: QuaternionHermitianEJA(2, field=RR)
2197 Euclidean Jordan algebra of dimension 6 over Real Field with
2198 53 bits of precision
2202 The dimension of this algebra is `2*n^2 - n`::
2204 sage: set_random_seed()
2205 sage: n_max = QuaternionHermitianEJA._max_random_instance_size()
2206 sage: n = ZZ.random_element(1, n_max)
2207 sage: J = QuaternionHermitianEJA(n)
2208 sage: J.dimension() == 2*(n^2) - n
2211 The Jordan multiplication is what we think it is::
2213 sage: set_random_seed()
2214 sage: J = QuaternionHermitianEJA.random_instance()
2215 sage: x,y = J.random_elements(2)
2216 sage: actual = (x*y).to_matrix()
2217 sage: X = x.to_matrix()
2218 sage: Y = y.to_matrix()
2219 sage: expected = (X*Y + Y*X)/2
2220 sage: actual == expected
2222 sage: J(expected) == x*y
2225 We can change the generator prefix::
2227 sage: QuaternionHermitianEJA(2, prefix='a').gens()
2228 (a0, a1, a2, a3, a4, a5)
2230 We can construct the (trivial) algebra of rank zero::
2232 sage: QuaternionHermitianEJA(0)
2233 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2237 def _denormalized_basis(cls
, n
):
2239 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
2241 Why do we embed these? Basically, because all of numerical
2242 linear algebra assumes that you're working with vectors consisting
2243 of `n` entries from a field and scalars from the same field. There's
2244 no way to tell SageMath that (for example) the vectors contain
2245 complex numbers, while the scalar field is real.
2249 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2253 sage: set_random_seed()
2254 sage: n = ZZ.random_element(1,5)
2255 sage: B = QuaternionHermitianEJA._denormalized_basis(n)
2256 sage: all( M.is_symmetric() for M in B )
2261 Q
= QuaternionAlgebra(QQ
,-1,-1)
2264 # This is like the symmetric case, but we need to be careful:
2266 # * We want conjugate-symmetry, not just symmetry.
2267 # * The diagonal will (as a result) be real.
2270 Eij
= matrix
.zero(Q
,n
)
2272 for j
in range(i
+1):
2276 Sij
= cls
.real_embed(Eij
)
2279 # The second, third, and fourth ones have a minus
2280 # because they're conjugated.
2281 # Eij = Eij + Eij.transpose()
2283 Sij_real
= cls
.real_embed(Eij
)
2285 # Eij = I*(Eij - Eij.transpose())
2288 Sij_I
= cls
.real_embed(Eij
)
2290 # Eij = J*(Eij - Eij.transpose())
2293 Sij_J
= cls
.real_embed(Eij
)
2295 # Eij = K*(Eij - Eij.transpose())
2298 Sij_K
= cls
.real_embed(Eij
)
2304 # Since we embedded these, we can drop back to the "field" that we
2305 # started with instead of the quaternion algebra "Q".
2306 return tuple( s
.change_ring(field
) for s
in S
)
2309 def __init__(self
, n
, **kwargs
):
2310 # We know this is a valid EJA, but will double-check
2311 # if the user passes check_axioms=True.
2312 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2314 super(QuaternionHermitianEJA
, self
).__init
__(self
._denormalized
_basis
(n
),
2315 self
.jordan_product
,
2316 self
.trace_inner_product
,
2318 # TODO: this could be factored out somehow, but is left here
2319 # because the MatrixEJA is not presently a subclass of the
2320 # FDEJA class that defines rank() and one().
2321 self
.rank
.set_cache(n
)
2322 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
2323 self
.one
.set_cache(self(idV
))
2327 def _max_random_instance_size():
2329 The maximum rank of a random QuaternionHermitianEJA.
2331 return 2 # Dimension 6
2334 def random_instance(cls
, **kwargs
):
2336 Return a random instance of this type of algebra.
2338 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2339 return cls(n
, **kwargs
)
2342 class HadamardEJA(ConcreteEJA
):
2344 Return the Euclidean Jordan Algebra corresponding to the set
2345 `R^n` under the Hadamard product.
2347 Note: this is nothing more than the Cartesian product of ``n``
2348 copies of the spin algebra. Once Cartesian product algebras
2349 are implemented, this can go.
2353 sage: from mjo.eja.eja_algebra import HadamardEJA
2357 This multiplication table can be verified by hand::
2359 sage: J = HadamardEJA(3)
2360 sage: e0,e1,e2 = J.gens()
2376 We can change the generator prefix::
2378 sage: HadamardEJA(3, prefix='r').gens()
2382 def __init__(self
, n
, **kwargs
):
2384 jordan_product
= lambda x
,y
: x
2385 inner_product
= lambda x
,y
: x
2387 def jordan_product(x
,y
):
2389 return P( xi
*yi
for (xi
,yi
) in zip(x
,y
) )
2391 def inner_product(x
,y
):
2394 # New defaults for keyword arguments. Don't orthonormalize
2395 # because our basis is already orthonormal with respect to our
2396 # inner-product. Don't check the axioms, because we know this
2397 # is a valid EJA... but do double-check if the user passes
2398 # check_axioms=True. Note: we DON'T override the "check_field"
2399 # default here, because the user can pass in a field!
2400 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2401 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2403 column_basis
= tuple( b
.column() for b
in FreeModule(ZZ
, n
).basis() )
2404 super().__init
__(column_basis
,
2409 self
.rank
.set_cache(n
)
2412 self
.one
.set_cache( self
.zero() )
2414 self
.one
.set_cache( sum(self
.gens()) )
2417 def _max_random_instance_size():
2419 The maximum dimension of a random HadamardEJA.
2424 def random_instance(cls
, **kwargs
):
2426 Return a random instance of this type of algebra.
2428 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2429 return cls(n
, **kwargs
)
2432 class BilinearFormEJA(ConcreteEJA
):
2434 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2435 with the half-trace inner product and jordan product ``x*y =
2436 (<Bx,y>,y_bar>, x0*y_bar + y0*x_bar)`` where `B = 1 \times B22` is
2437 a symmetric positive-definite "bilinear form" matrix. Its
2438 dimension is the size of `B`, and it has rank two in dimensions
2439 larger than two. It reduces to the ``JordanSpinEJA`` when `B` is
2440 the identity matrix of order ``n``.
2442 We insist that the one-by-one upper-left identity block of `B` be
2443 passed in as well so that we can be passed a matrix of size zero
2444 to construct a trivial algebra.
2448 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
2449 ....: JordanSpinEJA)
2453 When no bilinear form is specified, the identity matrix is used,
2454 and the resulting algebra is the Jordan spin algebra::
2456 sage: B = matrix.identity(AA,3)
2457 sage: J0 = BilinearFormEJA(B)
2458 sage: J1 = JordanSpinEJA(3)
2459 sage: J0.multiplication_table() == J0.multiplication_table()
2462 An error is raised if the matrix `B` does not correspond to a
2463 positive-definite bilinear form::
2465 sage: B = matrix.random(QQ,2,3)
2466 sage: J = BilinearFormEJA(B)
2467 Traceback (most recent call last):
2469 ValueError: bilinear form is not positive-definite
2470 sage: B = matrix.zero(QQ,3)
2471 sage: J = BilinearFormEJA(B)
2472 Traceback (most recent call last):
2474 ValueError: bilinear form is not positive-definite
2478 We can create a zero-dimensional algebra::
2480 sage: B = matrix.identity(AA,0)
2481 sage: J = BilinearFormEJA(B)
2485 We can check the multiplication condition given in the Jordan, von
2486 Neumann, and Wigner paper (and also discussed on my "On the
2487 symmetry..." paper). Note that this relies heavily on the standard
2488 choice of basis, as does anything utilizing the bilinear form
2489 matrix. We opt not to orthonormalize the basis, because if we
2490 did, we would have to normalize the `s_{i}` in a similar manner::
2492 sage: set_random_seed()
2493 sage: n = ZZ.random_element(5)
2494 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2495 sage: B11 = matrix.identity(QQ,1)
2496 sage: B22 = M.transpose()*M
2497 sage: B = block_matrix(2,2,[ [B11,0 ],
2499 sage: J = BilinearFormEJA(B, orthonormalize=False)
2500 sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
2501 sage: V = J.vector_space()
2502 sage: sis = [ J( V([0] + (M.inverse()*ei).list()).column() )
2503 ....: for ei in eis ]
2504 sage: actual = [ sis[i]*sis[j]
2505 ....: for i in range(n-1)
2506 ....: for j in range(n-1) ]
2507 sage: expected = [ J.one() if i == j else J.zero()
2508 ....: for i in range(n-1)
2509 ....: for j in range(n-1) ]
2510 sage: actual == expected
2514 def __init__(self
, B
, **kwargs
):
2515 # The matrix "B" is supplied by the user in most cases,
2516 # so it makes sense to check whether or not its positive-
2517 # definite unless we are specifically asked not to...
2518 if ("check_axioms" not in kwargs
) or kwargs
["check_axioms"]:
2519 if not B
.is_positive_definite():
2520 raise ValueError("bilinear form is not positive-definite")
2522 # However, all of the other data for this EJA is computed
2523 # by us in manner that guarantees the axioms are
2524 # satisfied. So, again, unless we are specifically asked to
2525 # verify things, we'll skip the rest of the checks.
2526 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2528 def inner_product(x
,y
):
2529 return (y
.T
*B
*x
)[0,0]
2531 def jordan_product(x
,y
):
2537 z0
= inner_product(y
,x
)
2538 zbar
= y0
*xbar
+ x0
*ybar
2539 return P([z0
] + zbar
.list())
2542 column_basis
= tuple( b
.column() for b
in FreeModule(ZZ
, n
).basis() )
2543 super(BilinearFormEJA
, self
).__init
__(column_basis
,
2548 # The rank of this algebra is two, unless we're in a
2549 # one-dimensional ambient space (because the rank is bounded
2550 # by the ambient dimension).
2551 self
.rank
.set_cache(min(n
,2))
2554 self
.one
.set_cache( self
.zero() )
2556 self
.one
.set_cache( self
.monomial(0) )
2559 def _max_random_instance_size():
2561 The maximum dimension of a random BilinearFormEJA.
2566 def random_instance(cls
, **kwargs
):
2568 Return a random instance of this algebra.
2570 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2572 B
= matrix
.identity(ZZ
, n
)
2573 return cls(B
, **kwargs
)
2575 B11
= matrix
.identity(ZZ
, 1)
2576 M
= matrix
.random(ZZ
, n
-1)
2577 I
= matrix
.identity(ZZ
, n
-1)
2579 while alpha
.is_zero():
2580 alpha
= ZZ
.random_element().abs()
2581 B22
= M
.transpose()*M
+ alpha
*I
2583 from sage
.matrix
.special
import block_matrix
2584 B
= block_matrix(2,2, [ [B11
, ZZ(0) ],
2587 return cls(B
, **kwargs
)
2590 class JordanSpinEJA(BilinearFormEJA
):
2592 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2593 with the usual inner product and jordan product ``x*y =
2594 (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2599 sage: from mjo.eja.eja_algebra import JordanSpinEJA
2603 This multiplication table can be verified by hand::
2605 sage: J = JordanSpinEJA(4)
2606 sage: e0,e1,e2,e3 = J.gens()
2622 We can change the generator prefix::
2624 sage: JordanSpinEJA(2, prefix='B').gens()
2629 Ensure that we have the usual inner product on `R^n`::
2631 sage: set_random_seed()
2632 sage: J = JordanSpinEJA.random_instance()
2633 sage: x,y = J.random_elements(2)
2634 sage: actual = x.inner_product(y)
2635 sage: expected = x.to_vector().inner_product(y.to_vector())
2636 sage: actual == expected
2640 def __init__(self
, n
, **kwargs
):
2641 # This is a special case of the BilinearFormEJA with the
2642 # identity matrix as its bilinear form.
2643 B
= matrix
.identity(ZZ
, n
)
2645 # Don't orthonormalize because our basis is already
2646 # orthonormal with respect to our inner-product.
2647 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2649 # But also don't pass check_field=False here, because the user
2650 # can pass in a field!
2651 super(JordanSpinEJA
, self
).__init
__(B
, **kwargs
)
2654 def _max_random_instance_size():
2656 The maximum dimension of a random JordanSpinEJA.
2661 def random_instance(cls
, **kwargs
):
2663 Return a random instance of this type of algebra.
2665 Needed here to override the implementation for ``BilinearFormEJA``.
2667 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2668 return cls(n
, **kwargs
)
2671 class TrivialEJA(ConcreteEJA
):
2673 The trivial Euclidean Jordan algebra consisting of only a zero element.
2677 sage: from mjo.eja.eja_algebra import TrivialEJA
2681 sage: J = TrivialEJA()
2688 sage: 7*J.one()*12*J.one()
2690 sage: J.one().inner_product(J.one())
2692 sage: J.one().norm()
2694 sage: J.one().subalgebra_generated_by()
2695 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2700 def __init__(self
, **kwargs
):
2701 jordan_product
= lambda x
,y
: x
2702 inner_product
= lambda x
,y
: 0
2705 # New defaults for keyword arguments
2706 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2707 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2709 super(TrivialEJA
, self
).__init
__(basis
,
2713 # The rank is zero using my definition, namely the dimension of the
2714 # largest subalgebra generated by any element.
2715 self
.rank
.set_cache(0)
2716 self
.one
.set_cache( self
.zero() )
2719 def random_instance(cls
, **kwargs
):
2720 # We don't take a "size" argument so the superclass method is
2721 # inappropriate for us.
2722 return cls(**kwargs
)
2725 class CartesianProductEJA(CombinatorialFreeModule_CartesianProduct
,
2726 FiniteDimensionalEJA
):
2728 The external (orthogonal) direct sum of two or more Euclidean
2729 Jordan algebras. Every Euclidean Jordan algebra decomposes into an
2730 orthogonal direct sum of simple Euclidean Jordan algebras which is
2731 then isometric to a Cartesian product, so no generality is lost by
2732 providing only this construction.
2736 sage: from mjo.eja.eja_algebra import (random_eja,
2737 ....: CartesianProductEJA,
2739 ....: JordanSpinEJA,
2740 ....: RealSymmetricEJA)
2744 The Jordan product is inherited from our factors and implemented by
2745 our CombinatorialFreeModule Cartesian product superclass::
2747 sage: set_random_seed()
2748 sage: J1 = HadamardEJA(2)
2749 sage: J2 = RealSymmetricEJA(2)
2750 sage: J = cartesian_product([J1,J2])
2751 sage: x,y = J.random_elements(2)
2755 The ability to retrieve the original factors is implemented by our
2756 CombinatorialFreeModule Cartesian product superclass::
2758 sage: J1 = HadamardEJA(2, field=QQ)
2759 sage: J2 = JordanSpinEJA(3, field=QQ)
2760 sage: J = cartesian_product([J1,J2])
2761 sage: J.cartesian_factors()
2762 (Euclidean Jordan algebra of dimension 2 over Rational Field,
2763 Euclidean Jordan algebra of dimension 3 over Rational Field)
2765 You can provide more than two factors::
2767 sage: J1 = HadamardEJA(2)
2768 sage: J2 = JordanSpinEJA(3)
2769 sage: J3 = RealSymmetricEJA(3)
2770 sage: cartesian_product([J1,J2,J3])
2771 Euclidean Jordan algebra of dimension 2 over Algebraic Real
2772 Field (+) Euclidean Jordan algebra of dimension 3 over Algebraic
2773 Real Field (+) Euclidean Jordan algebra of dimension 6 over
2774 Algebraic Real Field
2776 Rank is additive on a Cartesian product::
2778 sage: J1 = HadamardEJA(1)
2779 sage: J2 = RealSymmetricEJA(2)
2780 sage: J = cartesian_product([J1,J2])
2781 sage: J1.rank.clear_cache()
2782 sage: J2.rank.clear_cache()
2783 sage: J.rank.clear_cache()
2786 sage: J.rank() == J1.rank() + J2.rank()
2789 The same rank computation works over the rationals, with whatever
2792 sage: J1 = HadamardEJA(1, field=QQ, orthonormalize=False)
2793 sage: J2 = RealSymmetricEJA(2, field=QQ, orthonormalize=False)
2794 sage: J = cartesian_product([J1,J2])
2795 sage: J1.rank.clear_cache()
2796 sage: J2.rank.clear_cache()
2797 sage: J.rank.clear_cache()
2800 sage: J.rank() == J1.rank() + J2.rank()
2803 The product algebra will be associative if and only if all of its
2804 components are associative::
2806 sage: J1 = HadamardEJA(2)
2807 sage: J1.is_associative()
2809 sage: J2 = HadamardEJA(3)
2810 sage: J2.is_associative()
2812 sage: J3 = RealSymmetricEJA(3)
2813 sage: J3.is_associative()
2815 sage: CP1 = cartesian_product([J1,J2])
2816 sage: CP1.is_associative()
2818 sage: CP2 = cartesian_product([J1,J3])
2819 sage: CP2.is_associative()
2824 All factors must share the same base field::
2826 sage: J1 = HadamardEJA(2, field=QQ)
2827 sage: J2 = RealSymmetricEJA(2)
2828 sage: CartesianProductEJA((J1,J2))
2829 Traceback (most recent call last):
2831 ValueError: all factors must share the same base field
2833 The cached unit element is the same one that would be computed::
2835 sage: set_random_seed() # long time
2836 sage: J1 = random_eja() # long time
2837 sage: J2 = random_eja() # long time
2838 sage: J = cartesian_product([J1,J2]) # long time
2839 sage: actual = J.one() # long time
2840 sage: J.one.clear_cache() # long time
2841 sage: expected = J.one() # long time
2842 sage: actual == expected # long time
2846 def __init__(self
, algebras
, **kwargs
):
2847 CombinatorialFreeModule_CartesianProduct
.__init
__(self
,
2850 field
= algebras
[0].base_ring()
2851 if not all( J
.base_ring() == field
for J
in algebras
):
2852 raise ValueError("all factors must share the same base field")
2854 associative
= all( m
.is_associative() for m
in algebras
)
2856 # The definition of matrix_space() and self.basis() relies
2857 # only on the stuff in the CFM_CartesianProduct class, which
2858 # we've already initialized.
2859 Js
= self
.cartesian_factors()
2861 MS
= self
.matrix_space()
2863 MS(tuple( self
.cartesian_projection(i
)(b
).to_matrix()
2864 for i
in range(m
) ))
2865 for b
in self
.basis()
2868 # Define jordan/inner products that operate on that matrix_basis.
2869 def jordan_product(x
,y
):
2871 (Js
[i
](x
[i
])*Js
[i
](y
[i
])).to_matrix() for i
in range(m
)
2874 def inner_product(x
, y
):
2876 Js
[i
](x
[i
]).inner_product(Js
[i
](y
[i
])) for i
in range(m
)
2879 # There's no need to check the field since it already came
2880 # from an EJA. Likewise the axioms are guaranteed to be
2881 # satisfied, unless the guy writing this class sucks.
2883 # If you want the basis to be orthonormalized, orthonormalize
2885 FiniteDimensionalEJA
.__init
__(self
,
2890 orthonormalize
=False,
2891 associative
=associative
,
2892 cartesian_product
=True,
2896 ones
= tuple(J
.one() for J
in algebras
)
2897 self
.one
.set_cache(self
._cartesian
_product
_of
_elements
(ones
))
2898 self
.rank
.set_cache(sum(J
.rank() for J
in algebras
))
2900 def matrix_space(self
):
2902 Return the space that our matrix basis lives in as a Cartesian
2907 sage: from mjo.eja.eja_algebra import (HadamardEJA,
2908 ....: RealSymmetricEJA)
2912 sage: J1 = HadamardEJA(1)
2913 sage: J2 = RealSymmetricEJA(2)
2914 sage: J = cartesian_product([J1,J2])
2915 sage: J.matrix_space()
2916 The Cartesian product of (Full MatrixSpace of 1 by 1 dense
2917 matrices over Algebraic Real Field, Full MatrixSpace of 2
2918 by 2 dense matrices over Algebraic Real Field)
2921 from sage
.categories
.cartesian_product
import cartesian_product
2922 return cartesian_product( [J
.matrix_space()
2923 for J
in self
.cartesian_factors()] )
2926 def cartesian_projection(self
, i
):
2930 sage: from mjo.eja.eja_algebra import (random_eja,
2931 ....: JordanSpinEJA,
2933 ....: RealSymmetricEJA,
2934 ....: ComplexHermitianEJA)
2938 The projection morphisms are Euclidean Jordan algebra
2941 sage: J1 = HadamardEJA(2)
2942 sage: J2 = RealSymmetricEJA(2)
2943 sage: J = cartesian_product([J1,J2])
2944 sage: J.cartesian_projection(0)
2945 Linear operator between finite-dimensional Euclidean Jordan
2946 algebras represented by the matrix:
2949 Domain: Euclidean Jordan algebra of dimension 2 over Algebraic
2950 Real Field (+) Euclidean Jordan algebra of dimension 3 over
2951 Algebraic Real Field
2952 Codomain: Euclidean Jordan algebra of dimension 2 over Algebraic
2954 sage: J.cartesian_projection(1)
2955 Linear operator between finite-dimensional Euclidean Jordan
2956 algebras represented by the matrix:
2960 Domain: Euclidean Jordan algebra of dimension 2 over Algebraic
2961 Real Field (+) Euclidean Jordan algebra of dimension 3 over
2962 Algebraic Real Field
2963 Codomain: Euclidean Jordan algebra of dimension 3 over Algebraic
2966 The projections work the way you'd expect on the vector
2967 representation of an element::
2969 sage: J1 = JordanSpinEJA(2)
2970 sage: J2 = ComplexHermitianEJA(2)
2971 sage: J = cartesian_product([J1,J2])
2972 sage: pi_left = J.cartesian_projection(0)
2973 sage: pi_right = J.cartesian_projection(1)
2974 sage: pi_left(J.one()).to_vector()
2976 sage: pi_right(J.one()).to_vector()
2978 sage: J.one().to_vector()
2983 The answer never changes::
2985 sage: set_random_seed()
2986 sage: J1 = random_eja()
2987 sage: J2 = random_eja()
2988 sage: J = cartesian_product([J1,J2])
2989 sage: P0 = J.cartesian_projection(0)
2990 sage: P1 = J.cartesian_projection(0)
2995 Ji
= self
.cartesian_factors()[i
]
2996 # Requires the fix on Trac 31421/31422 to work!
2997 Pi
= super().cartesian_projection(i
)
2998 return FiniteDimensionalEJAOperator(self
,Ji
,Pi
.matrix())
3001 def cartesian_embedding(self
, i
):
3005 sage: from mjo.eja.eja_algebra import (random_eja,
3006 ....: JordanSpinEJA,
3008 ....: RealSymmetricEJA)
3012 The embedding morphisms are Euclidean Jordan algebra
3015 sage: J1 = HadamardEJA(2)
3016 sage: J2 = RealSymmetricEJA(2)
3017 sage: J = cartesian_product([J1,J2])
3018 sage: J.cartesian_embedding(0)
3019 Linear operator between finite-dimensional Euclidean Jordan
3020 algebras represented by the matrix:
3026 Domain: Euclidean Jordan algebra of dimension 2 over
3027 Algebraic Real Field
3028 Codomain: Euclidean Jordan algebra of dimension 2 over
3029 Algebraic Real Field (+) Euclidean Jordan algebra of
3030 dimension 3 over Algebraic Real Field
3031 sage: J.cartesian_embedding(1)
3032 Linear operator between finite-dimensional Euclidean Jordan
3033 algebras represented by the matrix:
3039 Domain: Euclidean Jordan algebra of dimension 3 over
3040 Algebraic Real Field
3041 Codomain: Euclidean Jordan algebra of dimension 2 over
3042 Algebraic Real Field (+) Euclidean Jordan algebra of
3043 dimension 3 over Algebraic Real Field
3045 The embeddings work the way you'd expect on the vector
3046 representation of an element::
3048 sage: J1 = JordanSpinEJA(3)
3049 sage: J2 = RealSymmetricEJA(2)
3050 sage: J = cartesian_product([J1,J2])
3051 sage: iota_left = J.cartesian_embedding(0)
3052 sage: iota_right = J.cartesian_embedding(1)
3053 sage: iota_left(J1.zero()) == J.zero()
3055 sage: iota_right(J2.zero()) == J.zero()
3057 sage: J1.one().to_vector()
3059 sage: iota_left(J1.one()).to_vector()
3061 sage: J2.one().to_vector()
3063 sage: iota_right(J2.one()).to_vector()
3065 sage: J.one().to_vector()
3070 The answer never changes::
3072 sage: set_random_seed()
3073 sage: J1 = random_eja()
3074 sage: J2 = random_eja()
3075 sage: J = cartesian_product([J1,J2])
3076 sage: E0 = J.cartesian_embedding(0)
3077 sage: E1 = J.cartesian_embedding(0)
3081 Composing a projection with the corresponding inclusion should
3082 produce the identity map, and mismatching them should produce
3085 sage: set_random_seed()
3086 sage: J1 = random_eja()
3087 sage: J2 = random_eja()
3088 sage: J = cartesian_product([J1,J2])
3089 sage: iota_left = J.cartesian_embedding(0)
3090 sage: iota_right = J.cartesian_embedding(1)
3091 sage: pi_left = J.cartesian_projection(0)
3092 sage: pi_right = J.cartesian_projection(1)
3093 sage: pi_left*iota_left == J1.one().operator()
3095 sage: pi_right*iota_right == J2.one().operator()
3097 sage: (pi_left*iota_right).is_zero()
3099 sage: (pi_right*iota_left).is_zero()
3103 Ji
= self
.cartesian_factors()[i
]
3104 # Requires the fix on Trac 31421/31422 to work!
3105 Ei
= super().cartesian_embedding(i
)
3106 return FiniteDimensionalEJAOperator(Ji
,self
,Ei
.matrix())
3109 def _element_constructor_(self
, elt
):
3111 Construct an element of this algebra from an ordered tuple.
3113 We just apply the element constructor from each of our factors
3114 to the corresponding component of the tuple, and package up
3119 sage: from mjo.eja.eja_algebra import (HadamardEJA,
3120 ....: RealSymmetricEJA)
3124 sage: J1 = HadamardEJA(3)
3125 sage: J2 = RealSymmetricEJA(2)
3126 sage: J = cartesian_product([J1,J2])
3127 sage: J( (J1.matrix_basis()[1], J2.matrix_basis()[2]) )
3130 m
= len(self
.cartesian_factors())
3132 z
= tuple( self
.cartesian_factors()[i
](elt
[i
]) for i
in range(m
) )
3133 return self
._cartesian
_product
_of
_elements
(z
)
3135 raise ValueError("not an element of this algebra")
3137 Element
= CartesianProductEJAElement
3140 FiniteDimensionalEJA
.CartesianProduct
= CartesianProductEJA
3142 random_eja
= ConcreteEJA
.random_instance
3143 #def random_eja(*args, **kwargs):
3144 # from sage.categories.cartesian_product import cartesian_product
3145 # J1 = HadamardEJA(1, **kwargs)
3146 # J2 = RealSymmetricEJA(2, **kwargs)
3147 # J = cartesian_product([J1,J2])