2 Representations and constructions for Euclidean Jordan algebras.
4 A Euclidean Jordan algebra is a Jordan algebra that has some
7 1. It is finite-dimensional.
8 2. Its scalar field is the real numbers.
9 3a. An inner product is defined on it, and...
10 3b. That inner product is compatible with the Jordan product
11 in the sense that `<x*y,z> = <y,x*z>` for all elements
12 `x,y,z` in the algebra.
14 Every Euclidean Jordan algebra is formally-real: for any two elements
15 `x` and `y` in the algebra, `x^{2} + y^{2} = 0` implies that `x = y =
16 0`. Conversely, every finite-dimensional formally-real Jordan algebra
17 can be made into a Euclidean Jordan algebra with an appropriate choice
20 Formally-real Jordan algebras were originally studied as a framework
21 for quantum mechanics. Today, Euclidean Jordan algebras are crucial in
22 symmetric cone optimization, since every symmetric cone arises as the
23 cone of squares in some Euclidean Jordan algebra.
25 It is known that every Euclidean Jordan algebra decomposes into an
26 orthogonal direct sum (essentially, a Cartesian product) of simple
27 algebras, and that moreover, up to Jordan-algebra isomorphism, there
28 are only five families of simple algebras. We provide constructions
29 for these simple algebras:
31 * :class:`BilinearFormEJA`
32 * :class:`RealSymmetricEJA`
33 * :class:`ComplexHermitianEJA`
34 * :class:`QuaternionHermitianEJA`
36 Missing from this list is the algebra of three-by-three octononion
37 Hermitian matrices, as there is (as of yet) no implementation of the
38 octonions in SageMath. In addition to these, we provide two other
39 example constructions,
41 * :class:`HadamardEJA`
44 The Jordan spin algebra is a bilinear form algebra where the bilinear
45 form is the identity. The Hadamard EJA is simply a Cartesian product
46 of one-dimensional spin algebras. And last but not least, the trivial
47 EJA is exactly what you think. Cartesian products of these are also
48 supported using the usual ``cartesian_product()`` function; as a
49 result, we support (up to isomorphism) all Euclidean Jordan algebras
50 that don't involve octonions.
54 sage: from mjo.eja.eja_algebra import random_eja
59 Euclidean Jordan algebra of dimension...
62 from itertools
import repeat
64 from sage
.algebras
.quatalg
.quaternion_algebra
import QuaternionAlgebra
65 from sage
.categories
.magmatic_algebras
import MagmaticAlgebras
66 from sage
.categories
.sets_cat
import cartesian_product
67 from sage
.combinat
.free_module
import (CombinatorialFreeModule
,
68 CombinatorialFreeModule_CartesianProduct
)
69 from sage
.matrix
.constructor
import matrix
70 from sage
.matrix
.matrix_space
import MatrixSpace
71 from sage
.misc
.cachefunc
import cached_method
72 from sage
.misc
.table
import table
73 from sage
.modules
.free_module
import FreeModule
, VectorSpace
74 from sage
.rings
.all
import (ZZ
, QQ
, AA
, QQbar
, RR
, RLF
, CLF
,
77 from mjo
.eja
.eja_element
import FiniteDimensionalEJAElement
78 from mjo
.eja
.eja_operator
import FiniteDimensionalEJAOperator
79 from mjo
.eja
.eja_utils
import _all2list
, _mat2vec
81 class FiniteDimensionalEJA(CombinatorialFreeModule
):
83 A finite-dimensional Euclidean Jordan algebra.
87 - ``basis`` -- a tuple; a tuple of basis elements in "matrix
88 form," which must be the same form as the arguments to
89 ``jordan_product`` and ``inner_product``. In reality, "matrix
90 form" can be either vectors, matrices, or a Cartesian product
91 (ordered tuple) of vectors or matrices. All of these would
92 ideally be vector spaces in sage with no special-casing
93 needed; but in reality we turn vectors into column-matrices
94 and Cartesian products `(a,b)` into column matrices
95 `(a,b)^{T}` after converting `a` and `b` themselves.
97 - ``jordan_product`` -- a function; afunction of two ``basis``
98 elements (in matrix form) that returns their jordan product,
99 also in matrix form; this will be applied to ``basis`` to
100 compute a multiplication table for the algebra.
102 - ``inner_product`` -- a function; a function of two ``basis``
103 elements (in matrix form) that returns their inner
104 product. This will be applied to ``basis`` to compute an
105 inner-product table (basically a matrix) for this algebra.
107 - ``field`` -- a subfield of the reals (default: ``AA``); the scalar
108 field for the algebra.
110 - ``orthonormalize`` -- boolean (default: ``True``); whether or
111 not to orthonormalize the basis. Doing so is expensive and
112 generally rules out using the rationals as your ``field``, but
113 is required for spectral decompositions.
116 Element
= FiniteDimensionalEJAElement
125 cartesian_product
=False,
130 # Keep track of whether or not the matrix basis consists of
131 # tuples, since we need special cases for them damned near
132 # everywhere. This is INDEPENDENT of whether or not the
133 # algebra is a cartesian product, since a subalgebra of a
134 # cartesian product will have a basis of tuples, but will not
135 # in general itself be a cartesian product algebra.
136 self
._matrix
_basis
_is
_cartesian
= False
139 if hasattr(basis
[0], 'cartesian_factors'):
140 self
._matrix
_basis
_is
_cartesian
= True
143 if not field
.is_subring(RR
):
144 # Note: this does return true for the real algebraic
145 # field, the rationals, and any quadratic field where
146 # we've specified a real embedding.
147 raise ValueError("scalar field is not real")
149 # If the basis given to us wasn't over the field that it's
150 # supposed to be over, fix that. Or, you know, crash.
151 if not cartesian_product
:
152 # The field for a cartesian product algebra comes from one
153 # of its factors and is the same for all factors, so
154 # there's no need to "reapply" it on product algebras.
155 if self
._matrix
_basis
_is
_cartesian
:
156 # OK since if n == 0, the basis does not consist of tuples.
157 P
= basis
[0].parent()
158 basis
= tuple( P(tuple(b_i
.change_ring(field
) for b_i
in b
))
161 basis
= tuple( b
.change_ring(field
) for b
in basis
)
165 # Check commutativity of the Jordan and inner-products.
166 # This has to be done before we build the multiplication
167 # and inner-product tables/matrices, because we take
168 # advantage of symmetry in the process.
169 if not all( jordan_product(bi
,bj
) == jordan_product(bj
,bi
)
172 raise ValueError("Jordan product is not commutative")
174 if not all( inner_product(bi
,bj
) == inner_product(bj
,bi
)
177 raise ValueError("inner-product is not commutative")
180 category
= MagmaticAlgebras(field
).FiniteDimensional()
181 category
= category
.WithBasis().Unital().Commutative()
184 # Element subalgebras can take advantage of this.
185 category
= category
.Associative()
186 if cartesian_product
:
187 category
= category
.CartesianProducts()
189 # Call the superclass constructor so that we can use its from_vector()
190 # method to build our multiplication table.
191 CombinatorialFreeModule
.__init
__(self
,
198 # Now comes all of the hard work. We'll be constructing an
199 # ambient vector space V that our (vectorized) basis lives in,
200 # as well as a subspace W of V spanned by those (vectorized)
201 # basis elements. The W-coordinates are the coefficients that
202 # we see in things like x = 1*e1 + 2*e2.
207 degree
= len(_all2list(basis
[0]))
209 # Build an ambient space that fits our matrix basis when
210 # written out as "long vectors."
211 V
= VectorSpace(field
, degree
)
213 # The matrix that will hole the orthonormal -> unorthonormal
214 # coordinate transformation.
215 self
._deortho
_matrix
= None
218 # Save a copy of the un-orthonormalized basis for later.
219 # Convert it to ambient V (vector) coordinates while we're
220 # at it, because we'd have to do it later anyway.
221 deortho_vector_basis
= tuple( V(_all2list(b
)) for b
in basis
)
223 from mjo
.eja
.eja_utils
import gram_schmidt
224 basis
= tuple(gram_schmidt(basis
, inner_product
))
226 # Save the (possibly orthonormalized) matrix basis for
228 self
._matrix
_basis
= basis
230 # Now create the vector space for the algebra, which will have
231 # its own set of non-ambient coordinates (in terms of the
233 vector_basis
= tuple( V(_all2list(b
)) for b
in basis
)
234 W
= V
.span_of_basis( vector_basis
, check
=check_axioms
)
237 # Now "W" is the vector space of our algebra coordinates. The
238 # variables "X1", "X2",... refer to the entries of vectors in
239 # W. Thus to convert back and forth between the orthonormal
240 # coordinates and the given ones, we need to stick the original
242 U
= V
.span_of_basis( deortho_vector_basis
, check
=check_axioms
)
243 self
._deortho
_matrix
= matrix( U
.coordinate_vector(q
)
244 for q
in vector_basis
)
247 # Now we actually compute the multiplication and inner-product
248 # tables/matrices using the possibly-orthonormalized basis.
249 self
._inner
_product
_matrix
= matrix
.identity(field
, n
)
250 self
._multiplication
_table
= [ [0 for j
in range(i
+1)]
253 # Note: the Jordan and inner-products are defined in terms
254 # of the ambient basis. It's important that their arguments
255 # are in ambient coordinates as well.
258 # ortho basis w.r.t. ambient coords
262 # The jordan product returns a matrixy answer, so we
263 # have to convert it to the algebra coordinates.
264 elt
= jordan_product(q_i
, q_j
)
265 elt
= W
.coordinate_vector(V(_all2list(elt
)))
266 self
._multiplication
_table
[i
][j
] = self
.from_vector(elt
)
268 if not orthonormalize
:
269 # If we're orthonormalizing the basis with respect
270 # to an inner-product, then the inner-product
271 # matrix with respect to the resulting basis is
272 # just going to be the identity.
273 ip
= inner_product(q_i
, q_j
)
274 self
._inner
_product
_matrix
[i
,j
] = ip
275 self
._inner
_product
_matrix
[j
,i
] = ip
277 self
._inner
_product
_matrix
._cache
= {'hermitian': True}
278 self
._inner
_product
_matrix
.set_immutable()
281 if not self
._is
_jordanian
():
282 raise ValueError("Jordan identity does not hold")
283 if not self
._inner
_product
_is
_associative
():
284 raise ValueError("inner product is not associative")
287 def _coerce_map_from_base_ring(self
):
289 Disable the map from the base ring into the algebra.
291 Performing a nonsense conversion like this automatically
292 is counterpedagogical. The fallback is to try the usual
293 element constructor, which should also fail.
297 sage: from mjo.eja.eja_algebra import random_eja
301 sage: set_random_seed()
302 sage: J = random_eja()
304 Traceback (most recent call last):
306 ValueError: not an element of this algebra
312 def product_on_basis(self
, i
, j
):
314 Returns the Jordan product of the `i` and `j`th basis elements.
316 This completely defines the Jordan product on the algebra, and
317 is used direclty by our superclass machinery to implement
322 sage: from mjo.eja.eja_algebra import random_eja
326 sage: set_random_seed()
327 sage: J = random_eja()
328 sage: n = J.dimension()
331 sage: ei_ej = J.zero()*J.zero()
333 ....: i = ZZ.random_element(n)
334 ....: j = ZZ.random_element(n)
335 ....: ei = J.gens()[i]
336 ....: ej = J.gens()[j]
337 ....: ei_ej = J.product_on_basis(i,j)
342 # We only stored the lower-triangular portion of the
343 # multiplication table.
345 return self
._multiplication
_table
[i
][j
]
347 return self
._multiplication
_table
[j
][i
]
349 def inner_product(self
, x
, y
):
351 The inner product associated with this Euclidean Jordan algebra.
353 Defaults to the trace inner product, but can be overridden by
354 subclasses if they are sure that the necessary properties are
359 sage: from mjo.eja.eja_algebra import (random_eja,
361 ....: BilinearFormEJA)
365 Our inner product is "associative," which means the following for
366 a symmetric bilinear form::
368 sage: set_random_seed()
369 sage: J = random_eja()
370 sage: x,y,z = J.random_elements(3)
371 sage: (x*y).inner_product(z) == y.inner_product(x*z)
376 Ensure that this is the usual inner product for the algebras
379 sage: set_random_seed()
380 sage: J = HadamardEJA.random_instance()
381 sage: x,y = J.random_elements(2)
382 sage: actual = x.inner_product(y)
383 sage: expected = x.to_vector().inner_product(y.to_vector())
384 sage: actual == expected
387 Ensure that this is one-half of the trace inner-product in a
388 BilinearFormEJA that isn't just the reals (when ``n`` isn't
389 one). This is in Faraut and Koranyi, and also my "On the
392 sage: set_random_seed()
393 sage: J = BilinearFormEJA.random_instance()
394 sage: n = J.dimension()
395 sage: x = J.random_element()
396 sage: y = J.random_element()
397 sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
401 B
= self
._inner
_product
_matrix
402 return (B
*x
.to_vector()).inner_product(y
.to_vector())
405 def is_associative(self
):
407 Return whether or not this algebra's Jordan product is associative.
411 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
415 sage: J = ComplexHermitianEJA(3, field=QQ, orthonormalize=False)
416 sage: J.is_associative()
418 sage: x = sum(J.gens())
419 sage: A = x.subalgebra_generated_by(orthonormalize=False)
420 sage: A.is_associative()
424 return "Associative" in self
.category().axioms()
426 def _is_commutative(self
):
428 Whether or not this algebra's multiplication table is commutative.
430 This method should of course always return ``True``, unless
431 this algebra was constructed with ``check_axioms=False`` and
432 passed an invalid multiplication table.
434 return all( self
.product_on_basis(i
,j
) == self
.product_on_basis(i
,j
)
435 for i
in range(self
.dimension())
436 for j
in range(self
.dimension()) )
438 def _is_jordanian(self
):
440 Whether or not this algebra's multiplication table respects the
441 Jordan identity `(x^{2})(xy) = x(x^{2}y)`.
443 We only check one arrangement of `x` and `y`, so for a
444 ``True`` result to be truly true, you should also check
445 :meth:`_is_commutative`. This method should of course always
446 return ``True``, unless this algebra was constructed with
447 ``check_axioms=False`` and passed an invalid multiplication table.
449 return all( (self
.gens()[i
]**2)*(self
.gens()[i
]*self
.gens()[j
])
451 (self
.gens()[i
])*((self
.gens()[i
]**2)*self
.gens()[j
])
452 for i
in range(self
.dimension())
453 for j
in range(self
.dimension()) )
455 def _inner_product_is_associative(self
):
457 Return whether or not this algebra's inner product `B` is
458 associative; that is, whether or not `B(xy,z) = B(x,yz)`.
460 This method should of course always return ``True``, unless
461 this algebra was constructed with ``check_axioms=False`` and
462 passed an invalid Jordan or inner-product.
466 # Used to check whether or not something is zero.
469 # This choice is sufficient to allow the construction of
470 # QuaternionHermitianEJA(2, field=RDF) with check_axioms=True.
473 for i
in range(self
.dimension()):
474 for j
in range(self
.dimension()):
475 for k
in range(self
.dimension()):
479 diff
= (x
*y
).inner_product(z
) - x
.inner_product(y
*z
)
481 if diff
.abs() > epsilon
:
486 def _element_constructor_(self
, elt
):
488 Construct an element of this algebra from its vector or matrix
491 This gets called only after the parent element _call_ method
492 fails to find a coercion for the argument.
496 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
498 ....: RealSymmetricEJA)
502 The identity in `S^n` is converted to the identity in the EJA::
504 sage: J = RealSymmetricEJA(3)
505 sage: I = matrix.identity(QQ,3)
506 sage: J(I) == J.one()
509 This skew-symmetric matrix can't be represented in the EJA::
511 sage: J = RealSymmetricEJA(3)
512 sage: A = matrix(QQ,3, lambda i,j: i-j)
514 Traceback (most recent call last):
516 ValueError: not an element of this algebra
518 Tuples work as well, provided that the matrix basis for the
519 algebra consists of them::
521 sage: J1 = HadamardEJA(3)
522 sage: J2 = RealSymmetricEJA(2)
523 sage: J = cartesian_product([J1,J2])
524 sage: J( (J1.matrix_basis()[1], J2.matrix_basis()[2]) )
529 Ensure that we can convert any element of the two non-matrix
530 simple algebras (whose matrix representations are columns)
531 back and forth faithfully::
533 sage: set_random_seed()
534 sage: J = HadamardEJA.random_instance()
535 sage: x = J.random_element()
536 sage: J(x.to_vector().column()) == x
538 sage: J = JordanSpinEJA.random_instance()
539 sage: x = J.random_element()
540 sage: J(x.to_vector().column()) == x
543 We cannot coerce elements between algebras just because their
544 matrix representations are compatible::
546 sage: J1 = HadamardEJA(3)
547 sage: J2 = JordanSpinEJA(3)
549 Traceback (most recent call last):
551 ValueError: not an element of this algebra
553 Traceback (most recent call last):
555 ValueError: not an element of this algebra
558 msg
= "not an element of this algebra"
559 if elt
in self
.base_ring():
560 # Ensure that no base ring -> algebra coercion is performed
561 # by this method. There's some stupidity in sage that would
562 # otherwise propagate to this method; for example, sage thinks
563 # that the integer 3 belongs to the space of 2-by-2 matrices.
564 raise ValueError(msg
)
567 # Try to convert a vector into a column-matrix...
569 except (AttributeError, TypeError):
570 # and ignore failure, because we weren't really expecting
571 # a vector as an argument anyway.
574 if elt
not in self
.matrix_space():
575 raise ValueError(msg
)
577 # Thanks for nothing! Matrix spaces aren't vector spaces in
578 # Sage, so we have to figure out its matrix-basis coordinates
579 # ourselves. We use the basis space's ring instead of the
580 # element's ring because the basis space might be an algebraic
581 # closure whereas the base ring of the 3-by-3 identity matrix
582 # could be QQ instead of QQbar.
584 # And, we also have to handle Cartesian product bases (when
585 # the matrix basis consists of tuples) here. The "good news"
586 # is that we're already converting everything to long vectors,
587 # and that strategy works for tuples as well.
589 # We pass check=False because the matrix basis is "guaranteed"
590 # to be linearly independent... right? Ha ha.
592 V
= VectorSpace(self
.base_ring(), len(elt
))
593 W
= V
.span_of_basis( (V(_all2list(s
)) for s
in self
.matrix_basis()),
597 coords
= W
.coordinate_vector(V(elt
))
598 except ArithmeticError: # vector is not in free module
599 raise ValueError(msg
)
601 return self
.from_vector(coords
)
605 Return a string representation of ``self``.
609 sage: from mjo.eja.eja_algebra import JordanSpinEJA
613 Ensure that it says what we think it says::
615 sage: JordanSpinEJA(2, field=AA)
616 Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
617 sage: JordanSpinEJA(3, field=RDF)
618 Euclidean Jordan algebra of dimension 3 over Real Double Field
621 fmt
= "Euclidean Jordan algebra of dimension {} over {}"
622 return fmt
.format(self
.dimension(), self
.base_ring())
626 def characteristic_polynomial_of(self
):
628 Return the algebra's "characteristic polynomial of" function,
629 which is itself a multivariate polynomial that, when evaluated
630 at the coordinates of some algebra element, returns that
631 element's characteristic polynomial.
633 The resulting polynomial has `n+1` variables, where `n` is the
634 dimension of this algebra. The first `n` variables correspond to
635 the coordinates of an algebra element: when evaluated at the
636 coordinates of an algebra element with respect to a certain
637 basis, the result is a univariate polynomial (in the one
638 remaining variable ``t``), namely the characteristic polynomial
643 sage: from mjo.eja.eja_algebra import JordanSpinEJA, TrivialEJA
647 The characteristic polynomial in the spin algebra is given in
648 Alizadeh, Example 11.11::
650 sage: J = JordanSpinEJA(3)
651 sage: p = J.characteristic_polynomial_of(); p
652 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
653 sage: xvec = J.one().to_vector()
657 By definition, the characteristic polynomial is a monic
658 degree-zero polynomial in a rank-zero algebra. Note that
659 Cayley-Hamilton is indeed satisfied since the polynomial
660 ``1`` evaluates to the identity element of the algebra on
663 sage: J = TrivialEJA()
664 sage: J.characteristic_polynomial_of()
671 # The list of coefficient polynomials a_0, a_1, a_2, ..., a_(r-1).
672 a
= self
._charpoly
_coefficients
()
674 # We go to a bit of trouble here to reorder the
675 # indeterminates, so that it's easier to evaluate the
676 # characteristic polynomial at x's coordinates and get back
677 # something in terms of t, which is what we want.
678 S
= PolynomialRing(self
.base_ring(),'t')
682 S
= PolynomialRing(S
, R
.variable_names())
685 return (t
**r
+ sum( a
[k
]*(t
**k
) for k
in range(r
) ))
687 def coordinate_polynomial_ring(self
):
689 The multivariate polynomial ring in which this algebra's
690 :meth:`characteristic_polynomial_of` lives.
694 sage: from mjo.eja.eja_algebra import (HadamardEJA,
695 ....: RealSymmetricEJA)
699 sage: J = HadamardEJA(2)
700 sage: J.coordinate_polynomial_ring()
701 Multivariate Polynomial Ring in X1, X2...
702 sage: J = RealSymmetricEJA(3,field=QQ,orthonormalize=False)
703 sage: J.coordinate_polynomial_ring()
704 Multivariate Polynomial Ring in X1, X2, X3, X4, X5, X6...
707 var_names
= tuple( "X%d" % z
for z
in range(1, self
.dimension()+1) )
708 return PolynomialRing(self
.base_ring(), var_names
)
710 def inner_product(self
, x
, y
):
712 The inner product associated with this Euclidean Jordan algebra.
714 Defaults to the trace inner product, but can be overridden by
715 subclasses if they are sure that the necessary properties are
720 sage: from mjo.eja.eja_algebra import (random_eja,
722 ....: BilinearFormEJA)
726 Our inner product is "associative," which means the following for
727 a symmetric bilinear form::
729 sage: set_random_seed()
730 sage: J = random_eja()
731 sage: x,y,z = J.random_elements(3)
732 sage: (x*y).inner_product(z) == y.inner_product(x*z)
737 Ensure that this is the usual inner product for the algebras
740 sage: set_random_seed()
741 sage: J = HadamardEJA.random_instance()
742 sage: x,y = J.random_elements(2)
743 sage: actual = x.inner_product(y)
744 sage: expected = x.to_vector().inner_product(y.to_vector())
745 sage: actual == expected
748 Ensure that this is one-half of the trace inner-product in a
749 BilinearFormEJA that isn't just the reals (when ``n`` isn't
750 one). This is in Faraut and Koranyi, and also my "On the
753 sage: set_random_seed()
754 sage: J = BilinearFormEJA.random_instance()
755 sage: n = J.dimension()
756 sage: x = J.random_element()
757 sage: y = J.random_element()
758 sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
761 B
= self
._inner
_product
_matrix
762 return (B
*x
.to_vector()).inner_product(y
.to_vector())
765 def is_trivial(self
):
767 Return whether or not this algebra is trivial.
769 A trivial algebra contains only the zero element.
773 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
778 sage: J = ComplexHermitianEJA(3)
784 sage: J = TrivialEJA()
789 return self
.dimension() == 0
792 def multiplication_table(self
):
794 Return a visual representation of this algebra's multiplication
795 table (on basis elements).
799 sage: from mjo.eja.eja_algebra import JordanSpinEJA
803 sage: J = JordanSpinEJA(4)
804 sage: J.multiplication_table()
805 +----++----+----+----+----+
806 | * || e0 | e1 | e2 | e3 |
807 +====++====+====+====+====+
808 | e0 || e0 | e1 | e2 | e3 |
809 +----++----+----+----+----+
810 | e1 || e1 | e0 | 0 | 0 |
811 +----++----+----+----+----+
812 | e2 || e2 | 0 | e0 | 0 |
813 +----++----+----+----+----+
814 | e3 || e3 | 0 | 0 | e0 |
815 +----++----+----+----+----+
819 # Prepend the header row.
820 M
= [["*"] + list(self
.gens())]
822 # And to each subsequent row, prepend an entry that belongs to
823 # the left-side "header column."
824 M
+= [ [self
.gens()[i
]] + [ self
.product_on_basis(i
,j
)
828 return table(M
, header_row
=True, header_column
=True, frame
=True)
831 def matrix_basis(self
):
833 Return an (often more natural) representation of this algebras
834 basis as an ordered tuple of matrices.
836 Every finite-dimensional Euclidean Jordan Algebra is a, up to
837 Jordan isomorphism, a direct sum of five simple
838 algebras---four of which comprise Hermitian matrices. And the
839 last type of algebra can of course be thought of as `n`-by-`1`
840 column matrices (ambiguusly called column vectors) to avoid
841 special cases. As a result, matrices (and column vectors) are
842 a natural representation format for Euclidean Jordan algebra
845 But, when we construct an algebra from a basis of matrices,
846 those matrix representations are lost in favor of coordinate
847 vectors *with respect to* that basis. We could eventually
848 convert back if we tried hard enough, but having the original
849 representations handy is valuable enough that we simply store
850 them and return them from this method.
852 Why implement this for non-matrix algebras? Avoiding special
853 cases for the :class:`BilinearFormEJA` pays with simplicity in
854 its own right. But mainly, we would like to be able to assume
855 that elements of a :class:`CartesianProductEJA` can be displayed
856 nicely, without having to have special classes for direct sums
857 one of whose components was a matrix algebra.
861 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
862 ....: RealSymmetricEJA)
866 sage: J = RealSymmetricEJA(2)
868 Finite family {0: e0, 1: e1, 2: e2}
869 sage: J.matrix_basis()
871 [1 0] [ 0 0.7071067811865475?] [0 0]
872 [0 0], [0.7071067811865475? 0], [0 1]
877 sage: J = JordanSpinEJA(2)
879 Finite family {0: e0, 1: e1}
880 sage: J.matrix_basis()
886 return self
._matrix
_basis
889 def matrix_space(self
):
891 Return the matrix space in which this algebra's elements live, if
892 we think of them as matrices (including column vectors of the
895 "By default" this will be an `n`-by-`1` column-matrix space,
896 except when the algebra is trivial. There it's `n`-by-`n`
897 (where `n` is zero), to ensure that two elements of the matrix
898 space (empty matrices) can be multiplied. For algebras of
899 matrices, this returns the space in which their
900 real embeddings live.
904 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
906 ....: QuaternionHermitianEJA,
911 By default, the matrix representation is just a column-matrix
912 equivalent to the vector representation::
914 sage: J = JordanSpinEJA(3)
915 sage: J.matrix_space()
916 Full MatrixSpace of 3 by 1 dense matrices over Algebraic
919 The matrix representation in the trivial algebra is
920 zero-by-zero instead of the usual `n`-by-one::
922 sage: J = TrivialEJA()
923 sage: J.matrix_space()
924 Full MatrixSpace of 0 by 0 dense matrices over Algebraic
927 The matrix space for complex/quaternion Hermitian matrix EJA
928 is the space in which their real-embeddings live, not the
929 original complex/quaternion matrix space::
931 sage: J = ComplexHermitianEJA(2,field=QQ,orthonormalize=False)
932 sage: J.matrix_space()
933 Full MatrixSpace of 4 by 4 dense matrices over Rational Field
934 sage: J = QuaternionHermitianEJA(1,field=QQ,orthonormalize=False)
935 sage: J.matrix_space()
936 Full MatrixSpace of 4 by 4 dense matrices over Rational Field
939 if self
.is_trivial():
940 return MatrixSpace(self
.base_ring(), 0)
942 return self
.matrix_basis()[0].parent()
948 Return the unit element of this algebra.
952 sage: from mjo.eja.eja_algebra import (HadamardEJA,
957 We can compute unit element in the Hadamard EJA::
959 sage: J = HadamardEJA(5)
961 e0 + e1 + e2 + e3 + e4
963 The unit element in the Hadamard EJA is inherited in the
964 subalgebras generated by its elements::
966 sage: J = HadamardEJA(5)
968 e0 + e1 + e2 + e3 + e4
969 sage: x = sum(J.gens())
970 sage: A = x.subalgebra_generated_by(orthonormalize=False)
973 sage: A.one().superalgebra_element()
974 e0 + e1 + e2 + e3 + e4
978 The identity element acts like the identity, regardless of
979 whether or not we orthonormalize::
981 sage: set_random_seed()
982 sage: J = random_eja()
983 sage: x = J.random_element()
984 sage: J.one()*x == x and x*J.one() == x
986 sage: A = x.subalgebra_generated_by()
987 sage: y = A.random_element()
988 sage: A.one()*y == y and y*A.one() == y
993 sage: set_random_seed()
994 sage: J = random_eja(field=QQ, orthonormalize=False)
995 sage: x = J.random_element()
996 sage: J.one()*x == x and x*J.one() == x
998 sage: A = x.subalgebra_generated_by(orthonormalize=False)
999 sage: y = A.random_element()
1000 sage: A.one()*y == y and y*A.one() == y
1003 The matrix of the unit element's operator is the identity,
1004 regardless of the base field and whether or not we
1007 sage: set_random_seed()
1008 sage: J = random_eja()
1009 sage: actual = J.one().operator().matrix()
1010 sage: expected = matrix.identity(J.base_ring(), J.dimension())
1011 sage: actual == expected
1013 sage: x = J.random_element()
1014 sage: A = x.subalgebra_generated_by()
1015 sage: actual = A.one().operator().matrix()
1016 sage: expected = matrix.identity(A.base_ring(), A.dimension())
1017 sage: actual == expected
1022 sage: set_random_seed()
1023 sage: J = random_eja(field=QQ, orthonormalize=False)
1024 sage: actual = J.one().operator().matrix()
1025 sage: expected = matrix.identity(J.base_ring(), J.dimension())
1026 sage: actual == expected
1028 sage: x = J.random_element()
1029 sage: A = x.subalgebra_generated_by(orthonormalize=False)
1030 sage: actual = A.one().operator().matrix()
1031 sage: expected = matrix.identity(A.base_ring(), A.dimension())
1032 sage: actual == expected
1035 Ensure that the cached unit element (often precomputed by
1036 hand) agrees with the computed one::
1038 sage: set_random_seed()
1039 sage: J = random_eja()
1040 sage: cached = J.one()
1041 sage: J.one.clear_cache()
1042 sage: J.one() == cached
1047 sage: set_random_seed()
1048 sage: J = random_eja(field=QQ, orthonormalize=False)
1049 sage: cached = J.one()
1050 sage: J.one.clear_cache()
1051 sage: J.one() == cached
1055 # We can brute-force compute the matrices of the operators
1056 # that correspond to the basis elements of this algebra.
1057 # If some linear combination of those basis elements is the
1058 # algebra identity, then the same linear combination of
1059 # their matrices has to be the identity matrix.
1061 # Of course, matrices aren't vectors in sage, so we have to
1062 # appeal to the "long vectors" isometry.
1063 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
1065 # Now we use basic linear algebra to find the coefficients,
1066 # of the matrices-as-vectors-linear-combination, which should
1067 # work for the original algebra basis too.
1068 A
= matrix(self
.base_ring(), oper_vecs
)
1070 # We used the isometry on the left-hand side already, but we
1071 # still need to do it for the right-hand side. Recall that we
1072 # wanted something that summed to the identity matrix.
1073 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
1075 # Now if there's an identity element in the algebra, this
1076 # should work. We solve on the left to avoid having to
1077 # transpose the matrix "A".
1078 return self
.from_vector(A
.solve_left(b
))
1081 def peirce_decomposition(self
, c
):
1083 The Peirce decomposition of this algebra relative to the
1086 In the future, this can be extended to a complete system of
1087 orthogonal idempotents.
1091 - ``c`` -- an idempotent of this algebra.
1095 A triple (J0, J5, J1) containing two subalgebras and one subspace
1098 - ``J0`` -- the algebra on the eigenspace of ``c.operator()``
1099 corresponding to the eigenvalue zero.
1101 - ``J5`` -- the eigenspace (NOT a subalgebra) of ``c.operator()``
1102 corresponding to the eigenvalue one-half.
1104 - ``J1`` -- the algebra on the eigenspace of ``c.operator()``
1105 corresponding to the eigenvalue one.
1107 These are the only possible eigenspaces for that operator, and this
1108 algebra is a direct sum of them. The spaces ``J0`` and ``J1`` are
1109 orthogonal, and are subalgebras of this algebra with the appropriate
1114 sage: from mjo.eja.eja_algebra import random_eja, RealSymmetricEJA
1118 The canonical example comes from the symmetric matrices, which
1119 decompose into diagonal and off-diagonal parts::
1121 sage: J = RealSymmetricEJA(3)
1122 sage: C = matrix(QQ, [ [1,0,0],
1126 sage: J0,J5,J1 = J.peirce_decomposition(c)
1128 Euclidean Jordan algebra of dimension 1...
1130 Vector space of degree 6 and dimension 2...
1132 Euclidean Jordan algebra of dimension 3...
1133 sage: J0.one().to_matrix()
1137 sage: orig_df = AA.options.display_format
1138 sage: AA.options.display_format = 'radical'
1139 sage: J.from_vector(J5.basis()[0]).to_matrix()
1143 sage: J.from_vector(J5.basis()[1]).to_matrix()
1147 sage: AA.options.display_format = orig_df
1148 sage: J1.one().to_matrix()
1155 Every algebra decomposes trivially with respect to its identity
1158 sage: set_random_seed()
1159 sage: J = random_eja()
1160 sage: J0,J5,J1 = J.peirce_decomposition(J.one())
1161 sage: J0.dimension() == 0 and J5.dimension() == 0
1163 sage: J1.superalgebra() == J and J1.dimension() == J.dimension()
1166 The decomposition is into eigenspaces, and its components are
1167 therefore necessarily orthogonal. Moreover, the identity
1168 elements in the two subalgebras are the projections onto their
1169 respective subspaces of the superalgebra's identity element::
1171 sage: set_random_seed()
1172 sage: J = random_eja()
1173 sage: x = J.random_element()
1174 sage: if not J.is_trivial():
1175 ....: while x.is_nilpotent():
1176 ....: x = J.random_element()
1177 sage: c = x.subalgebra_idempotent()
1178 sage: J0,J5,J1 = J.peirce_decomposition(c)
1180 sage: for (w,y,z) in zip(J0.basis(), J5.basis(), J1.basis()):
1181 ....: w = w.superalgebra_element()
1182 ....: y = J.from_vector(y)
1183 ....: z = z.superalgebra_element()
1184 ....: ipsum += w.inner_product(y).abs()
1185 ....: ipsum += w.inner_product(z).abs()
1186 ....: ipsum += y.inner_product(z).abs()
1189 sage: J1(c) == J1.one()
1191 sage: J0(J.one() - c) == J0.one()
1195 if not c
.is_idempotent():
1196 raise ValueError("element is not idempotent: %s" % c
)
1198 # Default these to what they should be if they turn out to be
1199 # trivial, because eigenspaces_left() won't return eigenvalues
1200 # corresponding to trivial spaces (e.g. it returns only the
1201 # eigenspace corresponding to lambda=1 if you take the
1202 # decomposition relative to the identity element).
1203 trivial
= self
.subalgebra(())
1204 J0
= trivial
# eigenvalue zero
1205 J5
= VectorSpace(self
.base_ring(), 0) # eigenvalue one-half
1206 J1
= trivial
# eigenvalue one
1208 for (eigval
, eigspace
) in c
.operator().matrix().right_eigenspaces():
1209 if eigval
== ~
(self
.base_ring()(2)):
1212 gens
= tuple( self
.from_vector(b
) for b
in eigspace
.basis() )
1213 subalg
= self
.subalgebra(gens
, check_axioms
=False)
1219 raise ValueError("unexpected eigenvalue: %s" % eigval
)
1224 def random_element(self
, thorough
=False):
1226 Return a random element of this algebra.
1228 Our algebra superclass method only returns a linear
1229 combination of at most two basis elements. We instead
1230 want the vector space "random element" method that
1231 returns a more diverse selection.
1235 - ``thorough`` -- (boolean; default False) whether or not we
1236 should generate irrational coefficients for the random
1237 element when our base ring is irrational; this slows the
1238 algebra operations to a crawl, but any truly random method
1242 # For a general base ring... maybe we can trust this to do the
1243 # right thing? Unlikely, but.
1244 V
= self
.vector_space()
1245 v
= V
.random_element()
1247 if self
.base_ring() is AA
:
1248 # The "random element" method of the algebraic reals is
1249 # stupid at the moment, and only returns integers between
1250 # -2 and 2, inclusive:
1252 # https://trac.sagemath.org/ticket/30875
1254 # Instead, we implement our own "random vector" method,
1255 # and then coerce that into the algebra. We use the vector
1256 # space degree here instead of the dimension because a
1257 # subalgebra could (for example) be spanned by only two
1258 # vectors, each with five coordinates. We need to
1259 # generate all five coordinates.
1261 v
*= QQbar
.random_element().real()
1263 v
*= QQ
.random_element()
1265 return self
.from_vector(V
.coordinate_vector(v
))
1267 def random_elements(self
, count
, thorough
=False):
1269 Return ``count`` random elements as a tuple.
1273 - ``thorough`` -- (boolean; default False) whether or not we
1274 should generate irrational coefficients for the random
1275 elements when our base ring is irrational; this slows the
1276 algebra operations to a crawl, but any truly random method
1281 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1285 sage: J = JordanSpinEJA(3)
1286 sage: x,y,z = J.random_elements(3)
1287 sage: all( [ x in J, y in J, z in J ])
1289 sage: len( J.random_elements(10) ) == 10
1293 return tuple( self
.random_element(thorough
)
1294 for idx
in range(count
) )
1298 def _charpoly_coefficients(self
):
1300 The `r` polynomial coefficients of the "characteristic polynomial
1305 sage: from mjo.eja.eja_algebra import random_eja
1309 The theory shows that these are all homogeneous polynomials of
1312 sage: set_random_seed()
1313 sage: J = random_eja()
1314 sage: all(p.is_homogeneous() for p in J._charpoly_coefficients())
1318 n
= self
.dimension()
1319 R
= self
.coordinate_polynomial_ring()
1321 F
= R
.fraction_field()
1324 # From a result in my book, these are the entries of the
1325 # basis representation of L_x.
1326 return sum( vars[k
]*self
.gens()[k
].operator().matrix()[i
,j
]
1329 L_x
= matrix(F
, n
, n
, L_x_i_j
)
1332 if self
.rank
.is_in_cache():
1334 # There's no need to pad the system with redundant
1335 # columns if we *know* they'll be redundant.
1338 # Compute an extra power in case the rank is equal to
1339 # the dimension (otherwise, we would stop at x^(r-1)).
1340 x_powers
= [ (L_x
**k
)*self
.one().to_vector()
1341 for k
in range(n
+1) ]
1342 A
= matrix
.column(F
, x_powers
[:n
])
1343 AE
= A
.extended_echelon_form()
1350 # The theory says that only the first "r" coefficients are
1351 # nonzero, and they actually live in the original polynomial
1352 # ring and not the fraction field. We negate them because in
1353 # the actual characteristic polynomial, they get moved to the
1354 # other side where x^r lives. We don't bother to trim A_rref
1355 # down to a square matrix and solve the resulting system,
1356 # because the upper-left r-by-r portion of A_rref is
1357 # guaranteed to be the identity matrix, so e.g.
1359 # A_rref.solve_right(Y)
1361 # would just be returning Y.
1362 return (-E
*b
)[:r
].change_ring(R
)
1367 Return the rank of this EJA.
1369 This is a cached method because we know the rank a priori for
1370 all of the algebras we can construct. Thus we can avoid the
1371 expensive ``_charpoly_coefficients()`` call unless we truly
1372 need to compute the whole characteristic polynomial.
1376 sage: from mjo.eja.eja_algebra import (HadamardEJA,
1377 ....: JordanSpinEJA,
1378 ....: RealSymmetricEJA,
1379 ....: ComplexHermitianEJA,
1380 ....: QuaternionHermitianEJA,
1385 The rank of the Jordan spin algebra is always two::
1387 sage: JordanSpinEJA(2).rank()
1389 sage: JordanSpinEJA(3).rank()
1391 sage: JordanSpinEJA(4).rank()
1394 The rank of the `n`-by-`n` Hermitian real, complex, or
1395 quaternion matrices is `n`::
1397 sage: RealSymmetricEJA(4).rank()
1399 sage: ComplexHermitianEJA(3).rank()
1401 sage: QuaternionHermitianEJA(2).rank()
1406 Ensure that every EJA that we know how to construct has a
1407 positive integer rank, unless the algebra is trivial in
1408 which case its rank will be zero::
1410 sage: set_random_seed()
1411 sage: J = random_eja()
1415 sage: r > 0 or (r == 0 and J.is_trivial())
1418 Ensure that computing the rank actually works, since the ranks
1419 of all simple algebras are known and will be cached by default::
1421 sage: set_random_seed() # long time
1422 sage: J = random_eja() # long time
1423 sage: cached = J.rank() # long time
1424 sage: J.rank.clear_cache() # long time
1425 sage: J.rank() == cached # long time
1429 return len(self
._charpoly
_coefficients
())
1432 def subalgebra(self
, basis
, **kwargs
):
1434 Create a subalgebra of this algebra from the given basis.
1436 from mjo
.eja
.eja_subalgebra
import FiniteDimensionalEJASubalgebra
1437 return FiniteDimensionalEJASubalgebra(self
, basis
, **kwargs
)
1440 def vector_space(self
):
1442 Return the vector space that underlies this algebra.
1446 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1450 sage: J = RealSymmetricEJA(2)
1451 sage: J.vector_space()
1452 Vector space of dimension 3 over...
1455 return self
.zero().to_vector().parent().ambient_vector_space()
1459 class RationalBasisEJA(FiniteDimensionalEJA
):
1461 New class for algebras whose supplied basis elements have all rational entries.
1465 sage: from mjo.eja.eja_algebra import BilinearFormEJA
1469 The supplied basis is orthonormalized by default::
1471 sage: B = matrix(QQ, [[1, 0, 0], [0, 25, -32], [0, -32, 41]])
1472 sage: J = BilinearFormEJA(B)
1473 sage: J.matrix_basis()
1490 # Abuse the check_field parameter to check that the entries of
1491 # out basis (in ambient coordinates) are in the field QQ.
1492 if not all( all(b_i
in QQ
for b_i
in b
.list()) for b
in basis
):
1493 raise TypeError("basis not rational")
1495 self
._rational
_algebra
= None
1497 # There's no point in constructing the extra algebra if this
1498 # one is already rational.
1500 # Note: the same Jordan and inner-products work here,
1501 # because they are necessarily defined with respect to
1502 # ambient coordinates and not any particular basis.
1503 self
._rational
_algebra
= FiniteDimensionalEJA(
1508 orthonormalize
=False,
1512 super().__init
__(basis
,
1516 check_field
=check_field
,
1520 def _charpoly_coefficients(self
):
1524 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
1525 ....: JordanSpinEJA)
1529 The base ring of the resulting polynomial coefficients is what
1530 it should be, and not the rationals (unless the algebra was
1531 already over the rationals)::
1533 sage: J = JordanSpinEJA(3)
1534 sage: J._charpoly_coefficients()
1535 (X1^2 - X2^2 - X3^2, -2*X1)
1536 sage: a0 = J._charpoly_coefficients()[0]
1538 Algebraic Real Field
1539 sage: a0.base_ring()
1540 Algebraic Real Field
1543 if self
._rational
_algebra
is None:
1544 # There's no need to construct *another* algebra over the
1545 # rationals if this one is already over the
1546 # rationals. Likewise, if we never orthonormalized our
1547 # basis, we might as well just use the given one.
1548 return super()._charpoly
_coefficients
()
1550 # Do the computation over the rationals. The answer will be
1551 # the same, because all we've done is a change of basis.
1552 # Then, change back from QQ to our real base ring
1553 a
= ( a_i
.change_ring(self
.base_ring())
1554 for a_i
in self
._rational
_algebra
._charpoly
_coefficients
() )
1556 if self
._deortho
_matrix
is None:
1557 # This can happen if our base ring was, say, AA and we
1558 # chose not to (or didn't need to) orthonormalize. It's
1559 # still faster to do the computations over QQ even if
1560 # the numbers in the boxes stay the same.
1563 # Otherwise, convert the coordinate variables back to the
1564 # deorthonormalized ones.
1565 R
= self
.coordinate_polynomial_ring()
1566 from sage
.modules
.free_module_element
import vector
1567 X
= vector(R
, R
.gens())
1568 BX
= self
._deortho
_matrix
*X
1570 subs_dict
= { X[i]: BX[i] for i in range(len(X)) }
1571 return tuple( a_i
.subs(subs_dict
) for a_i
in a
)
1573 class ConcreteEJA(RationalBasisEJA
):
1575 A class for the Euclidean Jordan algebras that we know by name.
1577 These are the Jordan algebras whose basis, multiplication table,
1578 rank, and so on are known a priori. More to the point, they are
1579 the Euclidean Jordan algebras for which we are able to conjure up
1580 a "random instance."
1584 sage: from mjo.eja.eja_algebra import ConcreteEJA
1588 Our basis is normalized with respect to the algebra's inner
1589 product, unless we specify otherwise::
1591 sage: set_random_seed()
1592 sage: J = ConcreteEJA.random_instance()
1593 sage: all( b.norm() == 1 for b in J.gens() )
1596 Since our basis is orthonormal with respect to the algebra's inner
1597 product, and since we know that this algebra is an EJA, any
1598 left-multiplication operator's matrix will be symmetric because
1599 natural->EJA basis representation is an isometry and within the
1600 EJA the operator is self-adjoint by the Jordan axiom::
1602 sage: set_random_seed()
1603 sage: J = ConcreteEJA.random_instance()
1604 sage: x = J.random_element()
1605 sage: x.operator().is_self_adjoint()
1610 def _max_random_instance_size():
1612 Return an integer "size" that is an upper bound on the size of
1613 this algebra when it is used in a random test
1614 case. Unfortunately, the term "size" is ambiguous -- when
1615 dealing with `R^n` under either the Hadamard or Jordan spin
1616 product, the "size" refers to the dimension `n`. When dealing
1617 with a matrix algebra (real symmetric or complex/quaternion
1618 Hermitian), it refers to the size of the matrix, which is far
1619 less than the dimension of the underlying vector space.
1621 This method must be implemented in each subclass.
1623 raise NotImplementedError
1626 def random_instance(cls
, *args
, **kwargs
):
1628 Return a random instance of this type of algebra.
1630 This method should be implemented in each subclass.
1632 from sage
.misc
.prandom
import choice
1633 eja_class
= choice(cls
.__subclasses
__())
1635 # These all bubble up to the RationalBasisEJA superclass
1636 # constructor, so any (kw)args valid there are also valid
1638 return eja_class
.random_instance(*args
, **kwargs
)
1643 def dimension_over_reals():
1645 The dimension of this matrix's base ring over the reals.
1647 The reals are dimension one over themselves, obviously; that's
1648 just `\mathbb{R}^{1}`. Likewise, the complex numbers `a + bi`
1649 have dimension two. Finally, the quaternions have dimension
1650 four over the reals.
1652 This is used to determine the size of the matrix returned from
1653 :meth:`real_embed`, among other things.
1655 raise NotImplementedError
1658 def real_embed(cls
,M
):
1660 Embed the matrix ``M`` into a space of real matrices.
1662 The matrix ``M`` can have entries in any field at the moment:
1663 the real numbers, complex numbers, or quaternions. And although
1664 they are not a field, we can probably support octonions at some
1665 point, too. This function returns a real matrix that "acts like"
1666 the original with respect to matrix multiplication; i.e.
1668 real_embed(M*N) = real_embed(M)*real_embed(N)
1671 if M
.ncols() != M
.nrows():
1672 raise ValueError("the matrix 'M' must be square")
1677 def real_unembed(cls
,M
):
1679 The inverse of :meth:`real_embed`.
1681 if M
.ncols() != M
.nrows():
1682 raise ValueError("the matrix 'M' must be square")
1683 if not ZZ(M
.nrows()).mod(cls
.dimension_over_reals()).is_zero():
1684 raise ValueError("the matrix 'M' must be a real embedding")
1688 def jordan_product(X
,Y
):
1689 return (X
*Y
+ Y
*X
)/2
1692 def trace_inner_product(cls
,X
,Y
):
1694 Compute the trace inner-product of two real-embeddings.
1698 sage: from mjo.eja.eja_algebra import (RealSymmetricEJA,
1699 ....: ComplexHermitianEJA,
1700 ....: QuaternionHermitianEJA)
1704 This gives the same answer as it would if we computed the trace
1705 from the unembedded (original) matrices::
1707 sage: set_random_seed()
1708 sage: J = RealSymmetricEJA.random_instance()
1709 sage: x,y = J.random_elements(2)
1710 sage: Xe = x.to_matrix()
1711 sage: Ye = y.to_matrix()
1712 sage: X = J.real_unembed(Xe)
1713 sage: Y = J.real_unembed(Ye)
1714 sage: expected = (X*Y).trace()
1715 sage: actual = J.trace_inner_product(Xe,Ye)
1716 sage: actual == expected
1721 sage: set_random_seed()
1722 sage: J = ComplexHermitianEJA.random_instance()
1723 sage: x,y = J.random_elements(2)
1724 sage: Xe = x.to_matrix()
1725 sage: Ye = y.to_matrix()
1726 sage: X = J.real_unembed(Xe)
1727 sage: Y = J.real_unembed(Ye)
1728 sage: expected = (X*Y).trace().real()
1729 sage: actual = J.trace_inner_product(Xe,Ye)
1730 sage: actual == expected
1735 sage: set_random_seed()
1736 sage: J = QuaternionHermitianEJA.random_instance()
1737 sage: x,y = J.random_elements(2)
1738 sage: Xe = x.to_matrix()
1739 sage: Ye = y.to_matrix()
1740 sage: X = J.real_unembed(Xe)
1741 sage: Y = J.real_unembed(Ye)
1742 sage: expected = (X*Y).trace().coefficient_tuple()[0]
1743 sage: actual = J.trace_inner_product(Xe,Ye)
1744 sage: actual == expected
1748 Xu
= cls
.real_unembed(X
)
1749 Yu
= cls
.real_unembed(Y
)
1750 tr
= (Xu
*Yu
).trace()
1753 # Works in QQ, AA, RDF, et cetera.
1755 except AttributeError:
1756 # A quaternion doesn't have a real() method, but does
1757 # have coefficient_tuple() method that returns the
1758 # coefficients of 1, i, j, and k -- in that order.
1759 return tr
.coefficient_tuple()[0]
1762 class RealMatrixEJA(MatrixEJA
):
1764 def dimension_over_reals():
1768 class RealSymmetricEJA(ConcreteEJA
, RealMatrixEJA
):
1770 The rank-n simple EJA consisting of real symmetric n-by-n
1771 matrices, the usual symmetric Jordan product, and the trace inner
1772 product. It has dimension `(n^2 + n)/2` over the reals.
1776 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1780 sage: J = RealSymmetricEJA(2)
1781 sage: e0, e1, e2 = J.gens()
1789 In theory, our "field" can be any subfield of the reals::
1791 sage: RealSymmetricEJA(2, field=RDF, check_axioms=True)
1792 Euclidean Jordan algebra of dimension 3 over Real Double Field
1793 sage: RealSymmetricEJA(2, field=RR, check_axioms=True)
1794 Euclidean Jordan algebra of dimension 3 over Real Field with
1795 53 bits of precision
1799 The dimension of this algebra is `(n^2 + n) / 2`::
1801 sage: set_random_seed()
1802 sage: n_max = RealSymmetricEJA._max_random_instance_size()
1803 sage: n = ZZ.random_element(1, n_max)
1804 sage: J = RealSymmetricEJA(n)
1805 sage: J.dimension() == (n^2 + n)/2
1808 The Jordan multiplication is what we think it is::
1810 sage: set_random_seed()
1811 sage: J = RealSymmetricEJA.random_instance()
1812 sage: x,y = J.random_elements(2)
1813 sage: actual = (x*y).to_matrix()
1814 sage: X = x.to_matrix()
1815 sage: Y = y.to_matrix()
1816 sage: expected = (X*Y + Y*X)/2
1817 sage: actual == expected
1819 sage: J(expected) == x*y
1822 We can change the generator prefix::
1824 sage: RealSymmetricEJA(3, prefix='q').gens()
1825 (q0, q1, q2, q3, q4, q5)
1827 We can construct the (trivial) algebra of rank zero::
1829 sage: RealSymmetricEJA(0)
1830 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1834 def _denormalized_basis(cls
, n
):
1836 Return a basis for the space of real symmetric n-by-n matrices.
1840 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1844 sage: set_random_seed()
1845 sage: n = ZZ.random_element(1,5)
1846 sage: B = RealSymmetricEJA._denormalized_basis(n)
1847 sage: all( M.is_symmetric() for M in B)
1851 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1855 for j
in range(i
+1):
1856 Eij
= matrix(ZZ
, n
, lambda k
,l
: k
==i
and l
==j
)
1860 Sij
= Eij
+ Eij
.transpose()
1866 def _max_random_instance_size():
1867 return 4 # Dimension 10
1870 def random_instance(cls
, **kwargs
):
1872 Return a random instance of this type of algebra.
1874 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1875 return cls(n
, **kwargs
)
1877 def __init__(self
, n
, **kwargs
):
1878 # We know this is a valid EJA, but will double-check
1879 # if the user passes check_axioms=True.
1880 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
1882 super(RealSymmetricEJA
, self
).__init
__(self
._denormalized
_basis
(n
),
1883 self
.jordan_product
,
1884 self
.trace_inner_product
,
1887 # TODO: this could be factored out somehow, but is left here
1888 # because the MatrixEJA is not presently a subclass of the
1889 # FDEJA class that defines rank() and one().
1890 self
.rank
.set_cache(n
)
1891 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
1892 self
.one
.set_cache(self(idV
))
1896 class ComplexMatrixEJA(MatrixEJA
):
1897 # A manual dictionary-cache for the complex_extension() method,
1898 # since apparently @classmethods can't also be @cached_methods.
1899 _complex_extension
= {}
1902 def complex_extension(cls
,field
):
1904 The complex field that we embed/unembed, as an extension
1905 of the given ``field``.
1907 if field
in cls
._complex
_extension
:
1908 return cls
._complex
_extension
[field
]
1910 # Sage doesn't know how to adjoin the complex "i" (the root of
1911 # x^2 + 1) to a field in a general way. Here, we just enumerate
1912 # all of the cases that I have cared to support so far.
1914 # Sage doesn't know how to embed AA into QQbar, i.e. how
1915 # to adjoin sqrt(-1) to AA.
1917 elif not field
.is_exact():
1919 F
= field
.complex_field()
1921 # Works for QQ and... maybe some other fields.
1922 R
= PolynomialRing(field
, 'z')
1924 F
= field
.extension(z
**2 + 1, 'I', embedding
=CLF(-1).sqrt())
1926 cls
._complex
_extension
[field
] = F
1930 def dimension_over_reals():
1934 def real_embed(cls
,M
):
1936 Embed the n-by-n complex matrix ``M`` into the space of real
1937 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1938 bi` to the block matrix ``[[a,b],[-b,a]]``.
1942 sage: from mjo.eja.eja_algebra import ComplexMatrixEJA
1946 sage: F = QuadraticField(-1, 'I')
1947 sage: x1 = F(4 - 2*i)
1948 sage: x2 = F(1 + 2*i)
1951 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1952 sage: ComplexMatrixEJA.real_embed(M)
1961 Embedding is a homomorphism (isomorphism, in fact)::
1963 sage: set_random_seed()
1964 sage: n = ZZ.random_element(3)
1965 sage: F = QuadraticField(-1, 'I')
1966 sage: X = random_matrix(F, n)
1967 sage: Y = random_matrix(F, n)
1968 sage: Xe = ComplexMatrixEJA.real_embed(X)
1969 sage: Ye = ComplexMatrixEJA.real_embed(Y)
1970 sage: XYe = ComplexMatrixEJA.real_embed(X*Y)
1975 super(ComplexMatrixEJA
,cls
).real_embed(M
)
1978 # We don't need any adjoined elements...
1979 field
= M
.base_ring().base_ring()
1985 blocks
.append(matrix(field
, 2, [ [ a
, b
],
1988 return matrix
.block(field
, n
, blocks
)
1992 def real_unembed(cls
,M
):
1994 The inverse of _embed_complex_matrix().
1998 sage: from mjo.eja.eja_algebra import ComplexMatrixEJA
2002 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
2003 ....: [-2, 1, -4, 3],
2004 ....: [ 9, 10, 11, 12],
2005 ....: [-10, 9, -12, 11] ])
2006 sage: ComplexMatrixEJA.real_unembed(A)
2008 [ 10*I + 9 12*I + 11]
2012 Unembedding is the inverse of embedding::
2014 sage: set_random_seed()
2015 sage: F = QuadraticField(-1, 'I')
2016 sage: M = random_matrix(F, 3)
2017 sage: Me = ComplexMatrixEJA.real_embed(M)
2018 sage: ComplexMatrixEJA.real_unembed(Me) == M
2022 super(ComplexMatrixEJA
,cls
).real_unembed(M
)
2024 d
= cls
.dimension_over_reals()
2025 F
= cls
.complex_extension(M
.base_ring())
2028 # Go top-left to bottom-right (reading order), converting every
2029 # 2-by-2 block we see to a single complex element.
2031 for k
in range(n
/d
):
2032 for j
in range(n
/d
):
2033 submat
= M
[d
*k
:d
*k
+d
,d
*j
:d
*j
+d
]
2034 if submat
[0,0] != submat
[1,1]:
2035 raise ValueError('bad on-diagonal submatrix')
2036 if submat
[0,1] != -submat
[1,0]:
2037 raise ValueError('bad off-diagonal submatrix')
2038 z
= submat
[0,0] + submat
[0,1]*i
2041 return matrix(F
, n
/d
, elements
)
2044 class ComplexHermitianEJA(ConcreteEJA
, ComplexMatrixEJA
):
2046 The rank-n simple EJA consisting of complex Hermitian n-by-n
2047 matrices over the real numbers, the usual symmetric Jordan product,
2048 and the real-part-of-trace inner product. It has dimension `n^2` over
2053 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
2057 In theory, our "field" can be any subfield of the reals::
2059 sage: ComplexHermitianEJA(2, field=RDF, check_axioms=True)
2060 Euclidean Jordan algebra of dimension 4 over Real Double Field
2061 sage: ComplexHermitianEJA(2, field=RR, check_axioms=True)
2062 Euclidean Jordan algebra of dimension 4 over Real Field with
2063 53 bits of precision
2067 The dimension of this algebra is `n^2`::
2069 sage: set_random_seed()
2070 sage: n_max = ComplexHermitianEJA._max_random_instance_size()
2071 sage: n = ZZ.random_element(1, n_max)
2072 sage: J = ComplexHermitianEJA(n)
2073 sage: J.dimension() == n^2
2076 The Jordan multiplication is what we think it is::
2078 sage: set_random_seed()
2079 sage: J = ComplexHermitianEJA.random_instance()
2080 sage: x,y = J.random_elements(2)
2081 sage: actual = (x*y).to_matrix()
2082 sage: X = x.to_matrix()
2083 sage: Y = y.to_matrix()
2084 sage: expected = (X*Y + Y*X)/2
2085 sage: actual == expected
2087 sage: J(expected) == x*y
2090 We can change the generator prefix::
2092 sage: ComplexHermitianEJA(2, prefix='z').gens()
2095 We can construct the (trivial) algebra of rank zero::
2097 sage: ComplexHermitianEJA(0)
2098 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2103 def _denormalized_basis(cls
, n
):
2105 Returns a basis for the space of complex Hermitian n-by-n matrices.
2107 Why do we embed these? Basically, because all of numerical linear
2108 algebra assumes that you're working with vectors consisting of `n`
2109 entries from a field and scalars from the same field. There's no way
2110 to tell SageMath that (for example) the vectors contain complex
2111 numbers, while the scalar field is real.
2115 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
2119 sage: set_random_seed()
2120 sage: n = ZZ.random_element(1,5)
2121 sage: B = ComplexHermitianEJA._denormalized_basis(n)
2122 sage: all( M.is_symmetric() for M in B)
2127 R
= PolynomialRing(field
, 'z')
2129 F
= field
.extension(z
**2 + 1, 'I')
2132 # This is like the symmetric case, but we need to be careful:
2134 # * We want conjugate-symmetry, not just symmetry.
2135 # * The diagonal will (as a result) be real.
2138 Eij
= matrix
.zero(F
,n
)
2140 for j
in range(i
+1):
2144 Sij
= cls
.real_embed(Eij
)
2147 # The second one has a minus because it's conjugated.
2148 Eij
[j
,i
] = 1 # Eij = Eij + Eij.transpose()
2149 Sij_real
= cls
.real_embed(Eij
)
2151 # Eij = I*Eij - I*Eij.transpose()
2154 Sij_imag
= cls
.real_embed(Eij
)
2160 # Since we embedded these, we can drop back to the "field" that we
2161 # started with instead of the complex extension "F".
2162 return tuple( s
.change_ring(field
) for s
in S
)
2165 def __init__(self
, n
, **kwargs
):
2166 # We know this is a valid EJA, but will double-check
2167 # if the user passes check_axioms=True.
2168 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2170 super(ComplexHermitianEJA
, self
).__init
__(self
._denormalized
_basis
(n
),
2171 self
.jordan_product
,
2172 self
.trace_inner_product
,
2174 # TODO: this could be factored out somehow, but is left here
2175 # because the MatrixEJA is not presently a subclass of the
2176 # FDEJA class that defines rank() and one().
2177 self
.rank
.set_cache(n
)
2178 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
2179 self
.one
.set_cache(self(idV
))
2182 def _max_random_instance_size():
2183 return 3 # Dimension 9
2186 def random_instance(cls
, **kwargs
):
2188 Return a random instance of this type of algebra.
2190 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2191 return cls(n
, **kwargs
)
2193 class QuaternionMatrixEJA(MatrixEJA
):
2195 # A manual dictionary-cache for the quaternion_extension() method,
2196 # since apparently @classmethods can't also be @cached_methods.
2197 _quaternion_extension
= {}
2200 def quaternion_extension(cls
,field
):
2202 The quaternion field that we embed/unembed, as an extension
2203 of the given ``field``.
2205 if field
in cls
._quaternion
_extension
:
2206 return cls
._quaternion
_extension
[field
]
2208 Q
= QuaternionAlgebra(field
,-1,-1)
2210 cls
._quaternion
_extension
[field
] = Q
2214 def dimension_over_reals():
2218 def real_embed(cls
,M
):
2220 Embed the n-by-n quaternion matrix ``M`` into the space of real
2221 matrices of size 4n-by-4n by first sending each quaternion entry `z
2222 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
2223 c+di],[-c + di, a-bi]]`, and then embedding those into a real
2228 sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA
2232 sage: Q = QuaternionAlgebra(QQ,-1,-1)
2233 sage: i,j,k = Q.gens()
2234 sage: x = 1 + 2*i + 3*j + 4*k
2235 sage: M = matrix(Q, 1, [[x]])
2236 sage: QuaternionMatrixEJA.real_embed(M)
2242 Embedding is a homomorphism (isomorphism, in fact)::
2244 sage: set_random_seed()
2245 sage: n = ZZ.random_element(2)
2246 sage: Q = QuaternionAlgebra(QQ,-1,-1)
2247 sage: X = random_matrix(Q, n)
2248 sage: Y = random_matrix(Q, n)
2249 sage: Xe = QuaternionMatrixEJA.real_embed(X)
2250 sage: Ye = QuaternionMatrixEJA.real_embed(Y)
2251 sage: XYe = QuaternionMatrixEJA.real_embed(X*Y)
2256 super(QuaternionMatrixEJA
,cls
).real_embed(M
)
2257 quaternions
= M
.base_ring()
2260 F
= QuadraticField(-1, 'I')
2265 t
= z
.coefficient_tuple()
2270 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
2271 [-c
+ d
*i
, a
- b
*i
]])
2272 realM
= ComplexMatrixEJA
.real_embed(cplxM
)
2273 blocks
.append(realM
)
2275 # We should have real entries by now, so use the realest field
2276 # we've got for the return value.
2277 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
2282 def real_unembed(cls
,M
):
2284 The inverse of _embed_quaternion_matrix().
2288 sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA
2292 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
2293 ....: [-2, 1, -4, 3],
2294 ....: [-3, 4, 1, -2],
2295 ....: [-4, -3, 2, 1]])
2296 sage: QuaternionMatrixEJA.real_unembed(M)
2297 [1 + 2*i + 3*j + 4*k]
2301 Unembedding is the inverse of embedding::
2303 sage: set_random_seed()
2304 sage: Q = QuaternionAlgebra(QQ, -1, -1)
2305 sage: M = random_matrix(Q, 3)
2306 sage: Me = QuaternionMatrixEJA.real_embed(M)
2307 sage: QuaternionMatrixEJA.real_unembed(Me) == M
2311 super(QuaternionMatrixEJA
,cls
).real_unembed(M
)
2313 d
= cls
.dimension_over_reals()
2315 # Use the base ring of the matrix to ensure that its entries can be
2316 # multiplied by elements of the quaternion algebra.
2317 Q
= cls
.quaternion_extension(M
.base_ring())
2320 # Go top-left to bottom-right (reading order), converting every
2321 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
2324 for l
in range(n
/d
):
2325 for m
in range(n
/d
):
2326 submat
= ComplexMatrixEJA
.real_unembed(
2327 M
[d
*l
:d
*l
+d
,d
*m
:d
*m
+d
] )
2328 if submat
[0,0] != submat
[1,1].conjugate():
2329 raise ValueError('bad on-diagonal submatrix')
2330 if submat
[0,1] != -submat
[1,0].conjugate():
2331 raise ValueError('bad off-diagonal submatrix')
2332 z
= submat
[0,0].real()
2333 z
+= submat
[0,0].imag()*i
2334 z
+= submat
[0,1].real()*j
2335 z
+= submat
[0,1].imag()*k
2338 return matrix(Q
, n
/d
, elements
)
2341 class QuaternionHermitianEJA(ConcreteEJA
, QuaternionMatrixEJA
):
2343 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
2344 matrices, the usual symmetric Jordan product, and the
2345 real-part-of-trace inner product. It has dimension `2n^2 - n` over
2350 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2354 In theory, our "field" can be any subfield of the reals::
2356 sage: QuaternionHermitianEJA(2, field=RDF, check_axioms=True)
2357 Euclidean Jordan algebra of dimension 6 over Real Double Field
2358 sage: QuaternionHermitianEJA(2, field=RR, check_axioms=True)
2359 Euclidean Jordan algebra of dimension 6 over Real Field with
2360 53 bits of precision
2364 The dimension of this algebra is `2*n^2 - n`::
2366 sage: set_random_seed()
2367 sage: n_max = QuaternionHermitianEJA._max_random_instance_size()
2368 sage: n = ZZ.random_element(1, n_max)
2369 sage: J = QuaternionHermitianEJA(n)
2370 sage: J.dimension() == 2*(n^2) - n
2373 The Jordan multiplication is what we think it is::
2375 sage: set_random_seed()
2376 sage: J = QuaternionHermitianEJA.random_instance()
2377 sage: x,y = J.random_elements(2)
2378 sage: actual = (x*y).to_matrix()
2379 sage: X = x.to_matrix()
2380 sage: Y = y.to_matrix()
2381 sage: expected = (X*Y + Y*X)/2
2382 sage: actual == expected
2384 sage: J(expected) == x*y
2387 We can change the generator prefix::
2389 sage: QuaternionHermitianEJA(2, prefix='a').gens()
2390 (a0, a1, a2, a3, a4, a5)
2392 We can construct the (trivial) algebra of rank zero::
2394 sage: QuaternionHermitianEJA(0)
2395 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2399 def _denormalized_basis(cls
, n
):
2401 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
2403 Why do we embed these? Basically, because all of numerical
2404 linear algebra assumes that you're working with vectors consisting
2405 of `n` entries from a field and scalars from the same field. There's
2406 no way to tell SageMath that (for example) the vectors contain
2407 complex numbers, while the scalar field is real.
2411 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2415 sage: set_random_seed()
2416 sage: n = ZZ.random_element(1,5)
2417 sage: B = QuaternionHermitianEJA._denormalized_basis(n)
2418 sage: all( M.is_symmetric() for M in B )
2423 Q
= QuaternionAlgebra(QQ
,-1,-1)
2426 # This is like the symmetric case, but we need to be careful:
2428 # * We want conjugate-symmetry, not just symmetry.
2429 # * The diagonal will (as a result) be real.
2432 Eij
= matrix
.zero(Q
,n
)
2434 for j
in range(i
+1):
2438 Sij
= cls
.real_embed(Eij
)
2441 # The second, third, and fourth ones have a minus
2442 # because they're conjugated.
2443 # Eij = Eij + Eij.transpose()
2445 Sij_real
= cls
.real_embed(Eij
)
2447 # Eij = I*(Eij - Eij.transpose())
2450 Sij_I
= cls
.real_embed(Eij
)
2452 # Eij = J*(Eij - Eij.transpose())
2455 Sij_J
= cls
.real_embed(Eij
)
2457 # Eij = K*(Eij - Eij.transpose())
2460 Sij_K
= cls
.real_embed(Eij
)
2466 # Since we embedded these, we can drop back to the "field" that we
2467 # started with instead of the quaternion algebra "Q".
2468 return tuple( s
.change_ring(field
) for s
in S
)
2471 def __init__(self
, n
, **kwargs
):
2472 # We know this is a valid EJA, but will double-check
2473 # if the user passes check_axioms=True.
2474 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2476 super(QuaternionHermitianEJA
, self
).__init
__(self
._denormalized
_basis
(n
),
2477 self
.jordan_product
,
2478 self
.trace_inner_product
,
2480 # TODO: this could be factored out somehow, but is left here
2481 # because the MatrixEJA is not presently a subclass of the
2482 # FDEJA class that defines rank() and one().
2483 self
.rank
.set_cache(n
)
2484 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
2485 self
.one
.set_cache(self(idV
))
2489 def _max_random_instance_size():
2491 The maximum rank of a random QuaternionHermitianEJA.
2493 return 2 # Dimension 6
2496 def random_instance(cls
, **kwargs
):
2498 Return a random instance of this type of algebra.
2500 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2501 return cls(n
, **kwargs
)
2504 class HadamardEJA(ConcreteEJA
):
2506 Return the Euclidean Jordan Algebra corresponding to the set
2507 `R^n` under the Hadamard product.
2509 Note: this is nothing more than the Cartesian product of ``n``
2510 copies of the spin algebra. Once Cartesian product algebras
2511 are implemented, this can go.
2515 sage: from mjo.eja.eja_algebra import HadamardEJA
2519 This multiplication table can be verified by hand::
2521 sage: J = HadamardEJA(3)
2522 sage: e0,e1,e2 = J.gens()
2538 We can change the generator prefix::
2540 sage: HadamardEJA(3, prefix='r').gens()
2544 def __init__(self
, n
, **kwargs
):
2546 jordan_product
= lambda x
,y
: x
2547 inner_product
= lambda x
,y
: x
2549 def jordan_product(x
,y
):
2551 return P( xi
*yi
for (xi
,yi
) in zip(x
,y
) )
2553 def inner_product(x
,y
):
2556 # New defaults for keyword arguments. Don't orthonormalize
2557 # because our basis is already orthonormal with respect to our
2558 # inner-product. Don't check the axioms, because we know this
2559 # is a valid EJA... but do double-check if the user passes
2560 # check_axioms=True. Note: we DON'T override the "check_field"
2561 # default here, because the user can pass in a field!
2562 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2563 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2565 column_basis
= tuple( b
.column() for b
in FreeModule(ZZ
, n
).basis() )
2566 super().__init
__(column_basis
,
2571 self
.rank
.set_cache(n
)
2574 self
.one
.set_cache( self
.zero() )
2576 self
.one
.set_cache( sum(self
.gens()) )
2579 def _max_random_instance_size():
2581 The maximum dimension of a random HadamardEJA.
2586 def random_instance(cls
, **kwargs
):
2588 Return a random instance of this type of algebra.
2590 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2591 return cls(n
, **kwargs
)
2594 class BilinearFormEJA(ConcreteEJA
):
2596 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2597 with the half-trace inner product and jordan product ``x*y =
2598 (<Bx,y>,y_bar>, x0*y_bar + y0*x_bar)`` where `B = 1 \times B22` is
2599 a symmetric positive-definite "bilinear form" matrix. Its
2600 dimension is the size of `B`, and it has rank two in dimensions
2601 larger than two. It reduces to the ``JordanSpinEJA`` when `B` is
2602 the identity matrix of order ``n``.
2604 We insist that the one-by-one upper-left identity block of `B` be
2605 passed in as well so that we can be passed a matrix of size zero
2606 to construct a trivial algebra.
2610 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
2611 ....: JordanSpinEJA)
2615 When no bilinear form is specified, the identity matrix is used,
2616 and the resulting algebra is the Jordan spin algebra::
2618 sage: B = matrix.identity(AA,3)
2619 sage: J0 = BilinearFormEJA(B)
2620 sage: J1 = JordanSpinEJA(3)
2621 sage: J0.multiplication_table() == J0.multiplication_table()
2624 An error is raised if the matrix `B` does not correspond to a
2625 positive-definite bilinear form::
2627 sage: B = matrix.random(QQ,2,3)
2628 sage: J = BilinearFormEJA(B)
2629 Traceback (most recent call last):
2631 ValueError: bilinear form is not positive-definite
2632 sage: B = matrix.zero(QQ,3)
2633 sage: J = BilinearFormEJA(B)
2634 Traceback (most recent call last):
2636 ValueError: bilinear form is not positive-definite
2640 We can create a zero-dimensional algebra::
2642 sage: B = matrix.identity(AA,0)
2643 sage: J = BilinearFormEJA(B)
2647 We can check the multiplication condition given in the Jordan, von
2648 Neumann, and Wigner paper (and also discussed on my "On the
2649 symmetry..." paper). Note that this relies heavily on the standard
2650 choice of basis, as does anything utilizing the bilinear form
2651 matrix. We opt not to orthonormalize the basis, because if we
2652 did, we would have to normalize the `s_{i}` in a similar manner::
2654 sage: set_random_seed()
2655 sage: n = ZZ.random_element(5)
2656 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2657 sage: B11 = matrix.identity(QQ,1)
2658 sage: B22 = M.transpose()*M
2659 sage: B = block_matrix(2,2,[ [B11,0 ],
2661 sage: J = BilinearFormEJA(B, orthonormalize=False)
2662 sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
2663 sage: V = J.vector_space()
2664 sage: sis = [ J( V([0] + (M.inverse()*ei).list()).column() )
2665 ....: for ei in eis ]
2666 sage: actual = [ sis[i]*sis[j]
2667 ....: for i in range(n-1)
2668 ....: for j in range(n-1) ]
2669 sage: expected = [ J.one() if i == j else J.zero()
2670 ....: for i in range(n-1)
2671 ....: for j in range(n-1) ]
2672 sage: actual == expected
2676 def __init__(self
, B
, **kwargs
):
2677 # The matrix "B" is supplied by the user in most cases,
2678 # so it makes sense to check whether or not its positive-
2679 # definite unless we are specifically asked not to...
2680 if ("check_axioms" not in kwargs
) or kwargs
["check_axioms"]:
2681 if not B
.is_positive_definite():
2682 raise ValueError("bilinear form is not positive-definite")
2684 # However, all of the other data for this EJA is computed
2685 # by us in manner that guarantees the axioms are
2686 # satisfied. So, again, unless we are specifically asked to
2687 # verify things, we'll skip the rest of the checks.
2688 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2690 def inner_product(x
,y
):
2691 return (y
.T
*B
*x
)[0,0]
2693 def jordan_product(x
,y
):
2699 z0
= inner_product(y
,x
)
2700 zbar
= y0
*xbar
+ x0
*ybar
2701 return P([z0
] + zbar
.list())
2704 column_basis
= tuple( b
.column() for b
in FreeModule(ZZ
, n
).basis() )
2705 super(BilinearFormEJA
, self
).__init
__(column_basis
,
2710 # The rank of this algebra is two, unless we're in a
2711 # one-dimensional ambient space (because the rank is bounded
2712 # by the ambient dimension).
2713 self
.rank
.set_cache(min(n
,2))
2716 self
.one
.set_cache( self
.zero() )
2718 self
.one
.set_cache( self
.monomial(0) )
2721 def _max_random_instance_size():
2723 The maximum dimension of a random BilinearFormEJA.
2728 def random_instance(cls
, **kwargs
):
2730 Return a random instance of this algebra.
2732 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2734 B
= matrix
.identity(ZZ
, n
)
2735 return cls(B
, **kwargs
)
2737 B11
= matrix
.identity(ZZ
, 1)
2738 M
= matrix
.random(ZZ
, n
-1)
2739 I
= matrix
.identity(ZZ
, n
-1)
2741 while alpha
.is_zero():
2742 alpha
= ZZ
.random_element().abs()
2743 B22
= M
.transpose()*M
+ alpha
*I
2745 from sage
.matrix
.special
import block_matrix
2746 B
= block_matrix(2,2, [ [B11
, ZZ(0) ],
2749 return cls(B
, **kwargs
)
2752 class JordanSpinEJA(BilinearFormEJA
):
2754 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2755 with the usual inner product and jordan product ``x*y =
2756 (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2761 sage: from mjo.eja.eja_algebra import JordanSpinEJA
2765 This multiplication table can be verified by hand::
2767 sage: J = JordanSpinEJA(4)
2768 sage: e0,e1,e2,e3 = J.gens()
2784 We can change the generator prefix::
2786 sage: JordanSpinEJA(2, prefix='B').gens()
2791 Ensure that we have the usual inner product on `R^n`::
2793 sage: set_random_seed()
2794 sage: J = JordanSpinEJA.random_instance()
2795 sage: x,y = J.random_elements(2)
2796 sage: actual = x.inner_product(y)
2797 sage: expected = x.to_vector().inner_product(y.to_vector())
2798 sage: actual == expected
2802 def __init__(self
, n
, **kwargs
):
2803 # This is a special case of the BilinearFormEJA with the
2804 # identity matrix as its bilinear form.
2805 B
= matrix
.identity(ZZ
, n
)
2807 # Don't orthonormalize because our basis is already
2808 # orthonormal with respect to our inner-product.
2809 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2811 # But also don't pass check_field=False here, because the user
2812 # can pass in a field!
2813 super(JordanSpinEJA
, self
).__init
__(B
, **kwargs
)
2816 def _max_random_instance_size():
2818 The maximum dimension of a random JordanSpinEJA.
2823 def random_instance(cls
, **kwargs
):
2825 Return a random instance of this type of algebra.
2827 Needed here to override the implementation for ``BilinearFormEJA``.
2829 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2830 return cls(n
, **kwargs
)
2833 class TrivialEJA(ConcreteEJA
):
2835 The trivial Euclidean Jordan algebra consisting of only a zero element.
2839 sage: from mjo.eja.eja_algebra import TrivialEJA
2843 sage: J = TrivialEJA()
2850 sage: 7*J.one()*12*J.one()
2852 sage: J.one().inner_product(J.one())
2854 sage: J.one().norm()
2856 sage: J.one().subalgebra_generated_by()
2857 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2862 def __init__(self
, **kwargs
):
2863 jordan_product
= lambda x
,y
: x
2864 inner_product
= lambda x
,y
: 0
2867 # New defaults for keyword arguments
2868 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2869 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2871 super(TrivialEJA
, self
).__init
__(basis
,
2875 # The rank is zero using my definition, namely the dimension of the
2876 # largest subalgebra generated by any element.
2877 self
.rank
.set_cache(0)
2878 self
.one
.set_cache( self
.zero() )
2881 def random_instance(cls
, **kwargs
):
2882 # We don't take a "size" argument so the superclass method is
2883 # inappropriate for us.
2884 return cls(**kwargs
)
2887 class CartesianProductEJA(CombinatorialFreeModule_CartesianProduct
,
2888 FiniteDimensionalEJA
):
2890 The external (orthogonal) direct sum of two or more Euclidean
2891 Jordan algebras. Every Euclidean Jordan algebra decomposes into an
2892 orthogonal direct sum of simple Euclidean Jordan algebras which is
2893 then isometric to a Cartesian product, so no generality is lost by
2894 providing only this construction.
2898 sage: from mjo.eja.eja_algebra import (random_eja,
2899 ....: CartesianProductEJA,
2901 ....: JordanSpinEJA,
2902 ....: RealSymmetricEJA)
2906 The Jordan product is inherited from our factors and implemented by
2907 our CombinatorialFreeModule Cartesian product superclass::
2909 sage: set_random_seed()
2910 sage: J1 = HadamardEJA(2)
2911 sage: J2 = RealSymmetricEJA(2)
2912 sage: J = cartesian_product([J1,J2])
2913 sage: x,y = J.random_elements(2)
2917 The ability to retrieve the original factors is implemented by our
2918 CombinatorialFreeModule Cartesian product superclass::
2920 sage: J1 = HadamardEJA(2, field=QQ)
2921 sage: J2 = JordanSpinEJA(3, field=QQ)
2922 sage: J = cartesian_product([J1,J2])
2923 sage: J.cartesian_factors()
2924 (Euclidean Jordan algebra of dimension 2 over Rational Field,
2925 Euclidean Jordan algebra of dimension 3 over Rational Field)
2927 You can provide more than two factors::
2929 sage: J1 = HadamardEJA(2)
2930 sage: J2 = JordanSpinEJA(3)
2931 sage: J3 = RealSymmetricEJA(3)
2932 sage: cartesian_product([J1,J2,J3])
2933 Euclidean Jordan algebra of dimension 2 over Algebraic Real
2934 Field (+) Euclidean Jordan algebra of dimension 3 over Algebraic
2935 Real Field (+) Euclidean Jordan algebra of dimension 6 over
2936 Algebraic Real Field
2938 Rank is additive on a Cartesian product::
2940 sage: J1 = HadamardEJA(1)
2941 sage: J2 = RealSymmetricEJA(2)
2942 sage: J = cartesian_product([J1,J2])
2943 sage: J1.rank.clear_cache()
2944 sage: J2.rank.clear_cache()
2945 sage: J.rank.clear_cache()
2948 sage: J.rank() == J1.rank() + J2.rank()
2951 The same rank computation works over the rationals, with whatever
2954 sage: J1 = HadamardEJA(1, field=QQ, orthonormalize=False)
2955 sage: J2 = RealSymmetricEJA(2, field=QQ, orthonormalize=False)
2956 sage: J = cartesian_product([J1,J2])
2957 sage: J1.rank.clear_cache()
2958 sage: J2.rank.clear_cache()
2959 sage: J.rank.clear_cache()
2962 sage: J.rank() == J1.rank() + J2.rank()
2965 The product algebra will be associative if and only if all of its
2966 components are associative::
2968 sage: J1 = HadamardEJA(2)
2969 sage: J1.is_associative()
2971 sage: J2 = HadamardEJA(3)
2972 sage: J2.is_associative()
2974 sage: J3 = RealSymmetricEJA(3)
2975 sage: J3.is_associative()
2977 sage: CP1 = cartesian_product([J1,J2])
2978 sage: CP1.is_associative()
2980 sage: CP2 = cartesian_product([J1,J3])
2981 sage: CP2.is_associative()
2986 All factors must share the same base field::
2988 sage: J1 = HadamardEJA(2, field=QQ)
2989 sage: J2 = RealSymmetricEJA(2)
2990 sage: CartesianProductEJA((J1,J2))
2991 Traceback (most recent call last):
2993 ValueError: all factors must share the same base field
2995 The cached unit element is the same one that would be computed::
2997 sage: set_random_seed() # long time
2998 sage: J1 = random_eja() # long time
2999 sage: J2 = random_eja() # long time
3000 sage: J = cartesian_product([J1,J2]) # long time
3001 sage: actual = J.one() # long time
3002 sage: J.one.clear_cache() # long time
3003 sage: expected = J.one() # long time
3004 sage: actual == expected # long time
3008 Element
= FiniteDimensionalEJAElement
3011 def __init__(self
, algebras
, **kwargs
):
3012 CombinatorialFreeModule_CartesianProduct
.__init
__(self
,
3015 field
= algebras
[0].base_ring()
3016 if not all( J
.base_ring() == field
for J
in algebras
):
3017 raise ValueError("all factors must share the same base field")
3019 associative
= all( m
.is_associative() for m
in algebras
)
3021 # The definition of matrix_space() and self.basis() relies
3022 # only on the stuff in the CFM_CartesianProduct class, which
3023 # we've already initialized.
3024 Js
= self
.cartesian_factors()
3026 MS
= self
.matrix_space()
3028 MS(tuple( self
.cartesian_projection(i
)(b
).to_matrix()
3029 for i
in range(m
) ))
3030 for b
in self
.basis()
3033 # Define jordan/inner products that operate on that matrix_basis.
3034 def jordan_product(x
,y
):
3036 (Js
[i
](x
[i
])*Js
[i
](y
[i
])).to_matrix() for i
in range(m
)
3039 def inner_product(x
, y
):
3041 Js
[i
](x
[i
]).inner_product(Js
[i
](y
[i
])) for i
in range(m
)
3044 # There's no need to check the field since it already came
3045 # from an EJA. Likewise the axioms are guaranteed to be
3046 # satisfied, unless the guy writing this class sucks.
3048 # If you want the basis to be orthonormalized, orthonormalize
3050 FiniteDimensionalEJA
.__init
__(self
,
3055 orthonormalize
=False,
3056 associative
=associative
,
3057 cartesian_product
=True,
3061 ones
= tuple(J
.one() for J
in algebras
)
3062 self
.one
.set_cache(self
._cartesian
_product
_of
_elements
(ones
))
3063 self
.rank
.set_cache(sum(J
.rank() for J
in algebras
))
3065 def matrix_space(self
):
3067 Return the space that our matrix basis lives in as a Cartesian
3072 sage: from mjo.eja.eja_algebra import (HadamardEJA,
3073 ....: RealSymmetricEJA)
3077 sage: J1 = HadamardEJA(1)
3078 sage: J2 = RealSymmetricEJA(2)
3079 sage: J = cartesian_product([J1,J2])
3080 sage: J.matrix_space()
3081 The Cartesian product of (Full MatrixSpace of 1 by 1 dense
3082 matrices over Algebraic Real Field, Full MatrixSpace of 2
3083 by 2 dense matrices over Algebraic Real Field)
3086 from sage
.categories
.cartesian_product
import cartesian_product
3087 return cartesian_product( [J
.matrix_space()
3088 for J
in self
.cartesian_factors()] )
3091 def cartesian_projection(self
, i
):
3095 sage: from mjo.eja.eja_algebra import (random_eja,
3096 ....: JordanSpinEJA,
3098 ....: RealSymmetricEJA,
3099 ....: ComplexHermitianEJA)
3103 The projection morphisms are Euclidean Jordan algebra
3106 sage: J1 = HadamardEJA(2)
3107 sage: J2 = RealSymmetricEJA(2)
3108 sage: J = cartesian_product([J1,J2])
3109 sage: J.cartesian_projection(0)
3110 Linear operator between finite-dimensional Euclidean Jordan
3111 algebras represented by the matrix:
3114 Domain: Euclidean Jordan algebra of dimension 2 over Algebraic
3115 Real Field (+) Euclidean Jordan algebra of dimension 3 over
3116 Algebraic Real Field
3117 Codomain: Euclidean Jordan algebra of dimension 2 over Algebraic
3119 sage: J.cartesian_projection(1)
3120 Linear operator between finite-dimensional Euclidean Jordan
3121 algebras represented by the matrix:
3125 Domain: Euclidean Jordan algebra of dimension 2 over Algebraic
3126 Real Field (+) Euclidean Jordan algebra of dimension 3 over
3127 Algebraic Real Field
3128 Codomain: Euclidean Jordan algebra of dimension 3 over Algebraic
3131 The projections work the way you'd expect on the vector
3132 representation of an element::
3134 sage: J1 = JordanSpinEJA(2)
3135 sage: J2 = ComplexHermitianEJA(2)
3136 sage: J = cartesian_product([J1,J2])
3137 sage: pi_left = J.cartesian_projection(0)
3138 sage: pi_right = J.cartesian_projection(1)
3139 sage: pi_left(J.one()).to_vector()
3141 sage: pi_right(J.one()).to_vector()
3143 sage: J.one().to_vector()
3148 The answer never changes::
3150 sage: set_random_seed()
3151 sage: J1 = random_eja()
3152 sage: J2 = random_eja()
3153 sage: J = cartesian_product([J1,J2])
3154 sage: P0 = J.cartesian_projection(0)
3155 sage: P1 = J.cartesian_projection(0)
3160 Ji
= self
.cartesian_factors()[i
]
3161 # Requires the fix on Trac 31421/31422 to work!
3162 Pi
= super().cartesian_projection(i
)
3163 return FiniteDimensionalEJAOperator(self
,Ji
,Pi
.matrix())
3166 def cartesian_embedding(self
, i
):
3170 sage: from mjo.eja.eja_algebra import (random_eja,
3171 ....: JordanSpinEJA,
3173 ....: RealSymmetricEJA)
3177 The embedding morphisms are Euclidean Jordan algebra
3180 sage: J1 = HadamardEJA(2)
3181 sage: J2 = RealSymmetricEJA(2)
3182 sage: J = cartesian_product([J1,J2])
3183 sage: J.cartesian_embedding(0)
3184 Linear operator between finite-dimensional Euclidean Jordan
3185 algebras represented by the matrix:
3191 Domain: Euclidean Jordan algebra of dimension 2 over
3192 Algebraic Real Field
3193 Codomain: Euclidean Jordan algebra of dimension 2 over
3194 Algebraic Real Field (+) Euclidean Jordan algebra of
3195 dimension 3 over Algebraic Real Field
3196 sage: J.cartesian_embedding(1)
3197 Linear operator between finite-dimensional Euclidean Jordan
3198 algebras represented by the matrix:
3204 Domain: Euclidean Jordan algebra of dimension 3 over
3205 Algebraic Real Field
3206 Codomain: Euclidean Jordan algebra of dimension 2 over
3207 Algebraic Real Field (+) Euclidean Jordan algebra of
3208 dimension 3 over Algebraic Real Field
3210 The embeddings work the way you'd expect on the vector
3211 representation of an element::
3213 sage: J1 = JordanSpinEJA(3)
3214 sage: J2 = RealSymmetricEJA(2)
3215 sage: J = cartesian_product([J1,J2])
3216 sage: iota_left = J.cartesian_embedding(0)
3217 sage: iota_right = J.cartesian_embedding(1)
3218 sage: iota_left(J1.zero()) == J.zero()
3220 sage: iota_right(J2.zero()) == J.zero()
3222 sage: J1.one().to_vector()
3224 sage: iota_left(J1.one()).to_vector()
3226 sage: J2.one().to_vector()
3228 sage: iota_right(J2.one()).to_vector()
3230 sage: J.one().to_vector()
3235 The answer never changes::
3237 sage: set_random_seed()
3238 sage: J1 = random_eja()
3239 sage: J2 = random_eja()
3240 sage: J = cartesian_product([J1,J2])
3241 sage: E0 = J.cartesian_embedding(0)
3242 sage: E1 = J.cartesian_embedding(0)
3246 Composing a projection with the corresponding inclusion should
3247 produce the identity map, and mismatching them should produce
3250 sage: set_random_seed()
3251 sage: J1 = random_eja()
3252 sage: J2 = random_eja()
3253 sage: J = cartesian_product([J1,J2])
3254 sage: iota_left = J.cartesian_embedding(0)
3255 sage: iota_right = J.cartesian_embedding(1)
3256 sage: pi_left = J.cartesian_projection(0)
3257 sage: pi_right = J.cartesian_projection(1)
3258 sage: pi_left*iota_left == J1.one().operator()
3260 sage: pi_right*iota_right == J2.one().operator()
3262 sage: (pi_left*iota_right).is_zero()
3264 sage: (pi_right*iota_left).is_zero()
3268 Ji
= self
.cartesian_factors()[i
]
3269 # Requires the fix on Trac 31421/31422 to work!
3270 Ei
= super().cartesian_embedding(i
)
3271 return FiniteDimensionalEJAOperator(Ji
,self
,Ei
.matrix())
3275 FiniteDimensionalEJA
.CartesianProduct
= CartesianProductEJA
3277 class RationalBasisCartesianProductEJA(CartesianProductEJA
,
3280 A separate class for products of algebras for which we know a
3285 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
3286 ....: RealSymmetricEJA)
3290 This gives us fast characteristic polynomial computations in
3291 product algebras, too::
3294 sage: J1 = JordanSpinEJA(2)
3295 sage: J2 = RealSymmetricEJA(3)
3296 sage: J = cartesian_product([J1,J2])
3297 sage: J.characteristic_polynomial_of().degree()
3303 def __init__(self
, algebras
, **kwargs
):
3304 CartesianProductEJA
.__init
__(self
, algebras
, **kwargs
)
3306 self
._rational
_algebra
= None
3307 if self
.vector_space().base_field() is not QQ
:
3308 self
._rational
_algebra
= cartesian_product([
3309 r
._rational
_algebra
for r
in algebras
3313 RationalBasisEJA
.CartesianProduct
= RationalBasisCartesianProductEJA
3315 random_eja
= ConcreteEJA
.random_instance