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_commutative(self
):
337 Whether or not this algebra's multiplication table is commutative.
339 This method should of course always return ``True``, unless
340 this algebra was constructed with ``check_axioms=False`` and
341 passed an invalid multiplication table.
343 return all( self
.product_on_basis(i
,j
) == self
.product_on_basis(i
,j
)
344 for i
in range(self
.dimension())
345 for j
in range(self
.dimension()) )
347 def _is_jordanian(self
):
349 Whether or not this algebra's multiplication table respects the
350 Jordan identity `(x^{2})(xy) = x(x^{2}y)`.
352 We only check one arrangement of `x` and `y`, so for a
353 ``True`` result to be truly true, you should also check
354 :meth:`_is_commutative`. This method should of course always
355 return ``True``, unless this algebra was constructed with
356 ``check_axioms=False`` and passed an invalid multiplication table.
358 return all( (self
.gens()[i
]**2)*(self
.gens()[i
]*self
.gens()[j
])
360 (self
.gens()[i
])*((self
.gens()[i
]**2)*self
.gens()[j
])
361 for i
in range(self
.dimension())
362 for j
in range(self
.dimension()) )
364 def _inner_product_is_associative(self
):
366 Return whether or not this algebra's inner product `B` is
367 associative; that is, whether or not `B(xy,z) = B(x,yz)`.
369 This method should of course always return ``True``, unless
370 this algebra was constructed with ``check_axioms=False`` and
371 passed an invalid Jordan or inner-product.
374 # Used to check whether or not something is zero in an inexact
375 # ring. This number is sufficient to allow the construction of
376 # QuaternionHermitianEJA(2, field=RDF) with check_axioms=True.
379 for i
in range(self
.dimension()):
380 for j
in range(self
.dimension()):
381 for k
in range(self
.dimension()):
385 diff
= (x
*y
).inner_product(z
) - x
.inner_product(y
*z
)
387 if self
.base_ring().is_exact():
391 if diff
.abs() > epsilon
:
396 def _element_constructor_(self
, elt
):
398 Construct an element of this algebra from its vector or matrix
401 This gets called only after the parent element _call_ method
402 fails to find a coercion for the argument.
406 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
408 ....: RealSymmetricEJA)
412 The identity in `S^n` is converted to the identity in the EJA::
414 sage: J = RealSymmetricEJA(3)
415 sage: I = matrix.identity(QQ,3)
416 sage: J(I) == J.one()
419 This skew-symmetric matrix can't be represented in the EJA::
421 sage: J = RealSymmetricEJA(3)
422 sage: A = matrix(QQ,3, lambda i,j: i-j)
424 Traceback (most recent call last):
426 ValueError: not an element of this algebra
430 Ensure that we can convert any element of the two non-matrix
431 simple algebras (whose matrix representations are columns)
432 back and forth faithfully::
434 sage: set_random_seed()
435 sage: J = HadamardEJA.random_instance()
436 sage: x = J.random_element()
437 sage: J(x.to_vector().column()) == x
439 sage: J = JordanSpinEJA.random_instance()
440 sage: x = J.random_element()
441 sage: J(x.to_vector().column()) == x
445 msg
= "not an element of this algebra"
447 # The superclass implementation of random_element()
448 # needs to be able to coerce "0" into the algebra.
450 elif elt
in self
.base_ring():
451 # Ensure that no base ring -> algebra coercion is performed
452 # by this method. There's some stupidity in sage that would
453 # otherwise propagate to this method; for example, sage thinks
454 # that the integer 3 belongs to the space of 2-by-2 matrices.
455 raise ValueError(msg
)
459 except (AttributeError, TypeError):
460 # Try to convert a vector into a column-matrix
463 if elt
not in self
.matrix_space():
464 raise ValueError(msg
)
466 # Thanks for nothing! Matrix spaces aren't vector spaces in
467 # Sage, so we have to figure out its matrix-basis coordinates
468 # ourselves. We use the basis space's ring instead of the
469 # element's ring because the basis space might be an algebraic
470 # closure whereas the base ring of the 3-by-3 identity matrix
471 # could be QQ instead of QQbar.
473 # We pass check=False because the matrix basis is "guaranteed"
474 # to be linearly independent... right? Ha ha.
475 V
= VectorSpace(self
.base_ring(), elt
.nrows()*elt
.ncols())
476 W
= V
.span_of_basis( (_mat2vec(s
) for s
in self
.matrix_basis()),
480 coords
= W
.coordinate_vector(_mat2vec(elt
))
481 except ArithmeticError: # vector is not in free module
482 raise ValueError(msg
)
484 return self
.from_vector(coords
)
488 Return a string representation of ``self``.
492 sage: from mjo.eja.eja_algebra import JordanSpinEJA
496 Ensure that it says what we think it says::
498 sage: JordanSpinEJA(2, field=AA)
499 Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
500 sage: JordanSpinEJA(3, field=RDF)
501 Euclidean Jordan algebra of dimension 3 over Real Double Field
504 fmt
= "Euclidean Jordan algebra of dimension {} over {}"
505 return fmt
.format(self
.dimension(), self
.base_ring())
509 def characteristic_polynomial_of(self
):
511 Return the algebra's "characteristic polynomial of" function,
512 which is itself a multivariate polynomial that, when evaluated
513 at the coordinates of some algebra element, returns that
514 element's characteristic polynomial.
516 The resulting polynomial has `n+1` variables, where `n` is the
517 dimension of this algebra. The first `n` variables correspond to
518 the coordinates of an algebra element: when evaluated at the
519 coordinates of an algebra element with respect to a certain
520 basis, the result is a univariate polynomial (in the one
521 remaining variable ``t``), namely the characteristic polynomial
526 sage: from mjo.eja.eja_algebra import JordanSpinEJA, TrivialEJA
530 The characteristic polynomial in the spin algebra is given in
531 Alizadeh, Example 11.11::
533 sage: J = JordanSpinEJA(3)
534 sage: p = J.characteristic_polynomial_of(); p
535 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
536 sage: xvec = J.one().to_vector()
540 By definition, the characteristic polynomial is a monic
541 degree-zero polynomial in a rank-zero algebra. Note that
542 Cayley-Hamilton is indeed satisfied since the polynomial
543 ``1`` evaluates to the identity element of the algebra on
546 sage: J = TrivialEJA()
547 sage: J.characteristic_polynomial_of()
554 # The list of coefficient polynomials a_0, a_1, a_2, ..., a_(r-1).
555 a
= self
._charpoly
_coefficients
()
557 # We go to a bit of trouble here to reorder the
558 # indeterminates, so that it's easier to evaluate the
559 # characteristic polynomial at x's coordinates and get back
560 # something in terms of t, which is what we want.
561 S
= PolynomialRing(self
.base_ring(),'t')
565 S
= PolynomialRing(S
, R
.variable_names())
568 return (t
**r
+ sum( a
[k
]*(t
**k
) for k
in range(r
) ))
570 def coordinate_polynomial_ring(self
):
572 The multivariate polynomial ring in which this algebra's
573 :meth:`characteristic_polynomial_of` lives.
577 sage: from mjo.eja.eja_algebra import (HadamardEJA,
578 ....: RealSymmetricEJA)
582 sage: J = HadamardEJA(2)
583 sage: J.coordinate_polynomial_ring()
584 Multivariate Polynomial Ring in X1, X2...
585 sage: J = RealSymmetricEJA(3,field=QQ,orthonormalize=False)
586 sage: J.coordinate_polynomial_ring()
587 Multivariate Polynomial Ring in X1, X2, X3, X4, X5, X6...
590 var_names
= tuple( "X%d" % z
for z
in range(1, self
.dimension()+1) )
591 return PolynomialRing(self
.base_ring(), var_names
)
593 def inner_product(self
, x
, y
):
595 The inner product associated with this Euclidean Jordan algebra.
597 Defaults to the trace inner product, but can be overridden by
598 subclasses if they are sure that the necessary properties are
603 sage: from mjo.eja.eja_algebra import (random_eja,
605 ....: BilinearFormEJA)
609 Our inner product is "associative," which means the following for
610 a symmetric bilinear form::
612 sage: set_random_seed()
613 sage: J = random_eja()
614 sage: x,y,z = J.random_elements(3)
615 sage: (x*y).inner_product(z) == y.inner_product(x*z)
620 Ensure that this is the usual inner product for the algebras
623 sage: set_random_seed()
624 sage: J = HadamardEJA.random_instance()
625 sage: x,y = J.random_elements(2)
626 sage: actual = x.inner_product(y)
627 sage: expected = x.to_vector().inner_product(y.to_vector())
628 sage: actual == expected
631 Ensure that this is one-half of the trace inner-product in a
632 BilinearFormEJA that isn't just the reals (when ``n`` isn't
633 one). This is in Faraut and Koranyi, and also my "On the
636 sage: set_random_seed()
637 sage: J = BilinearFormEJA.random_instance()
638 sage: n = J.dimension()
639 sage: x = J.random_element()
640 sage: y = J.random_element()
641 sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
644 B
= self
._inner
_product
_matrix
645 return (B
*x
.to_vector()).inner_product(y
.to_vector())
648 def is_trivial(self
):
650 Return whether or not this algebra is trivial.
652 A trivial algebra contains only the zero element.
656 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
661 sage: J = ComplexHermitianEJA(3)
667 sage: J = TrivialEJA()
672 return self
.dimension() == 0
675 def multiplication_table(self
):
677 Return a visual representation of this algebra's multiplication
678 table (on basis elements).
682 sage: from mjo.eja.eja_algebra import JordanSpinEJA
686 sage: J = JordanSpinEJA(4)
687 sage: J.multiplication_table()
688 +----++----+----+----+----+
689 | * || e0 | e1 | e2 | e3 |
690 +====++====+====+====+====+
691 | e0 || e0 | e1 | e2 | e3 |
692 +----++----+----+----+----+
693 | e1 || e1 | e0 | 0 | 0 |
694 +----++----+----+----+----+
695 | e2 || e2 | 0 | e0 | 0 |
696 +----++----+----+----+----+
697 | e3 || e3 | 0 | 0 | e0 |
698 +----++----+----+----+----+
702 # Prepend the header row.
703 M
= [["*"] + list(self
.gens())]
705 # And to each subsequent row, prepend an entry that belongs to
706 # the left-side "header column."
707 M
+= [ [self
.gens()[i
]] + [ self
.product_on_basis(i
,j
)
711 return table(M
, header_row
=True, header_column
=True, frame
=True)
714 def matrix_basis(self
):
716 Return an (often more natural) representation of this algebras
717 basis as an ordered tuple of matrices.
719 Every finite-dimensional Euclidean Jordan Algebra is a, up to
720 Jordan isomorphism, a direct sum of five simple
721 algebras---four of which comprise Hermitian matrices. And the
722 last type of algebra can of course be thought of as `n`-by-`1`
723 column matrices (ambiguusly called column vectors) to avoid
724 special cases. As a result, matrices (and column vectors) are
725 a natural representation format for Euclidean Jordan algebra
728 But, when we construct an algebra from a basis of matrices,
729 those matrix representations are lost in favor of coordinate
730 vectors *with respect to* that basis. We could eventually
731 convert back if we tried hard enough, but having the original
732 representations handy is valuable enough that we simply store
733 them and return them from this method.
735 Why implement this for non-matrix algebras? Avoiding special
736 cases for the :class:`BilinearFormEJA` pays with simplicity in
737 its own right. But mainly, we would like to be able to assume
738 that elements of a :class:`CartesianProductEJA` can be displayed
739 nicely, without having to have special classes for direct sums
740 one of whose components was a matrix algebra.
744 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
745 ....: RealSymmetricEJA)
749 sage: J = RealSymmetricEJA(2)
751 Finite family {0: e0, 1: e1, 2: e2}
752 sage: J.matrix_basis()
754 [1 0] [ 0 0.7071067811865475?] [0 0]
755 [0 0], [0.7071067811865475? 0], [0 1]
760 sage: J = JordanSpinEJA(2)
762 Finite family {0: e0, 1: e1}
763 sage: J.matrix_basis()
769 return self
._matrix
_basis
772 def matrix_space(self
):
774 Return the matrix space in which this algebra's elements live, if
775 we think of them as matrices (including column vectors of the
778 Generally this will be an `n`-by-`1` column-vector space,
779 except when the algebra is trivial. There it's `n`-by-`n`
780 (where `n` is zero), to ensure that two elements of the matrix
781 space (empty matrices) can be multiplied.
783 Matrix algebras override this with something more useful.
785 if self
.is_trivial():
786 return MatrixSpace(self
.base_ring(), 0)
788 return self
.matrix_basis()[0].parent()
794 Return the unit element of this algebra.
798 sage: from mjo.eja.eja_algebra import (HadamardEJA,
803 We can compute unit element in the Hadamard EJA::
805 sage: J = HadamardEJA(5)
807 e0 + e1 + e2 + e3 + e4
809 The unit element in the Hadamard EJA is inherited in the
810 subalgebras generated by its elements::
812 sage: J = HadamardEJA(5)
814 e0 + e1 + e2 + e3 + e4
815 sage: x = sum(J.gens())
816 sage: A = x.subalgebra_generated_by(orthonormalize=False)
819 sage: A.one().superalgebra_element()
820 e0 + e1 + e2 + e3 + e4
824 The identity element acts like the identity, regardless of
825 whether or not we orthonormalize::
827 sage: set_random_seed()
828 sage: J = random_eja()
829 sage: x = J.random_element()
830 sage: J.one()*x == x and x*J.one() == x
832 sage: A = x.subalgebra_generated_by()
833 sage: y = A.random_element()
834 sage: A.one()*y == y and y*A.one() == y
839 sage: set_random_seed()
840 sage: J = random_eja(field=QQ, orthonormalize=False)
841 sage: x = J.random_element()
842 sage: J.one()*x == x and x*J.one() == x
844 sage: A = x.subalgebra_generated_by(orthonormalize=False)
845 sage: y = A.random_element()
846 sage: A.one()*y == y and y*A.one() == y
849 The matrix of the unit element's operator is the identity,
850 regardless of the base field and whether or not we
853 sage: set_random_seed()
854 sage: J = random_eja()
855 sage: actual = J.one().operator().matrix()
856 sage: expected = matrix.identity(J.base_ring(), J.dimension())
857 sage: actual == expected
859 sage: x = J.random_element()
860 sage: A = x.subalgebra_generated_by()
861 sage: actual = A.one().operator().matrix()
862 sage: expected = matrix.identity(A.base_ring(), A.dimension())
863 sage: actual == expected
868 sage: set_random_seed()
869 sage: J = random_eja(field=QQ, orthonormalize=False)
870 sage: actual = J.one().operator().matrix()
871 sage: expected = matrix.identity(J.base_ring(), J.dimension())
872 sage: actual == expected
874 sage: x = J.random_element()
875 sage: A = x.subalgebra_generated_by(orthonormalize=False)
876 sage: actual = A.one().operator().matrix()
877 sage: expected = matrix.identity(A.base_ring(), A.dimension())
878 sage: actual == expected
881 Ensure that the cached unit element (often precomputed by
882 hand) agrees with the computed one::
884 sage: set_random_seed()
885 sage: J = random_eja()
886 sage: cached = J.one()
887 sage: J.one.clear_cache()
888 sage: J.one() == cached
893 sage: set_random_seed()
894 sage: J = random_eja(field=QQ, orthonormalize=False)
895 sage: cached = J.one()
896 sage: J.one.clear_cache()
897 sage: J.one() == cached
901 # We can brute-force compute the matrices of the operators
902 # that correspond to the basis elements of this algebra.
903 # If some linear combination of those basis elements is the
904 # algebra identity, then the same linear combination of
905 # their matrices has to be the identity matrix.
907 # Of course, matrices aren't vectors in sage, so we have to
908 # appeal to the "long vectors" isometry.
909 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
911 # Now we use basic linear algebra to find the coefficients,
912 # of the matrices-as-vectors-linear-combination, which should
913 # work for the original algebra basis too.
914 A
= matrix(self
.base_ring(), oper_vecs
)
916 # We used the isometry on the left-hand side already, but we
917 # still need to do it for the right-hand side. Recall that we
918 # wanted something that summed to the identity matrix.
919 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
921 # Now if there's an identity element in the algebra, this
922 # should work. We solve on the left to avoid having to
923 # transpose the matrix "A".
924 return self
.from_vector(A
.solve_left(b
))
927 def peirce_decomposition(self
, c
):
929 The Peirce decomposition of this algebra relative to the
932 In the future, this can be extended to a complete system of
933 orthogonal idempotents.
937 - ``c`` -- an idempotent of this algebra.
941 A triple (J0, J5, J1) containing two subalgebras and one subspace
944 - ``J0`` -- the algebra on the eigenspace of ``c.operator()``
945 corresponding to the eigenvalue zero.
947 - ``J5`` -- the eigenspace (NOT a subalgebra) of ``c.operator()``
948 corresponding to the eigenvalue one-half.
950 - ``J1`` -- the algebra on the eigenspace of ``c.operator()``
951 corresponding to the eigenvalue one.
953 These are the only possible eigenspaces for that operator, and this
954 algebra is a direct sum of them. The spaces ``J0`` and ``J1`` are
955 orthogonal, and are subalgebras of this algebra with the appropriate
960 sage: from mjo.eja.eja_algebra import random_eja, RealSymmetricEJA
964 The canonical example comes from the symmetric matrices, which
965 decompose into diagonal and off-diagonal parts::
967 sage: J = RealSymmetricEJA(3)
968 sage: C = matrix(QQ, [ [1,0,0],
972 sage: J0,J5,J1 = J.peirce_decomposition(c)
974 Euclidean Jordan algebra of dimension 1...
976 Vector space of degree 6 and dimension 2...
978 Euclidean Jordan algebra of dimension 3...
979 sage: J0.one().to_matrix()
983 sage: orig_df = AA.options.display_format
984 sage: AA.options.display_format = 'radical'
985 sage: J.from_vector(J5.basis()[0]).to_matrix()
989 sage: J.from_vector(J5.basis()[1]).to_matrix()
993 sage: AA.options.display_format = orig_df
994 sage: J1.one().to_matrix()
1001 Every algebra decomposes trivially with respect to its identity
1004 sage: set_random_seed()
1005 sage: J = random_eja()
1006 sage: J0,J5,J1 = J.peirce_decomposition(J.one())
1007 sage: J0.dimension() == 0 and J5.dimension() == 0
1009 sage: J1.superalgebra() == J and J1.dimension() == J.dimension()
1012 The decomposition is into eigenspaces, and its components are
1013 therefore necessarily orthogonal. Moreover, the identity
1014 elements in the two subalgebras are the projections onto their
1015 respective subspaces of the superalgebra's identity element::
1017 sage: set_random_seed()
1018 sage: J = random_eja()
1019 sage: x = J.random_element()
1020 sage: if not J.is_trivial():
1021 ....: while x.is_nilpotent():
1022 ....: x = J.random_element()
1023 sage: c = x.subalgebra_idempotent()
1024 sage: J0,J5,J1 = J.peirce_decomposition(c)
1026 sage: for (w,y,z) in zip(J0.basis(), J5.basis(), J1.basis()):
1027 ....: w = w.superalgebra_element()
1028 ....: y = J.from_vector(y)
1029 ....: z = z.superalgebra_element()
1030 ....: ipsum += w.inner_product(y).abs()
1031 ....: ipsum += w.inner_product(z).abs()
1032 ....: ipsum += y.inner_product(z).abs()
1035 sage: J1(c) == J1.one()
1037 sage: J0(J.one() - c) == J0.one()
1041 if not c
.is_idempotent():
1042 raise ValueError("element is not idempotent: %s" % c
)
1044 from mjo
.eja
.eja_subalgebra
import FiniteDimensionalEJASubalgebra
1046 # Default these to what they should be if they turn out to be
1047 # trivial, because eigenspaces_left() won't return eigenvalues
1048 # corresponding to trivial spaces (e.g. it returns only the
1049 # eigenspace corresponding to lambda=1 if you take the
1050 # decomposition relative to the identity element).
1051 trivial
= FiniteDimensionalEJASubalgebra(self
, ())
1052 J0
= trivial
# eigenvalue zero
1053 J5
= VectorSpace(self
.base_ring(), 0) # eigenvalue one-half
1054 J1
= trivial
# eigenvalue one
1056 for (eigval
, eigspace
) in c
.operator().matrix().right_eigenspaces():
1057 if eigval
== ~
(self
.base_ring()(2)):
1060 gens
= tuple( self
.from_vector(b
) for b
in eigspace
.basis() )
1061 subalg
= FiniteDimensionalEJASubalgebra(self
,
1069 raise ValueError("unexpected eigenvalue: %s" % eigval
)
1074 def random_element(self
, thorough
=False):
1076 Return a random element of this algebra.
1078 Our algebra superclass method only returns a linear
1079 combination of at most two basis elements. We instead
1080 want the vector space "random element" method that
1081 returns a more diverse selection.
1085 - ``thorough`` -- (boolean; default False) whether or not we
1086 should generate irrational coefficients for the random
1087 element when our base ring is irrational; this slows the
1088 algebra operations to a crawl, but any truly random method
1092 # For a general base ring... maybe we can trust this to do the
1093 # right thing? Unlikely, but.
1094 V
= self
.vector_space()
1095 v
= V
.random_element()
1097 if self
.base_ring() is AA
:
1098 # The "random element" method of the algebraic reals is
1099 # stupid at the moment, and only returns integers between
1100 # -2 and 2, inclusive:
1102 # https://trac.sagemath.org/ticket/30875
1104 # Instead, we implement our own "random vector" method,
1105 # and then coerce that into the algebra. We use the vector
1106 # space degree here instead of the dimension because a
1107 # subalgebra could (for example) be spanned by only two
1108 # vectors, each with five coordinates. We need to
1109 # generate all five coordinates.
1111 v
*= QQbar
.random_element().real()
1113 v
*= QQ
.random_element()
1115 return self
.from_vector(V
.coordinate_vector(v
))
1117 def random_elements(self
, count
, thorough
=False):
1119 Return ``count`` random elements as a tuple.
1123 - ``thorough`` -- (boolean; default False) whether or not we
1124 should generate irrational coefficients for the random
1125 elements when our base ring is irrational; this slows the
1126 algebra operations to a crawl, but any truly random method
1131 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1135 sage: J = JordanSpinEJA(3)
1136 sage: x,y,z = J.random_elements(3)
1137 sage: all( [ x in J, y in J, z in J ])
1139 sage: len( J.random_elements(10) ) == 10
1143 return tuple( self
.random_element(thorough
)
1144 for idx
in range(count
) )
1148 def _charpoly_coefficients(self
):
1150 The `r` polynomial coefficients of the "characteristic polynomial
1155 sage: from mjo.eja.eja_algebra import random_eja
1159 The theory shows that these are all homogeneous polynomials of
1162 sage: set_random_seed()
1163 sage: J = random_eja()
1164 sage: all(p.is_homogeneous() for p in J._charpoly_coefficients())
1168 n
= self
.dimension()
1169 R
= self
.coordinate_polynomial_ring()
1171 F
= R
.fraction_field()
1174 # From a result in my book, these are the entries of the
1175 # basis representation of L_x.
1176 return sum( vars[k
]*self
.gens()[k
].operator().matrix()[i
,j
]
1179 L_x
= matrix(F
, n
, n
, L_x_i_j
)
1182 if self
.rank
.is_in_cache():
1184 # There's no need to pad the system with redundant
1185 # columns if we *know* they'll be redundant.
1188 # Compute an extra power in case the rank is equal to
1189 # the dimension (otherwise, we would stop at x^(r-1)).
1190 x_powers
= [ (L_x
**k
)*self
.one().to_vector()
1191 for k
in range(n
+1) ]
1192 A
= matrix
.column(F
, x_powers
[:n
])
1193 AE
= A
.extended_echelon_form()
1200 # The theory says that only the first "r" coefficients are
1201 # nonzero, and they actually live in the original polynomial
1202 # ring and not the fraction field. We negate them because in
1203 # the actual characteristic polynomial, they get moved to the
1204 # other side where x^r lives. We don't bother to trim A_rref
1205 # down to a square matrix and solve the resulting system,
1206 # because the upper-left r-by-r portion of A_rref is
1207 # guaranteed to be the identity matrix, so e.g.
1209 # A_rref.solve_right(Y)
1211 # would just be returning Y.
1212 return (-E
*b
)[:r
].change_ring(R
)
1217 Return the rank of this EJA.
1219 This is a cached method because we know the rank a priori for
1220 all of the algebras we can construct. Thus we can avoid the
1221 expensive ``_charpoly_coefficients()`` call unless we truly
1222 need to compute the whole characteristic polynomial.
1226 sage: from mjo.eja.eja_algebra import (HadamardEJA,
1227 ....: JordanSpinEJA,
1228 ....: RealSymmetricEJA,
1229 ....: ComplexHermitianEJA,
1230 ....: QuaternionHermitianEJA,
1235 The rank of the Jordan spin algebra is always two::
1237 sage: JordanSpinEJA(2).rank()
1239 sage: JordanSpinEJA(3).rank()
1241 sage: JordanSpinEJA(4).rank()
1244 The rank of the `n`-by-`n` Hermitian real, complex, or
1245 quaternion matrices is `n`::
1247 sage: RealSymmetricEJA(4).rank()
1249 sage: ComplexHermitianEJA(3).rank()
1251 sage: QuaternionHermitianEJA(2).rank()
1256 Ensure that every EJA that we know how to construct has a
1257 positive integer rank, unless the algebra is trivial in
1258 which case its rank will be zero::
1260 sage: set_random_seed()
1261 sage: J = random_eja()
1265 sage: r > 0 or (r == 0 and J.is_trivial())
1268 Ensure that computing the rank actually works, since the ranks
1269 of all simple algebras are known and will be cached by default::
1271 sage: set_random_seed() # long time
1272 sage: J = random_eja() # long time
1273 sage: cached = J.rank() # long time
1274 sage: J.rank.clear_cache() # long time
1275 sage: J.rank() == cached # long time
1279 return len(self
._charpoly
_coefficients
())
1282 def vector_space(self
):
1284 Return the vector space that underlies this algebra.
1288 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1292 sage: J = RealSymmetricEJA(2)
1293 sage: J.vector_space()
1294 Vector space of dimension 3 over...
1297 return self
.zero().to_vector().parent().ambient_vector_space()
1300 Element
= FiniteDimensionalEJAElement
1302 class RationalBasisEJA(FiniteDimensionalEJA
):
1304 New class for algebras whose supplied basis elements have all rational entries.
1308 sage: from mjo.eja.eja_algebra import BilinearFormEJA
1312 The supplied basis is orthonormalized by default::
1314 sage: B = matrix(QQ, [[1, 0, 0], [0, 25, -32], [0, -32, 41]])
1315 sage: J = BilinearFormEJA(B)
1316 sage: J.matrix_basis()
1333 # Abuse the check_field parameter to check that the entries of
1334 # out basis (in ambient coordinates) are in the field QQ.
1335 if not all( all(b_i
in QQ
for b_i
in b
.list()) for b
in basis
):
1336 raise TypeError("basis not rational")
1338 self
._rational
_algebra
= None
1340 # There's no point in constructing the extra algebra if this
1341 # one is already rational.
1343 # Note: the same Jordan and inner-products work here,
1344 # because they are necessarily defined with respect to
1345 # ambient coordinates and not any particular basis.
1346 self
._rational
_algebra
= FiniteDimensionalEJA(
1351 orthonormalize
=False,
1355 super().__init
__(basis
,
1359 check_field
=check_field
,
1363 def _charpoly_coefficients(self
):
1367 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
1368 ....: JordanSpinEJA)
1372 The base ring of the resulting polynomial coefficients is what
1373 it should be, and not the rationals (unless the algebra was
1374 already over the rationals)::
1376 sage: J = JordanSpinEJA(3)
1377 sage: J._charpoly_coefficients()
1378 (X1^2 - X2^2 - X3^2, -2*X1)
1379 sage: a0 = J._charpoly_coefficients()[0]
1381 Algebraic Real Field
1382 sage: a0.base_ring()
1383 Algebraic Real Field
1386 if self
._rational
_algebra
is None:
1387 # There's no need to construct *another* algebra over the
1388 # rationals if this one is already over the
1389 # rationals. Likewise, if we never orthonormalized our
1390 # basis, we might as well just use the given one.
1391 return super()._charpoly
_coefficients
()
1393 # Do the computation over the rationals. The answer will be
1394 # the same, because all we've done is a change of basis.
1395 # Then, change back from QQ to our real base ring
1396 a
= ( a_i
.change_ring(self
.base_ring())
1397 for a_i
in self
._rational
_algebra
._charpoly
_coefficients
() )
1399 if self
._deortho
_matrix
is None:
1400 # This can happen if our base ring was, say, AA and we
1401 # chose not to (or didn't need to) orthonormalize. It's
1402 # still faster to do the computations over QQ even if
1403 # the numbers in the boxes stay the same.
1406 # Otherwise, convert the coordinate variables back to the
1407 # deorthonormalized ones.
1408 R
= self
.coordinate_polynomial_ring()
1409 from sage
.modules
.free_module_element
import vector
1410 X
= vector(R
, R
.gens())
1411 BX
= self
._deortho
_matrix
*X
1413 subs_dict
= { X[i]: BX[i] for i in range(len(X)) }
1414 return tuple( a_i
.subs(subs_dict
) for a_i
in a
)
1416 class ConcreteEJA(RationalBasisEJA
):
1418 A class for the Euclidean Jordan algebras that we know by name.
1420 These are the Jordan algebras whose basis, multiplication table,
1421 rank, and so on are known a priori. More to the point, they are
1422 the Euclidean Jordan algebras for which we are able to conjure up
1423 a "random instance."
1427 sage: from mjo.eja.eja_algebra import ConcreteEJA
1431 Our basis is normalized with respect to the algebra's inner
1432 product, unless we specify otherwise::
1434 sage: set_random_seed()
1435 sage: J = ConcreteEJA.random_instance()
1436 sage: all( b.norm() == 1 for b in J.gens() )
1439 Since our basis is orthonormal with respect to the algebra's inner
1440 product, and since we know that this algebra is an EJA, any
1441 left-multiplication operator's matrix will be symmetric because
1442 natural->EJA basis representation is an isometry and within the
1443 EJA the operator is self-adjoint by the Jordan axiom::
1445 sage: set_random_seed()
1446 sage: J = ConcreteEJA.random_instance()
1447 sage: x = J.random_element()
1448 sage: x.operator().is_self_adjoint()
1453 def _max_random_instance_size():
1455 Return an integer "size" that is an upper bound on the size of
1456 this algebra when it is used in a random test
1457 case. Unfortunately, the term "size" is ambiguous -- when
1458 dealing with `R^n` under either the Hadamard or Jordan spin
1459 product, the "size" refers to the dimension `n`. When dealing
1460 with a matrix algebra (real symmetric or complex/quaternion
1461 Hermitian), it refers to the size of the matrix, which is far
1462 less than the dimension of the underlying vector space.
1464 This method must be implemented in each subclass.
1466 raise NotImplementedError
1469 def random_instance(cls
, *args
, **kwargs
):
1471 Return a random instance of this type of algebra.
1473 This method should be implemented in each subclass.
1475 from sage
.misc
.prandom
import choice
1476 eja_class
= choice(cls
.__subclasses
__())
1478 # These all bubble up to the RationalBasisEJA superclass
1479 # constructor, so any (kw)args valid there are also valid
1481 return eja_class
.random_instance(*args
, **kwargs
)
1486 def dimension_over_reals():
1488 The dimension of this matrix's base ring over the reals.
1490 The reals are dimension one over themselves, obviously; that's
1491 just `\mathbb{R}^{1}`. Likewise, the complex numbers `a + bi`
1492 have dimension two. Finally, the quaternions have dimension
1493 four over the reals.
1495 This is used to determine the size of the matrix returned from
1496 :meth:`real_embed`, among other things.
1498 raise NotImplementedError
1501 def real_embed(cls
,M
):
1503 Embed the matrix ``M`` into a space of real matrices.
1505 The matrix ``M`` can have entries in any field at the moment:
1506 the real numbers, complex numbers, or quaternions. And although
1507 they are not a field, we can probably support octonions at some
1508 point, too. This function returns a real matrix that "acts like"
1509 the original with respect to matrix multiplication; i.e.
1511 real_embed(M*N) = real_embed(M)*real_embed(N)
1514 if M
.ncols() != M
.nrows():
1515 raise ValueError("the matrix 'M' must be square")
1520 def real_unembed(cls
,M
):
1522 The inverse of :meth:`real_embed`.
1524 if M
.ncols() != M
.nrows():
1525 raise ValueError("the matrix 'M' must be square")
1526 if not ZZ(M
.nrows()).mod(cls
.dimension_over_reals()).is_zero():
1527 raise ValueError("the matrix 'M' must be a real embedding")
1531 def jordan_product(X
,Y
):
1532 return (X
*Y
+ Y
*X
)/2
1535 def trace_inner_product(cls
,X
,Y
):
1537 Compute the trace inner-product of two real-embeddings.
1541 sage: from mjo.eja.eja_algebra import (RealSymmetricEJA,
1542 ....: ComplexHermitianEJA,
1543 ....: QuaternionHermitianEJA)
1547 This gives the same answer as it would if we computed the trace
1548 from the unembedded (original) matrices::
1550 sage: set_random_seed()
1551 sage: J = RealSymmetricEJA.random_instance()
1552 sage: x,y = J.random_elements(2)
1553 sage: Xe = x.to_matrix()
1554 sage: Ye = y.to_matrix()
1555 sage: X = J.real_unembed(Xe)
1556 sage: Y = J.real_unembed(Ye)
1557 sage: expected = (X*Y).trace()
1558 sage: actual = J.trace_inner_product(Xe,Ye)
1559 sage: actual == expected
1564 sage: set_random_seed()
1565 sage: J = ComplexHermitianEJA.random_instance()
1566 sage: x,y = J.random_elements(2)
1567 sage: Xe = x.to_matrix()
1568 sage: Ye = y.to_matrix()
1569 sage: X = J.real_unembed(Xe)
1570 sage: Y = J.real_unembed(Ye)
1571 sage: expected = (X*Y).trace().real()
1572 sage: actual = J.trace_inner_product(Xe,Ye)
1573 sage: actual == expected
1578 sage: set_random_seed()
1579 sage: J = QuaternionHermitianEJA.random_instance()
1580 sage: x,y = J.random_elements(2)
1581 sage: Xe = x.to_matrix()
1582 sage: Ye = y.to_matrix()
1583 sage: X = J.real_unembed(Xe)
1584 sage: Y = J.real_unembed(Ye)
1585 sage: expected = (X*Y).trace().coefficient_tuple()[0]
1586 sage: actual = J.trace_inner_product(Xe,Ye)
1587 sage: actual == expected
1591 Xu
= cls
.real_unembed(X
)
1592 Yu
= cls
.real_unembed(Y
)
1593 tr
= (Xu
*Yu
).trace()
1596 # Works in QQ, AA, RDF, et cetera.
1598 except AttributeError:
1599 # A quaternion doesn't have a real() method, but does
1600 # have coefficient_tuple() method that returns the
1601 # coefficients of 1, i, j, and k -- in that order.
1602 return tr
.coefficient_tuple()[0]
1605 class RealMatrixEJA(MatrixEJA
):
1607 def dimension_over_reals():
1611 class RealSymmetricEJA(ConcreteEJA
, RealMatrixEJA
):
1613 The rank-n simple EJA consisting of real symmetric n-by-n
1614 matrices, the usual symmetric Jordan product, and the trace inner
1615 product. It has dimension `(n^2 + n)/2` over the reals.
1619 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1623 sage: J = RealSymmetricEJA(2)
1624 sage: e0, e1, e2 = J.gens()
1632 In theory, our "field" can be any subfield of the reals::
1634 sage: RealSymmetricEJA(2, field=RDF)
1635 Euclidean Jordan algebra of dimension 3 over Real Double Field
1636 sage: RealSymmetricEJA(2, field=RR)
1637 Euclidean Jordan algebra of dimension 3 over Real Field with
1638 53 bits of precision
1642 The dimension of this algebra is `(n^2 + n) / 2`::
1644 sage: set_random_seed()
1645 sage: n_max = RealSymmetricEJA._max_random_instance_size()
1646 sage: n = ZZ.random_element(1, n_max)
1647 sage: J = RealSymmetricEJA(n)
1648 sage: J.dimension() == (n^2 + n)/2
1651 The Jordan multiplication is what we think it is::
1653 sage: set_random_seed()
1654 sage: J = RealSymmetricEJA.random_instance()
1655 sage: x,y = J.random_elements(2)
1656 sage: actual = (x*y).to_matrix()
1657 sage: X = x.to_matrix()
1658 sage: Y = y.to_matrix()
1659 sage: expected = (X*Y + Y*X)/2
1660 sage: actual == expected
1662 sage: J(expected) == x*y
1665 We can change the generator prefix::
1667 sage: RealSymmetricEJA(3, prefix='q').gens()
1668 (q0, q1, q2, q3, q4, q5)
1670 We can construct the (trivial) algebra of rank zero::
1672 sage: RealSymmetricEJA(0)
1673 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1677 def _denormalized_basis(cls
, n
):
1679 Return a basis for the space of real symmetric n-by-n matrices.
1683 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1687 sage: set_random_seed()
1688 sage: n = ZZ.random_element(1,5)
1689 sage: B = RealSymmetricEJA._denormalized_basis(n)
1690 sage: all( M.is_symmetric() for M in B)
1694 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1698 for j
in range(i
+1):
1699 Eij
= matrix(ZZ
, n
, lambda k
,l
: k
==i
and l
==j
)
1703 Sij
= Eij
+ Eij
.transpose()
1709 def _max_random_instance_size():
1710 return 4 # Dimension 10
1713 def random_instance(cls
, **kwargs
):
1715 Return a random instance of this type of algebra.
1717 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1718 return cls(n
, **kwargs
)
1720 def __init__(self
, n
, **kwargs
):
1721 # We know this is a valid EJA, but will double-check
1722 # if the user passes check_axioms=True.
1723 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
1725 super(RealSymmetricEJA
, self
).__init
__(self
._denormalized
_basis
(n
),
1726 self
.jordan_product
,
1727 self
.trace_inner_product
,
1730 # TODO: this could be factored out somehow, but is left here
1731 # because the MatrixEJA is not presently a subclass of the
1732 # FDEJA class that defines rank() and one().
1733 self
.rank
.set_cache(n
)
1734 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
1735 self
.one
.set_cache(self(idV
))
1739 class ComplexMatrixEJA(MatrixEJA
):
1740 # A manual dictionary-cache for the complex_extension() method,
1741 # since apparently @classmethods can't also be @cached_methods.
1742 _complex_extension
= {}
1745 def complex_extension(cls
,field
):
1747 The complex field that we embed/unembed, as an extension
1748 of the given ``field``.
1750 if field
in cls
._complex
_extension
:
1751 return cls
._complex
_extension
[field
]
1753 # Sage doesn't know how to adjoin the complex "i" (the root of
1754 # x^2 + 1) to a field in a general way. Here, we just enumerate
1755 # all of the cases that I have cared to support so far.
1757 # Sage doesn't know how to embed AA into QQbar, i.e. how
1758 # to adjoin sqrt(-1) to AA.
1760 elif not field
.is_exact():
1762 F
= field
.complex_field()
1764 # Works for QQ and... maybe some other fields.
1765 R
= PolynomialRing(field
, 'z')
1767 F
= field
.extension(z
**2 + 1, 'I', embedding
=CLF(-1).sqrt())
1769 cls
._complex
_extension
[field
] = F
1773 def dimension_over_reals():
1777 def real_embed(cls
,M
):
1779 Embed the n-by-n complex matrix ``M`` into the space of real
1780 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1781 bi` to the block matrix ``[[a,b],[-b,a]]``.
1785 sage: from mjo.eja.eja_algebra import ComplexMatrixEJA
1789 sage: F = QuadraticField(-1, 'I')
1790 sage: x1 = F(4 - 2*i)
1791 sage: x2 = F(1 + 2*i)
1794 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1795 sage: ComplexMatrixEJA.real_embed(M)
1804 Embedding is a homomorphism (isomorphism, in fact)::
1806 sage: set_random_seed()
1807 sage: n = ZZ.random_element(3)
1808 sage: F = QuadraticField(-1, 'I')
1809 sage: X = random_matrix(F, n)
1810 sage: Y = random_matrix(F, n)
1811 sage: Xe = ComplexMatrixEJA.real_embed(X)
1812 sage: Ye = ComplexMatrixEJA.real_embed(Y)
1813 sage: XYe = ComplexMatrixEJA.real_embed(X*Y)
1818 super(ComplexMatrixEJA
,cls
).real_embed(M
)
1821 # We don't need any adjoined elements...
1822 field
= M
.base_ring().base_ring()
1828 blocks
.append(matrix(field
, 2, [ [ a
, b
],
1831 return matrix
.block(field
, n
, blocks
)
1835 def real_unembed(cls
,M
):
1837 The inverse of _embed_complex_matrix().
1841 sage: from mjo.eja.eja_algebra import ComplexMatrixEJA
1845 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1846 ....: [-2, 1, -4, 3],
1847 ....: [ 9, 10, 11, 12],
1848 ....: [-10, 9, -12, 11] ])
1849 sage: ComplexMatrixEJA.real_unembed(A)
1851 [ 10*I + 9 12*I + 11]
1855 Unembedding is the inverse of embedding::
1857 sage: set_random_seed()
1858 sage: F = QuadraticField(-1, 'I')
1859 sage: M = random_matrix(F, 3)
1860 sage: Me = ComplexMatrixEJA.real_embed(M)
1861 sage: ComplexMatrixEJA.real_unembed(Me) == M
1865 super(ComplexMatrixEJA
,cls
).real_unembed(M
)
1867 d
= cls
.dimension_over_reals()
1868 F
= cls
.complex_extension(M
.base_ring())
1871 # Go top-left to bottom-right (reading order), converting every
1872 # 2-by-2 block we see to a single complex element.
1874 for k
in range(n
/d
):
1875 for j
in range(n
/d
):
1876 submat
= M
[d
*k
:d
*k
+d
,d
*j
:d
*j
+d
]
1877 if submat
[0,0] != submat
[1,1]:
1878 raise ValueError('bad on-diagonal submatrix')
1879 if submat
[0,1] != -submat
[1,0]:
1880 raise ValueError('bad off-diagonal submatrix')
1881 z
= submat
[0,0] + submat
[0,1]*i
1884 return matrix(F
, n
/d
, elements
)
1887 class ComplexHermitianEJA(ConcreteEJA
, ComplexMatrixEJA
):
1889 The rank-n simple EJA consisting of complex Hermitian n-by-n
1890 matrices over the real numbers, the usual symmetric Jordan product,
1891 and the real-part-of-trace inner product. It has dimension `n^2` over
1896 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1900 In theory, our "field" can be any subfield of the reals::
1902 sage: ComplexHermitianEJA(2, field=RDF)
1903 Euclidean Jordan algebra of dimension 4 over Real Double Field
1904 sage: ComplexHermitianEJA(2, field=RR)
1905 Euclidean Jordan algebra of dimension 4 over Real Field with
1906 53 bits of precision
1910 The dimension of this algebra is `n^2`::
1912 sage: set_random_seed()
1913 sage: n_max = ComplexHermitianEJA._max_random_instance_size()
1914 sage: n = ZZ.random_element(1, n_max)
1915 sage: J = ComplexHermitianEJA(n)
1916 sage: J.dimension() == n^2
1919 The Jordan multiplication is what we think it is::
1921 sage: set_random_seed()
1922 sage: J = ComplexHermitianEJA.random_instance()
1923 sage: x,y = J.random_elements(2)
1924 sage: actual = (x*y).to_matrix()
1925 sage: X = x.to_matrix()
1926 sage: Y = y.to_matrix()
1927 sage: expected = (X*Y + Y*X)/2
1928 sage: actual == expected
1930 sage: J(expected) == x*y
1933 We can change the generator prefix::
1935 sage: ComplexHermitianEJA(2, prefix='z').gens()
1938 We can construct the (trivial) algebra of rank zero::
1940 sage: ComplexHermitianEJA(0)
1941 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1946 def _denormalized_basis(cls
, n
):
1948 Returns a basis for the space of complex Hermitian n-by-n matrices.
1950 Why do we embed these? Basically, because all of numerical linear
1951 algebra assumes that you're working with vectors consisting of `n`
1952 entries from a field and scalars from the same field. There's no way
1953 to tell SageMath that (for example) the vectors contain complex
1954 numbers, while the scalar field is real.
1958 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1962 sage: set_random_seed()
1963 sage: n = ZZ.random_element(1,5)
1964 sage: B = ComplexHermitianEJA._denormalized_basis(n)
1965 sage: all( M.is_symmetric() for M in B)
1970 R
= PolynomialRing(field
, 'z')
1972 F
= field
.extension(z
**2 + 1, 'I')
1975 # This is like the symmetric case, but we need to be careful:
1977 # * We want conjugate-symmetry, not just symmetry.
1978 # * The diagonal will (as a result) be real.
1981 Eij
= matrix
.zero(F
,n
)
1983 for j
in range(i
+1):
1987 Sij
= cls
.real_embed(Eij
)
1990 # The second one has a minus because it's conjugated.
1991 Eij
[j
,i
] = 1 # Eij = Eij + Eij.transpose()
1992 Sij_real
= cls
.real_embed(Eij
)
1994 # Eij = I*Eij - I*Eij.transpose()
1997 Sij_imag
= cls
.real_embed(Eij
)
2003 # Since we embedded these, we can drop back to the "field" that we
2004 # started with instead of the complex extension "F".
2005 return tuple( s
.change_ring(field
) for s
in S
)
2008 def __init__(self
, n
, **kwargs
):
2009 # We know this is a valid EJA, but will double-check
2010 # if the user passes check_axioms=True.
2011 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2013 super(ComplexHermitianEJA
, self
).__init
__(self
._denormalized
_basis
(n
),
2014 self
.jordan_product
,
2015 self
.trace_inner_product
,
2017 # TODO: this could be factored out somehow, but is left here
2018 # because the MatrixEJA is not presently a subclass of the
2019 # FDEJA class that defines rank() and one().
2020 self
.rank
.set_cache(n
)
2021 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
2022 self
.one
.set_cache(self(idV
))
2025 def _max_random_instance_size():
2026 return 3 # Dimension 9
2029 def random_instance(cls
, **kwargs
):
2031 Return a random instance of this type of algebra.
2033 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2034 return cls(n
, **kwargs
)
2036 class QuaternionMatrixEJA(MatrixEJA
):
2038 # A manual dictionary-cache for the quaternion_extension() method,
2039 # since apparently @classmethods can't also be @cached_methods.
2040 _quaternion_extension
= {}
2043 def quaternion_extension(cls
,field
):
2045 The quaternion field that we embed/unembed, as an extension
2046 of the given ``field``.
2048 if field
in cls
._quaternion
_extension
:
2049 return cls
._quaternion
_extension
[field
]
2051 Q
= QuaternionAlgebra(field
,-1,-1)
2053 cls
._quaternion
_extension
[field
] = Q
2057 def dimension_over_reals():
2061 def real_embed(cls
,M
):
2063 Embed the n-by-n quaternion matrix ``M`` into the space of real
2064 matrices of size 4n-by-4n by first sending each quaternion entry `z
2065 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
2066 c+di],[-c + di, a-bi]]`, and then embedding those into a real
2071 sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA
2075 sage: Q = QuaternionAlgebra(QQ,-1,-1)
2076 sage: i,j,k = Q.gens()
2077 sage: x = 1 + 2*i + 3*j + 4*k
2078 sage: M = matrix(Q, 1, [[x]])
2079 sage: QuaternionMatrixEJA.real_embed(M)
2085 Embedding is a homomorphism (isomorphism, in fact)::
2087 sage: set_random_seed()
2088 sage: n = ZZ.random_element(2)
2089 sage: Q = QuaternionAlgebra(QQ,-1,-1)
2090 sage: X = random_matrix(Q, n)
2091 sage: Y = random_matrix(Q, n)
2092 sage: Xe = QuaternionMatrixEJA.real_embed(X)
2093 sage: Ye = QuaternionMatrixEJA.real_embed(Y)
2094 sage: XYe = QuaternionMatrixEJA.real_embed(X*Y)
2099 super(QuaternionMatrixEJA
,cls
).real_embed(M
)
2100 quaternions
= M
.base_ring()
2103 F
= QuadraticField(-1, 'I')
2108 t
= z
.coefficient_tuple()
2113 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
2114 [-c
+ d
*i
, a
- b
*i
]])
2115 realM
= ComplexMatrixEJA
.real_embed(cplxM
)
2116 blocks
.append(realM
)
2118 # We should have real entries by now, so use the realest field
2119 # we've got for the return value.
2120 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
2125 def real_unembed(cls
,M
):
2127 The inverse of _embed_quaternion_matrix().
2131 sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA
2135 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
2136 ....: [-2, 1, -4, 3],
2137 ....: [-3, 4, 1, -2],
2138 ....: [-4, -3, 2, 1]])
2139 sage: QuaternionMatrixEJA.real_unembed(M)
2140 [1 + 2*i + 3*j + 4*k]
2144 Unembedding is the inverse of embedding::
2146 sage: set_random_seed()
2147 sage: Q = QuaternionAlgebra(QQ, -1, -1)
2148 sage: M = random_matrix(Q, 3)
2149 sage: Me = QuaternionMatrixEJA.real_embed(M)
2150 sage: QuaternionMatrixEJA.real_unembed(Me) == M
2154 super(QuaternionMatrixEJA
,cls
).real_unembed(M
)
2156 d
= cls
.dimension_over_reals()
2158 # Use the base ring of the matrix to ensure that its entries can be
2159 # multiplied by elements of the quaternion algebra.
2160 Q
= cls
.quaternion_extension(M
.base_ring())
2163 # Go top-left to bottom-right (reading order), converting every
2164 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
2167 for l
in range(n
/d
):
2168 for m
in range(n
/d
):
2169 submat
= ComplexMatrixEJA
.real_unembed(
2170 M
[d
*l
:d
*l
+d
,d
*m
:d
*m
+d
] )
2171 if submat
[0,0] != submat
[1,1].conjugate():
2172 raise ValueError('bad on-diagonal submatrix')
2173 if submat
[0,1] != -submat
[1,0].conjugate():
2174 raise ValueError('bad off-diagonal submatrix')
2175 z
= submat
[0,0].real()
2176 z
+= submat
[0,0].imag()*i
2177 z
+= submat
[0,1].real()*j
2178 z
+= submat
[0,1].imag()*k
2181 return matrix(Q
, n
/d
, elements
)
2184 class QuaternionHermitianEJA(ConcreteEJA
, QuaternionMatrixEJA
):
2186 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
2187 matrices, the usual symmetric Jordan product, and the
2188 real-part-of-trace inner product. It has dimension `2n^2 - n` over
2193 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2197 In theory, our "field" can be any subfield of the reals::
2199 sage: QuaternionHermitianEJA(2, field=RDF)
2200 Euclidean Jordan algebra of dimension 6 over Real Double Field
2201 sage: QuaternionHermitianEJA(2, field=RR)
2202 Euclidean Jordan algebra of dimension 6 over Real Field with
2203 53 bits of precision
2207 The dimension of this algebra is `2*n^2 - n`::
2209 sage: set_random_seed()
2210 sage: n_max = QuaternionHermitianEJA._max_random_instance_size()
2211 sage: n = ZZ.random_element(1, n_max)
2212 sage: J = QuaternionHermitianEJA(n)
2213 sage: J.dimension() == 2*(n^2) - n
2216 The Jordan multiplication is what we think it is::
2218 sage: set_random_seed()
2219 sage: J = QuaternionHermitianEJA.random_instance()
2220 sage: x,y = J.random_elements(2)
2221 sage: actual = (x*y).to_matrix()
2222 sage: X = x.to_matrix()
2223 sage: Y = y.to_matrix()
2224 sage: expected = (X*Y + Y*X)/2
2225 sage: actual == expected
2227 sage: J(expected) == x*y
2230 We can change the generator prefix::
2232 sage: QuaternionHermitianEJA(2, prefix='a').gens()
2233 (a0, a1, a2, a3, a4, a5)
2235 We can construct the (trivial) algebra of rank zero::
2237 sage: QuaternionHermitianEJA(0)
2238 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2242 def _denormalized_basis(cls
, n
):
2244 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
2246 Why do we embed these? Basically, because all of numerical
2247 linear algebra assumes that you're working with vectors consisting
2248 of `n` entries from a field and scalars from the same field. There's
2249 no way to tell SageMath that (for example) the vectors contain
2250 complex numbers, while the scalar field is real.
2254 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2258 sage: set_random_seed()
2259 sage: n = ZZ.random_element(1,5)
2260 sage: B = QuaternionHermitianEJA._denormalized_basis(n)
2261 sage: all( M.is_symmetric() for M in B )
2266 Q
= QuaternionAlgebra(QQ
,-1,-1)
2269 # This is like the symmetric case, but we need to be careful:
2271 # * We want conjugate-symmetry, not just symmetry.
2272 # * The diagonal will (as a result) be real.
2275 Eij
= matrix
.zero(Q
,n
)
2277 for j
in range(i
+1):
2281 Sij
= cls
.real_embed(Eij
)
2284 # The second, third, and fourth ones have a minus
2285 # because they're conjugated.
2286 # Eij = Eij + Eij.transpose()
2288 Sij_real
= cls
.real_embed(Eij
)
2290 # Eij = I*(Eij - Eij.transpose())
2293 Sij_I
= cls
.real_embed(Eij
)
2295 # Eij = J*(Eij - Eij.transpose())
2298 Sij_J
= cls
.real_embed(Eij
)
2300 # Eij = K*(Eij - Eij.transpose())
2303 Sij_K
= cls
.real_embed(Eij
)
2309 # Since we embedded these, we can drop back to the "field" that we
2310 # started with instead of the quaternion algebra "Q".
2311 return tuple( s
.change_ring(field
) for s
in S
)
2314 def __init__(self
, n
, **kwargs
):
2315 # We know this is a valid EJA, but will double-check
2316 # if the user passes check_axioms=True.
2317 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2319 super(QuaternionHermitianEJA
, self
).__init
__(self
._denormalized
_basis
(n
),
2320 self
.jordan_product
,
2321 self
.trace_inner_product
,
2323 # TODO: this could be factored out somehow, but is left here
2324 # because the MatrixEJA is not presently a subclass of the
2325 # FDEJA class that defines rank() and one().
2326 self
.rank
.set_cache(n
)
2327 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
2328 self
.one
.set_cache(self(idV
))
2332 def _max_random_instance_size():
2334 The maximum rank of a random QuaternionHermitianEJA.
2336 return 2 # Dimension 6
2339 def random_instance(cls
, **kwargs
):
2341 Return a random instance of this type of algebra.
2343 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2344 return cls(n
, **kwargs
)
2347 class HadamardEJA(ConcreteEJA
):
2349 Return the Euclidean Jordan Algebra corresponding to the set
2350 `R^n` under the Hadamard product.
2352 Note: this is nothing more than the Cartesian product of ``n``
2353 copies of the spin algebra. Once Cartesian product algebras
2354 are implemented, this can go.
2358 sage: from mjo.eja.eja_algebra import HadamardEJA
2362 This multiplication table can be verified by hand::
2364 sage: J = HadamardEJA(3)
2365 sage: e0,e1,e2 = J.gens()
2381 We can change the generator prefix::
2383 sage: HadamardEJA(3, prefix='r').gens()
2387 def __init__(self
, n
, **kwargs
):
2389 jordan_product
= lambda x
,y
: x
2390 inner_product
= lambda x
,y
: x
2392 def jordan_product(x
,y
):
2394 return P( xi
*yi
for (xi
,yi
) in zip(x
,y
) )
2396 def inner_product(x
,y
):
2399 # New defaults for keyword arguments. Don't orthonormalize
2400 # because our basis is already orthonormal with respect to our
2401 # inner-product. Don't check the axioms, because we know this
2402 # is a valid EJA... but do double-check if the user passes
2403 # check_axioms=True. Note: we DON'T override the "check_field"
2404 # default here, because the user can pass in a field!
2405 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2406 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2408 column_basis
= tuple( b
.column() for b
in FreeModule(ZZ
, n
).basis() )
2409 super().__init
__(column_basis
,
2414 self
.rank
.set_cache(n
)
2417 self
.one
.set_cache( self
.zero() )
2419 self
.one
.set_cache( sum(self
.gens()) )
2422 def _max_random_instance_size():
2424 The maximum dimension of a random HadamardEJA.
2429 def random_instance(cls
, **kwargs
):
2431 Return a random instance of this type of algebra.
2433 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2434 return cls(n
, **kwargs
)
2437 class BilinearFormEJA(ConcreteEJA
):
2439 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2440 with the half-trace inner product and jordan product ``x*y =
2441 (<Bx,y>,y_bar>, x0*y_bar + y0*x_bar)`` where `B = 1 \times B22` is
2442 a symmetric positive-definite "bilinear form" matrix. Its
2443 dimension is the size of `B`, and it has rank two in dimensions
2444 larger than two. It reduces to the ``JordanSpinEJA`` when `B` is
2445 the identity matrix of order ``n``.
2447 We insist that the one-by-one upper-left identity block of `B` be
2448 passed in as well so that we can be passed a matrix of size zero
2449 to construct a trivial algebra.
2453 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
2454 ....: JordanSpinEJA)
2458 When no bilinear form is specified, the identity matrix is used,
2459 and the resulting algebra is the Jordan spin algebra::
2461 sage: B = matrix.identity(AA,3)
2462 sage: J0 = BilinearFormEJA(B)
2463 sage: J1 = JordanSpinEJA(3)
2464 sage: J0.multiplication_table() == J0.multiplication_table()
2467 An error is raised if the matrix `B` does not correspond to a
2468 positive-definite bilinear form::
2470 sage: B = matrix.random(QQ,2,3)
2471 sage: J = BilinearFormEJA(B)
2472 Traceback (most recent call last):
2474 ValueError: bilinear form is not positive-definite
2475 sage: B = matrix.zero(QQ,3)
2476 sage: J = BilinearFormEJA(B)
2477 Traceback (most recent call last):
2479 ValueError: bilinear form is not positive-definite
2483 We can create a zero-dimensional algebra::
2485 sage: B = matrix.identity(AA,0)
2486 sage: J = BilinearFormEJA(B)
2490 We can check the multiplication condition given in the Jordan, von
2491 Neumann, and Wigner paper (and also discussed on my "On the
2492 symmetry..." paper). Note that this relies heavily on the standard
2493 choice of basis, as does anything utilizing the bilinear form
2494 matrix. We opt not to orthonormalize the basis, because if we
2495 did, we would have to normalize the `s_{i}` in a similar manner::
2497 sage: set_random_seed()
2498 sage: n = ZZ.random_element(5)
2499 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2500 sage: B11 = matrix.identity(QQ,1)
2501 sage: B22 = M.transpose()*M
2502 sage: B = block_matrix(2,2,[ [B11,0 ],
2504 sage: J = BilinearFormEJA(B, orthonormalize=False)
2505 sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
2506 sage: V = J.vector_space()
2507 sage: sis = [ J( V([0] + (M.inverse()*ei).list()).column() )
2508 ....: for ei in eis ]
2509 sage: actual = [ sis[i]*sis[j]
2510 ....: for i in range(n-1)
2511 ....: for j in range(n-1) ]
2512 sage: expected = [ J.one() if i == j else J.zero()
2513 ....: for i in range(n-1)
2514 ....: for j in range(n-1) ]
2515 sage: actual == expected
2519 def __init__(self
, B
, **kwargs
):
2520 # The matrix "B" is supplied by the user in most cases,
2521 # so it makes sense to check whether or not its positive-
2522 # definite unless we are specifically asked not to...
2523 if ("check_axioms" not in kwargs
) or kwargs
["check_axioms"]:
2524 if not B
.is_positive_definite():
2525 raise ValueError("bilinear form is not positive-definite")
2527 # However, all of the other data for this EJA is computed
2528 # by us in manner that guarantees the axioms are
2529 # satisfied. So, again, unless we are specifically asked to
2530 # verify things, we'll skip the rest of the checks.
2531 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2533 def inner_product(x
,y
):
2534 return (y
.T
*B
*x
)[0,0]
2536 def jordan_product(x
,y
):
2542 z0
= inner_product(y
,x
)
2543 zbar
= y0
*xbar
+ x0
*ybar
2544 return P([z0
] + zbar
.list())
2547 column_basis
= tuple( b
.column() for b
in FreeModule(ZZ
, n
).basis() )
2548 super(BilinearFormEJA
, self
).__init
__(column_basis
,
2553 # The rank of this algebra is two, unless we're in a
2554 # one-dimensional ambient space (because the rank is bounded
2555 # by the ambient dimension).
2556 self
.rank
.set_cache(min(n
,2))
2559 self
.one
.set_cache( self
.zero() )
2561 self
.one
.set_cache( self
.monomial(0) )
2564 def _max_random_instance_size():
2566 The maximum dimension of a random BilinearFormEJA.
2571 def random_instance(cls
, **kwargs
):
2573 Return a random instance of this algebra.
2575 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2577 B
= matrix
.identity(ZZ
, n
)
2578 return cls(B
, **kwargs
)
2580 B11
= matrix
.identity(ZZ
, 1)
2581 M
= matrix
.random(ZZ
, n
-1)
2582 I
= matrix
.identity(ZZ
, n
-1)
2584 while alpha
.is_zero():
2585 alpha
= ZZ
.random_element().abs()
2586 B22
= M
.transpose()*M
+ alpha
*I
2588 from sage
.matrix
.special
import block_matrix
2589 B
= block_matrix(2,2, [ [B11
, ZZ(0) ],
2592 return cls(B
, **kwargs
)
2595 class JordanSpinEJA(BilinearFormEJA
):
2597 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2598 with the usual inner product and jordan product ``x*y =
2599 (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2604 sage: from mjo.eja.eja_algebra import JordanSpinEJA
2608 This multiplication table can be verified by hand::
2610 sage: J = JordanSpinEJA(4)
2611 sage: e0,e1,e2,e3 = J.gens()
2627 We can change the generator prefix::
2629 sage: JordanSpinEJA(2, prefix='B').gens()
2634 Ensure that we have the usual inner product on `R^n`::
2636 sage: set_random_seed()
2637 sage: J = JordanSpinEJA.random_instance()
2638 sage: x,y = J.random_elements(2)
2639 sage: actual = x.inner_product(y)
2640 sage: expected = x.to_vector().inner_product(y.to_vector())
2641 sage: actual == expected
2645 def __init__(self
, n
, **kwargs
):
2646 # This is a special case of the BilinearFormEJA with the
2647 # identity matrix as its bilinear form.
2648 B
= matrix
.identity(ZZ
, n
)
2650 # Don't orthonormalize because our basis is already
2651 # orthonormal with respect to our inner-product.
2652 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2654 # But also don't pass check_field=False here, because the user
2655 # can pass in a field!
2656 super(JordanSpinEJA
, self
).__init
__(B
, **kwargs
)
2659 def _max_random_instance_size():
2661 The maximum dimension of a random JordanSpinEJA.
2666 def random_instance(cls
, **kwargs
):
2668 Return a random instance of this type of algebra.
2670 Needed here to override the implementation for ``BilinearFormEJA``.
2672 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2673 return cls(n
, **kwargs
)
2676 class TrivialEJA(ConcreteEJA
):
2678 The trivial Euclidean Jordan algebra consisting of only a zero element.
2682 sage: from mjo.eja.eja_algebra import TrivialEJA
2686 sage: J = TrivialEJA()
2693 sage: 7*J.one()*12*J.one()
2695 sage: J.one().inner_product(J.one())
2697 sage: J.one().norm()
2699 sage: J.one().subalgebra_generated_by()
2700 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2705 def __init__(self
, **kwargs
):
2706 jordan_product
= lambda x
,y
: x
2707 inner_product
= lambda x
,y
: 0
2710 # New defaults for keyword arguments
2711 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2712 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2714 super(TrivialEJA
, self
).__init
__(basis
,
2718 # The rank is zero using my definition, namely the dimension of the
2719 # largest subalgebra generated by any element.
2720 self
.rank
.set_cache(0)
2721 self
.one
.set_cache( self
.zero() )
2724 def random_instance(cls
, **kwargs
):
2725 # We don't take a "size" argument so the superclass method is
2726 # inappropriate for us.
2727 return cls(**kwargs
)
2730 class CartesianProductEJA(CombinatorialFreeModule_CartesianProduct
,
2731 FiniteDimensionalEJA
):
2733 The external (orthogonal) direct sum of two or more Euclidean
2734 Jordan algebras. Every Euclidean Jordan algebra decomposes into an
2735 orthogonal direct sum of simple Euclidean Jordan algebras which is
2736 then isometric to a Cartesian product, so no generality is lost by
2737 providing only this construction.
2741 sage: from mjo.eja.eja_algebra import (random_eja,
2742 ....: CartesianProductEJA,
2744 ....: JordanSpinEJA,
2745 ....: RealSymmetricEJA)
2749 The Jordan product is inherited from our factors and implemented by
2750 our CombinatorialFreeModule Cartesian product superclass::
2752 sage: set_random_seed()
2753 sage: J1 = HadamardEJA(2)
2754 sage: J2 = RealSymmetricEJA(2)
2755 sage: J = cartesian_product([J1,J2])
2756 sage: x,y = J.random_elements(2)
2760 The ability to retrieve the original factors is implemented by our
2761 CombinatorialFreeModule Cartesian product superclass::
2763 sage: J1 = HadamardEJA(2, field=QQ)
2764 sage: J2 = JordanSpinEJA(3, field=QQ)
2765 sage: J = cartesian_product([J1,J2])
2766 sage: J.cartesian_factors()
2767 (Euclidean Jordan algebra of dimension 2 over Rational Field,
2768 Euclidean Jordan algebra of dimension 3 over Rational Field)
2770 You can provide more than two factors::
2772 sage: J1 = HadamardEJA(2)
2773 sage: J2 = JordanSpinEJA(3)
2774 sage: J3 = RealSymmetricEJA(3)
2775 sage: cartesian_product([J1,J2,J3])
2776 Euclidean Jordan algebra of dimension 2 over Algebraic Real
2777 Field (+) Euclidean Jordan algebra of dimension 3 over Algebraic
2778 Real Field (+) Euclidean Jordan algebra of dimension 6 over
2779 Algebraic Real Field
2781 Rank is additive on a Cartesian product::
2783 sage: J1 = HadamardEJA(1)
2784 sage: J2 = RealSymmetricEJA(2)
2785 sage: J = cartesian_product([J1,J2])
2786 sage: J1.rank.clear_cache()
2787 sage: J2.rank.clear_cache()
2788 sage: J.rank.clear_cache()
2791 sage: J.rank() == J1.rank() + J2.rank()
2794 The same rank computation works over the rationals, with whatever
2797 sage: J1 = HadamardEJA(1, field=QQ, orthonormalize=False)
2798 sage: J2 = RealSymmetricEJA(2, field=QQ, orthonormalize=False)
2799 sage: J = cartesian_product([J1,J2])
2800 sage: J1.rank.clear_cache()
2801 sage: J2.rank.clear_cache()
2802 sage: J.rank.clear_cache()
2805 sage: J.rank() == J1.rank() + J2.rank()
2808 The product algebra will be associative if and only if all of its
2809 components are associative::
2811 sage: J1 = HadamardEJA(2)
2812 sage: J1.is_associative()
2814 sage: J2 = HadamardEJA(3)
2815 sage: J2.is_associative()
2817 sage: J3 = RealSymmetricEJA(3)
2818 sage: J3.is_associative()
2820 sage: CP1 = cartesian_product([J1,J2])
2821 sage: CP1.is_associative()
2823 sage: CP2 = cartesian_product([J1,J3])
2824 sage: CP2.is_associative()
2829 All factors must share the same base field::
2831 sage: J1 = HadamardEJA(2, field=QQ)
2832 sage: J2 = RealSymmetricEJA(2)
2833 sage: CartesianProductEJA((J1,J2))
2834 Traceback (most recent call last):
2836 ValueError: all factors must share the same base field
2838 The "cached" Jordan and inner products are the componentwise
2841 sage: set_random_seed()
2842 sage: J1 = random_eja()
2843 sage: J2 = random_eja()
2844 sage: J = cartesian_product([J1,J2])
2845 sage: x,y = J.random_elements(2)
2846 sage: x*y == J.cartesian_jordan_product(x,y)
2848 sage: x.inner_product(y) == J.cartesian_inner_product(x,y)
2851 The cached unit element is the same one that would be computed::
2853 sage: set_random_seed() # long time
2854 sage: J1 = random_eja() # long time
2855 sage: J2 = random_eja() # long time
2856 sage: J = cartesian_product([J1,J2]) # long time
2857 sage: actual = J.one() # long time
2858 sage: J.one.clear_cache() # long time
2859 sage: expected = J.one() # long time
2860 sage: actual == expected # long time
2864 def __init__(self
, algebras
, **kwargs
):
2865 CombinatorialFreeModule_CartesianProduct
.__init
__(self
,
2868 field
= algebras
[0].base_ring()
2869 if not all( J
.base_ring() == field
for J
in algebras
):
2870 raise ValueError("all factors must share the same base field")
2872 associative
= all( m
.is_associative() for m
in algebras
)
2874 # The definition of matrix_space() and self.basis() relies
2875 # only on the stuff in the CFM_CartesianProduct class, which
2876 # we've already initialized.
2877 Js
= self
.cartesian_factors()
2879 MS
= self
.matrix_space()
2881 MS(tuple( self
.cartesian_projection(i
)(b
).to_matrix()
2882 for i
in range(m
) ))
2883 for b
in self
.basis()
2886 # Define jordan/inner products that operate on that matrix_basis.
2887 def jordan_product(x
,y
):
2889 (Js
[i
](x
[i
])*Js
[i
](y
[i
])).to_matrix() for i
in range(m
)
2892 def inner_product(x
, y
):
2894 Js
[i
](x
[i
]).inner_product(Js
[i
](y
[i
])) for i
in range(m
)
2897 # There's no need to check the field since it already came
2898 # from an EJA. Likewise the axioms are guaranteed to be
2899 # satisfied, unless the guy writing this class sucks.
2901 # If you want the basis to be orthonormalized, orthonormalize
2903 FiniteDimensionalEJA
.__init
__(self
,
2908 orthonormalize
=False,
2909 associative
=associative
,
2910 cartesian_product
=True,
2914 ones
= tuple(J
.one() for J
in algebras
)
2915 self
.one
.set_cache(self
._cartesian
_product
_of
_elements
(ones
))
2916 self
.rank
.set_cache(sum(J
.rank() for J
in algebras
))
2918 def matrix_space(self
):
2920 Return the space that our matrix basis lives in as a Cartesian
2925 sage: from mjo.eja.eja_algebra import (HadamardEJA,
2926 ....: RealSymmetricEJA)
2930 sage: J1 = HadamardEJA(1)
2931 sage: J2 = RealSymmetricEJA(2)
2932 sage: J = cartesian_product([J1,J2])
2933 sage: J.matrix_space()
2934 The Cartesian product of (Full MatrixSpace of 1 by 1 dense
2935 matrices over Algebraic Real Field, Full MatrixSpace of 2
2936 by 2 dense matrices over Algebraic Real Field)
2939 from sage
.categories
.cartesian_product
import cartesian_product
2940 return cartesian_product( [J
.matrix_space()
2941 for J
in self
.cartesian_factors()] )
2944 def cartesian_projection(self
, i
):
2948 sage: from mjo.eja.eja_algebra import (random_eja,
2949 ....: JordanSpinEJA,
2951 ....: RealSymmetricEJA,
2952 ....: ComplexHermitianEJA)
2956 The projection morphisms are Euclidean Jordan algebra
2959 sage: J1 = HadamardEJA(2)
2960 sage: J2 = RealSymmetricEJA(2)
2961 sage: J = cartesian_product([J1,J2])
2962 sage: J.cartesian_projection(0)
2963 Linear operator between finite-dimensional Euclidean Jordan
2964 algebras represented by the matrix:
2967 Domain: Euclidean Jordan algebra of dimension 2 over Algebraic
2968 Real Field (+) Euclidean Jordan algebra of dimension 3 over
2969 Algebraic Real Field
2970 Codomain: Euclidean Jordan algebra of dimension 2 over Algebraic
2972 sage: J.cartesian_projection(1)
2973 Linear operator between finite-dimensional Euclidean Jordan
2974 algebras represented by the matrix:
2978 Domain: Euclidean Jordan algebra of dimension 2 over Algebraic
2979 Real Field (+) Euclidean Jordan algebra of dimension 3 over
2980 Algebraic Real Field
2981 Codomain: Euclidean Jordan algebra of dimension 3 over Algebraic
2984 The projections work the way you'd expect on the vector
2985 representation of an element::
2987 sage: J1 = JordanSpinEJA(2)
2988 sage: J2 = ComplexHermitianEJA(2)
2989 sage: J = cartesian_product([J1,J2])
2990 sage: pi_left = J.cartesian_projection(0)
2991 sage: pi_right = J.cartesian_projection(1)
2992 sage: pi_left(J.one()).to_vector()
2994 sage: pi_right(J.one()).to_vector()
2996 sage: J.one().to_vector()
3001 The answer never changes::
3003 sage: set_random_seed()
3004 sage: J1 = random_eja()
3005 sage: J2 = random_eja()
3006 sage: J = cartesian_product([J1,J2])
3007 sage: P0 = J.cartesian_projection(0)
3008 sage: P1 = J.cartesian_projection(0)
3013 Ji
= self
.cartesian_factors()[i
]
3014 # Requires the fix on Trac 31421/31422 to work!
3015 Pi
= super().cartesian_projection(i
)
3016 return FiniteDimensionalEJAOperator(self
,Ji
,Pi
.matrix())
3019 def cartesian_embedding(self
, i
):
3023 sage: from mjo.eja.eja_algebra import (random_eja,
3024 ....: JordanSpinEJA,
3026 ....: RealSymmetricEJA)
3030 The embedding morphisms are Euclidean Jordan algebra
3033 sage: J1 = HadamardEJA(2)
3034 sage: J2 = RealSymmetricEJA(2)
3035 sage: J = cartesian_product([J1,J2])
3036 sage: J.cartesian_embedding(0)
3037 Linear operator between finite-dimensional Euclidean Jordan
3038 algebras represented by the matrix:
3044 Domain: Euclidean Jordan algebra of dimension 2 over
3045 Algebraic Real Field
3046 Codomain: Euclidean Jordan algebra of dimension 2 over
3047 Algebraic Real Field (+) Euclidean Jordan algebra of
3048 dimension 3 over Algebraic Real Field
3049 sage: J.cartesian_embedding(1)
3050 Linear operator between finite-dimensional Euclidean Jordan
3051 algebras represented by the matrix:
3057 Domain: Euclidean Jordan algebra of dimension 3 over
3058 Algebraic Real Field
3059 Codomain: Euclidean Jordan algebra of dimension 2 over
3060 Algebraic Real Field (+) Euclidean Jordan algebra of
3061 dimension 3 over Algebraic Real Field
3063 The embeddings work the way you'd expect on the vector
3064 representation of an element::
3066 sage: J1 = JordanSpinEJA(3)
3067 sage: J2 = RealSymmetricEJA(2)
3068 sage: J = cartesian_product([J1,J2])
3069 sage: iota_left = J.cartesian_embedding(0)
3070 sage: iota_right = J.cartesian_embedding(1)
3071 sage: iota_left(J1.zero()) == J.zero()
3073 sage: iota_right(J2.zero()) == J.zero()
3075 sage: J1.one().to_vector()
3077 sage: iota_left(J1.one()).to_vector()
3079 sage: J2.one().to_vector()
3081 sage: iota_right(J2.one()).to_vector()
3083 sage: J.one().to_vector()
3088 The answer never changes::
3090 sage: set_random_seed()
3091 sage: J1 = random_eja()
3092 sage: J2 = random_eja()
3093 sage: J = cartesian_product([J1,J2])
3094 sage: E0 = J.cartesian_embedding(0)
3095 sage: E1 = J.cartesian_embedding(0)
3099 Composing a projection with the corresponding inclusion should
3100 produce the identity map, and mismatching them should produce
3103 sage: set_random_seed()
3104 sage: J1 = random_eja()
3105 sage: J2 = random_eja()
3106 sage: J = cartesian_product([J1,J2])
3107 sage: iota_left = J.cartesian_embedding(0)
3108 sage: iota_right = J.cartesian_embedding(1)
3109 sage: pi_left = J.cartesian_projection(0)
3110 sage: pi_right = J.cartesian_projection(1)
3111 sage: pi_left*iota_left == J1.one().operator()
3113 sage: pi_right*iota_right == J2.one().operator()
3115 sage: (pi_left*iota_right).is_zero()
3117 sage: (pi_right*iota_left).is_zero()
3121 Ji
= self
.cartesian_factors()[i
]
3122 # Requires the fix on Trac 31421/31422 to work!
3123 Ei
= super().cartesian_embedding(i
)
3124 return FiniteDimensionalEJAOperator(Ji
,self
,Ei
.matrix())
3127 def cartesian_jordan_product(self
, x
, y
):
3129 The componentwise Jordan product.
3131 We project ``x`` and ``y`` onto our factors, and add up the
3132 Jordan products from the subalgebras. This may still be useful
3133 after (if) the default Jordan product in the Cartesian product
3134 algebra is overridden.
3138 sage: from mjo.eja.eja_algebra import (HadamardEJA,
3139 ....: JordanSpinEJA)
3143 sage: J1 = HadamardEJA(3)
3144 sage: J2 = JordanSpinEJA(3)
3145 sage: J = cartesian_product([J1,J2])
3146 sage: x1 = J1.from_vector(vector(QQ,(1,2,1)))
3147 sage: y1 = J1.from_vector(vector(QQ,(1,0,2)))
3148 sage: x2 = J2.from_vector(vector(QQ,(1,2,3)))
3149 sage: y2 = J2.from_vector(vector(QQ,(1,1,1)))
3150 sage: z1 = J.from_vector(vector(QQ,(1,2,1,1,2,3)))
3151 sage: z2 = J.from_vector(vector(QQ,(1,0,2,1,1,1)))
3152 sage: (x1*y1).to_vector()
3154 sage: (x2*y2).to_vector()
3156 sage: J.cartesian_jordan_product(z1,z2).to_vector()
3160 m
= len(self
.cartesian_factors())
3161 projections
= ( self
.cartesian_projection(i
) for i
in range(m
) )
3162 products
= ( P(x
)*P(y
) for P
in projections
)
3163 return self
._cartesian
_product
_of
_elements
(tuple(products
))
3165 def cartesian_inner_product(self
, x
, y
):
3167 The standard componentwise Cartesian inner-product.
3169 We project ``x`` and ``y`` onto our factors, and add up the
3170 inner-products from the subalgebras. This may still be useful
3171 after (if) the default inner product in the Cartesian product
3172 algebra is overridden.
3176 sage: from mjo.eja.eja_algebra import (HadamardEJA,
3177 ....: QuaternionHermitianEJA)
3181 sage: J1 = HadamardEJA(3,field=QQ)
3182 sage: J2 = QuaternionHermitianEJA(2,field=QQ,orthonormalize=False)
3183 sage: J = cartesian_product([J1,J2])
3188 sage: x1.inner_product(x2)
3190 sage: y1.inner_product(y2)
3192 sage: z1 = J._cartesian_product_of_elements((x1,y1))
3193 sage: z2 = J._cartesian_product_of_elements((x2,y2))
3194 sage: J.cartesian_inner_product(z1,z2)
3198 m
= len(self
.cartesian_factors())
3199 projections
= ( self
.cartesian_projection(i
) for i
in range(m
) )
3200 return sum( P(x
).inner_product(P(y
)) for P
in projections
)
3203 def _element_constructor_(self
, elt
):
3205 Construct an element of this algebra from an ordered tuple.
3207 We just apply the element constructor from each of our factors
3208 to the corresponding component of the tuple, and package up
3213 sage: from mjo.eja.eja_algebra import (HadamardEJA,
3214 ....: RealSymmetricEJA)
3218 sage: J1 = HadamardEJA(3)
3219 sage: J2 = RealSymmetricEJA(2)
3220 sage: J = cartesian_product([J1,J2])
3221 sage: J( (J1.matrix_basis()[1], J2.matrix_basis()[2]) )
3224 m
= len(self
.cartesian_factors())
3226 z
= tuple( self
.cartesian_factors()[i
](elt
[i
]) for i
in range(m
) )
3227 return self
._cartesian
_product
_of
_elements
(z
)
3229 raise ValueError("not an element of this algebra")
3231 Element
= CartesianProductEJAElement
3234 FiniteDimensionalEJA
.CartesianProduct
= CartesianProductEJA
3236 random_eja
= ConcreteEJA
.random_instance
3237 #def random_eja(*args, **kwargs):
3238 # from sage.categories.cartesian_product import cartesian_product
3239 # J1 = HadamardEJA(1, **kwargs)
3240 # J2 = RealSymmetricEJA(2, **kwargs)
3241 # J = cartesian_product([J1,J2])