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()
183 # Element subalgebras can take advantage of this.
184 category
= category
.Associative()
185 if cartesian_product
:
186 category
= category
.CartesianProducts()
188 # Call the superclass constructor so that we can use its from_vector()
189 # method to build our multiplication table.
190 CombinatorialFreeModule
.__init
__(self
,
197 # Now comes all of the hard work. We'll be constructing an
198 # ambient vector space V that our (vectorized) basis lives in,
199 # as well as a subspace W of V spanned by those (vectorized)
200 # basis elements. The W-coordinates are the coefficients that
201 # we see in things like x = 1*e1 + 2*e2.
206 degree
= len(_all2list(basis
[0]))
208 # Build an ambient space that fits our matrix basis when
209 # written out as "long vectors."
210 V
= VectorSpace(field
, degree
)
212 # The matrix that will hole the orthonormal -> unorthonormal
213 # coordinate transformation.
214 self
._deortho
_matrix
= None
217 # Save a copy of the un-orthonormalized basis for later.
218 # Convert it to ambient V (vector) coordinates while we're
219 # at it, because we'd have to do it later anyway.
220 deortho_vector_basis
= tuple( V(_all2list(b
)) for b
in basis
)
222 from mjo
.eja
.eja_utils
import gram_schmidt
223 basis
= tuple(gram_schmidt(basis
, inner_product
))
225 # Save the (possibly orthonormalized) matrix basis for
227 self
._matrix
_basis
= basis
229 # Now create the vector space for the algebra, which will have
230 # its own set of non-ambient coordinates (in terms of the
232 vector_basis
= tuple( V(_all2list(b
)) for b
in basis
)
233 W
= V
.span_of_basis( vector_basis
, check
=check_axioms
)
236 # Now "W" is the vector space of our algebra coordinates. The
237 # variables "X1", "X2",... refer to the entries of vectors in
238 # W. Thus to convert back and forth between the orthonormal
239 # coordinates and the given ones, we need to stick the original
241 U
= V
.span_of_basis( deortho_vector_basis
, check
=check_axioms
)
242 self
._deortho
_matrix
= matrix( U
.coordinate_vector(q
)
243 for q
in vector_basis
)
246 # Now we actually compute the multiplication and inner-product
247 # tables/matrices using the possibly-orthonormalized basis.
248 self
._inner
_product
_matrix
= matrix
.identity(field
, n
)
249 self
._multiplication
_table
= [ [0 for j
in range(i
+1)]
252 # Note: the Jordan and inner-products are defined in terms
253 # of the ambient basis. It's important that their arguments
254 # are in ambient coordinates as well.
257 # ortho basis w.r.t. ambient coords
261 # The jordan product returns a matrixy answer, so we
262 # have to convert it to the algebra coordinates.
263 elt
= jordan_product(q_i
, q_j
)
264 elt
= W
.coordinate_vector(V(_all2list(elt
)))
265 self
._multiplication
_table
[i
][j
] = self
.from_vector(elt
)
267 if not orthonormalize
:
268 # If we're orthonormalizing the basis with respect
269 # to an inner-product, then the inner-product
270 # matrix with respect to the resulting basis is
271 # just going to be the identity.
272 ip
= inner_product(q_i
, q_j
)
273 self
._inner
_product
_matrix
[i
,j
] = ip
274 self
._inner
_product
_matrix
[j
,i
] = ip
276 self
._inner
_product
_matrix
._cache
= {'hermitian': True}
277 self
._inner
_product
_matrix
.set_immutable()
280 if not self
._is
_jordanian
():
281 raise ValueError("Jordan identity does not hold")
282 if not self
._inner
_product
_is
_associative
():
283 raise ValueError("inner product is not associative")
286 def _coerce_map_from_base_ring(self
):
288 Disable the map from the base ring into the algebra.
290 Performing a nonsense conversion like this automatically
291 is counterpedagogical. The fallback is to try the usual
292 element constructor, which should also fail.
296 sage: from mjo.eja.eja_algebra import random_eja
300 sage: set_random_seed()
301 sage: J = random_eja()
303 Traceback (most recent call last):
305 ValueError: not an element of this algebra
311 def product_on_basis(self
, i
, j
):
313 Returns the Jordan product of the `i` and `j`th basis elements.
315 This completely defines the Jordan product on the algebra, and
316 is used direclty by our superclass machinery to implement
321 sage: from mjo.eja.eja_algebra import random_eja
325 sage: set_random_seed()
326 sage: J = random_eja()
327 sage: n = J.dimension()
330 sage: ei_ej = J.zero()*J.zero()
332 ....: i = ZZ.random_element(n)
333 ....: j = ZZ.random_element(n)
334 ....: ei = J.gens()[i]
335 ....: ej = J.gens()[j]
336 ....: ei_ej = J.product_on_basis(i,j)
341 # We only stored the lower-triangular portion of the
342 # multiplication table.
344 return self
._multiplication
_table
[i
][j
]
346 return self
._multiplication
_table
[j
][i
]
348 def inner_product(self
, x
, y
):
350 The inner product associated with this Euclidean Jordan algebra.
352 Defaults to the trace inner product, but can be overridden by
353 subclasses if they are sure that the necessary properties are
358 sage: from mjo.eja.eja_algebra import (random_eja,
360 ....: BilinearFormEJA)
364 Our inner product is "associative," which means the following for
365 a symmetric bilinear form::
367 sage: set_random_seed()
368 sage: J = random_eja()
369 sage: x,y,z = J.random_elements(3)
370 sage: (x*y).inner_product(z) == y.inner_product(x*z)
375 Ensure that this is the usual inner product for the algebras
378 sage: set_random_seed()
379 sage: J = HadamardEJA.random_instance()
380 sage: x,y = J.random_elements(2)
381 sage: actual = x.inner_product(y)
382 sage: expected = x.to_vector().inner_product(y.to_vector())
383 sage: actual == expected
386 Ensure that this is one-half of the trace inner-product in a
387 BilinearFormEJA that isn't just the reals (when ``n`` isn't
388 one). This is in Faraut and Koranyi, and also my "On the
391 sage: set_random_seed()
392 sage: J = BilinearFormEJA.random_instance()
393 sage: n = J.dimension()
394 sage: x = J.random_element()
395 sage: y = J.random_element()
396 sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
400 B
= self
._inner
_product
_matrix
401 return (B
*x
.to_vector()).inner_product(y
.to_vector())
404 def is_associative(self
):
406 Return whether or not this algebra's Jordan product is associative.
410 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
414 sage: J = ComplexHermitianEJA(3, field=QQ, orthonormalize=False)
415 sage: J.is_associative()
417 sage: x = sum(J.gens())
418 sage: A = x.subalgebra_generated_by(orthonormalize=False)
419 sage: A.is_associative()
423 return "Associative" in self
.category().axioms()
425 def _is_jordanian(self
):
427 Whether or not this algebra's multiplication table respects the
428 Jordan identity `(x^{2})(xy) = x(x^{2}y)`.
430 We only check one arrangement of `x` and `y`, so for a
431 ``True`` result to be truly true, you should also check
432 :meth:`is_commutative`. This method should of course always
433 return ``True``, unless this algebra was constructed with
434 ``check_axioms=False`` and passed an invalid multiplication table.
436 return all( (self
.gens()[i
]**2)*(self
.gens()[i
]*self
.gens()[j
])
438 (self
.gens()[i
])*((self
.gens()[i
]**2)*self
.gens()[j
])
439 for i
in range(self
.dimension())
440 for j
in range(self
.dimension()) )
442 def _inner_product_is_associative(self
):
444 Return whether or not this algebra's inner product `B` is
445 associative; that is, whether or not `B(xy,z) = B(x,yz)`.
447 This method should of course always return ``True``, unless
448 this algebra was constructed with ``check_axioms=False`` and
449 passed an invalid Jordan or inner-product.
452 # Used to check whether or not something is zero in an inexact
453 # ring. This number is sufficient to allow the construction of
454 # QuaternionHermitianEJA(2, field=RDF) with check_axioms=True.
457 for i
in range(self
.dimension()):
458 for j
in range(self
.dimension()):
459 for k
in range(self
.dimension()):
463 diff
= (x
*y
).inner_product(z
) - x
.inner_product(y
*z
)
465 if self
.base_ring().is_exact():
469 if diff
.abs() > epsilon
:
474 def _element_constructor_(self
, elt
):
476 Construct an element of this algebra from its vector or matrix
479 This gets called only after the parent element _call_ method
480 fails to find a coercion for the argument.
484 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
486 ....: RealSymmetricEJA)
490 The identity in `S^n` is converted to the identity in the EJA::
492 sage: J = RealSymmetricEJA(3)
493 sage: I = matrix.identity(QQ,3)
494 sage: J(I) == J.one()
497 This skew-symmetric matrix can't be represented in the EJA::
499 sage: J = RealSymmetricEJA(3)
500 sage: A = matrix(QQ,3, lambda i,j: i-j)
502 Traceback (most recent call last):
504 ValueError: not an element of this algebra
506 Tuples work as well, provided that the matrix basis for the
507 algebra consists of them::
509 sage: J1 = HadamardEJA(3)
510 sage: J2 = RealSymmetricEJA(2)
511 sage: J = cartesian_product([J1,J2])
512 sage: J( (J1.matrix_basis()[1], J2.matrix_basis()[2]) )
517 Ensure that we can convert any element of the two non-matrix
518 simple algebras (whose matrix representations are columns)
519 back and forth faithfully::
521 sage: set_random_seed()
522 sage: J = HadamardEJA.random_instance()
523 sage: x = J.random_element()
524 sage: J(x.to_vector().column()) == x
526 sage: J = JordanSpinEJA.random_instance()
527 sage: x = J.random_element()
528 sage: J(x.to_vector().column()) == x
531 We cannot coerce elements between algebras just because their
532 matrix representations are compatible::
534 sage: J1 = HadamardEJA(3)
535 sage: J2 = JordanSpinEJA(3)
537 Traceback (most recent call last):
539 ValueError: not an element of this algebra
541 Traceback (most recent call last):
543 ValueError: not an element of this algebra
546 msg
= "not an element of this algebra"
547 if elt
in self
.base_ring():
548 # Ensure that no base ring -> algebra coercion is performed
549 # by this method. There's some stupidity in sage that would
550 # otherwise propagate to this method; for example, sage thinks
551 # that the integer 3 belongs to the space of 2-by-2 matrices.
552 raise ValueError(msg
)
555 # Try to convert a vector into a column-matrix...
557 except (AttributeError, TypeError):
558 # and ignore failure, because we weren't really expecting
559 # a vector as an argument anyway.
562 if elt
not in self
.matrix_space():
563 raise ValueError(msg
)
565 # Thanks for nothing! Matrix spaces aren't vector spaces in
566 # Sage, so we have to figure out its matrix-basis coordinates
567 # ourselves. We use the basis space's ring instead of the
568 # element's ring because the basis space might be an algebraic
569 # closure whereas the base ring of the 3-by-3 identity matrix
570 # could be QQ instead of QQbar.
572 # And, we also have to handle Cartesian product bases (when
573 # the matrix basis consists of tuples) here. The "good news"
574 # is that we're already converting everything to long vectors,
575 # and that strategy works for tuples as well.
577 # We pass check=False because the matrix basis is "guaranteed"
578 # to be linearly independent... right? Ha ha.
580 V
= VectorSpace(self
.base_ring(), len(elt
))
581 W
= V
.span_of_basis( (V(_all2list(s
)) for s
in self
.matrix_basis()),
585 coords
= W
.coordinate_vector(V(elt
))
586 except ArithmeticError: # vector is not in free module
587 raise ValueError(msg
)
589 return self
.from_vector(coords
)
593 Return a string representation of ``self``.
597 sage: from mjo.eja.eja_algebra import JordanSpinEJA
601 Ensure that it says what we think it says::
603 sage: JordanSpinEJA(2, field=AA)
604 Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
605 sage: JordanSpinEJA(3, field=RDF)
606 Euclidean Jordan algebra of dimension 3 over Real Double Field
609 fmt
= "Euclidean Jordan algebra of dimension {} over {}"
610 return fmt
.format(self
.dimension(), self
.base_ring())
614 def characteristic_polynomial_of(self
):
616 Return the algebra's "characteristic polynomial of" function,
617 which is itself a multivariate polynomial that, when evaluated
618 at the coordinates of some algebra element, returns that
619 element's characteristic polynomial.
621 The resulting polynomial has `n+1` variables, where `n` is the
622 dimension of this algebra. The first `n` variables correspond to
623 the coordinates of an algebra element: when evaluated at the
624 coordinates of an algebra element with respect to a certain
625 basis, the result is a univariate polynomial (in the one
626 remaining variable ``t``), namely the characteristic polynomial
631 sage: from mjo.eja.eja_algebra import JordanSpinEJA, TrivialEJA
635 The characteristic polynomial in the spin algebra is given in
636 Alizadeh, Example 11.11::
638 sage: J = JordanSpinEJA(3)
639 sage: p = J.characteristic_polynomial_of(); p
640 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
641 sage: xvec = J.one().to_vector()
645 By definition, the characteristic polynomial is a monic
646 degree-zero polynomial in a rank-zero algebra. Note that
647 Cayley-Hamilton is indeed satisfied since the polynomial
648 ``1`` evaluates to the identity element of the algebra on
651 sage: J = TrivialEJA()
652 sage: J.characteristic_polynomial_of()
659 # The list of coefficient polynomials a_0, a_1, a_2, ..., a_(r-1).
660 a
= self
._charpoly
_coefficients
()
662 # We go to a bit of trouble here to reorder the
663 # indeterminates, so that it's easier to evaluate the
664 # characteristic polynomial at x's coordinates and get back
665 # something in terms of t, which is what we want.
666 S
= PolynomialRing(self
.base_ring(),'t')
670 S
= PolynomialRing(S
, R
.variable_names())
673 return (t
**r
+ sum( a
[k
]*(t
**k
) for k
in range(r
) ))
675 def coordinate_polynomial_ring(self
):
677 The multivariate polynomial ring in which this algebra's
678 :meth:`characteristic_polynomial_of` lives.
682 sage: from mjo.eja.eja_algebra import (HadamardEJA,
683 ....: RealSymmetricEJA)
687 sage: J = HadamardEJA(2)
688 sage: J.coordinate_polynomial_ring()
689 Multivariate Polynomial Ring in X1, X2...
690 sage: J = RealSymmetricEJA(3,field=QQ,orthonormalize=False)
691 sage: J.coordinate_polynomial_ring()
692 Multivariate Polynomial Ring in X1, X2, X3, X4, X5, X6...
695 var_names
= tuple( "X%d" % z
for z
in range(1, self
.dimension()+1) )
696 return PolynomialRing(self
.base_ring(), var_names
)
698 def inner_product(self
, x
, y
):
700 The inner product associated with this Euclidean Jordan algebra.
702 Defaults to the trace inner product, but can be overridden by
703 subclasses if they are sure that the necessary properties are
708 sage: from mjo.eja.eja_algebra import (random_eja,
710 ....: BilinearFormEJA)
714 Our inner product is "associative," which means the following for
715 a symmetric bilinear form::
717 sage: set_random_seed()
718 sage: J = random_eja()
719 sage: x,y,z = J.random_elements(3)
720 sage: (x*y).inner_product(z) == y.inner_product(x*z)
725 Ensure that this is the usual inner product for the algebras
728 sage: set_random_seed()
729 sage: J = HadamardEJA.random_instance()
730 sage: x,y = J.random_elements(2)
731 sage: actual = x.inner_product(y)
732 sage: expected = x.to_vector().inner_product(y.to_vector())
733 sage: actual == expected
736 Ensure that this is one-half of the trace inner-product in a
737 BilinearFormEJA that isn't just the reals (when ``n`` isn't
738 one). This is in Faraut and Koranyi, and also my "On the
741 sage: set_random_seed()
742 sage: J = BilinearFormEJA.random_instance()
743 sage: n = J.dimension()
744 sage: x = J.random_element()
745 sage: y = J.random_element()
746 sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
749 B
= self
._inner
_product
_matrix
750 return (B
*x
.to_vector()).inner_product(y
.to_vector())
753 def is_trivial(self
):
755 Return whether or not this algebra is trivial.
757 A trivial algebra contains only the zero element.
761 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
766 sage: J = ComplexHermitianEJA(3)
772 sage: J = TrivialEJA()
777 return self
.dimension() == 0
780 def multiplication_table(self
):
782 Return a visual representation of this algebra's multiplication
783 table (on basis elements).
787 sage: from mjo.eja.eja_algebra import JordanSpinEJA
791 sage: J = JordanSpinEJA(4)
792 sage: J.multiplication_table()
793 +----++----+----+----+----+
794 | * || e0 | e1 | e2 | e3 |
795 +====++====+====+====+====+
796 | e0 || e0 | e1 | e2 | e3 |
797 +----++----+----+----+----+
798 | e1 || e1 | e0 | 0 | 0 |
799 +----++----+----+----+----+
800 | e2 || e2 | 0 | e0 | 0 |
801 +----++----+----+----+----+
802 | e3 || e3 | 0 | 0 | e0 |
803 +----++----+----+----+----+
807 # Prepend the header row.
808 M
= [["*"] + list(self
.gens())]
810 # And to each subsequent row, prepend an entry that belongs to
811 # the left-side "header column."
812 M
+= [ [self
.gens()[i
]] + [ self
.product_on_basis(i
,j
)
816 return table(M
, header_row
=True, header_column
=True, frame
=True)
819 def matrix_basis(self
):
821 Return an (often more natural) representation of this algebras
822 basis as an ordered tuple of matrices.
824 Every finite-dimensional Euclidean Jordan Algebra is a, up to
825 Jordan isomorphism, a direct sum of five simple
826 algebras---four of which comprise Hermitian matrices. And the
827 last type of algebra can of course be thought of as `n`-by-`1`
828 column matrices (ambiguusly called column vectors) to avoid
829 special cases. As a result, matrices (and column vectors) are
830 a natural representation format for Euclidean Jordan algebra
833 But, when we construct an algebra from a basis of matrices,
834 those matrix representations are lost in favor of coordinate
835 vectors *with respect to* that basis. We could eventually
836 convert back if we tried hard enough, but having the original
837 representations handy is valuable enough that we simply store
838 them and return them from this method.
840 Why implement this for non-matrix algebras? Avoiding special
841 cases for the :class:`BilinearFormEJA` pays with simplicity in
842 its own right. But mainly, we would like to be able to assume
843 that elements of a :class:`CartesianProductEJA` can be displayed
844 nicely, without having to have special classes for direct sums
845 one of whose components was a matrix algebra.
849 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
850 ....: RealSymmetricEJA)
854 sage: J = RealSymmetricEJA(2)
856 Finite family {0: e0, 1: e1, 2: e2}
857 sage: J.matrix_basis()
859 [1 0] [ 0 0.7071067811865475?] [0 0]
860 [0 0], [0.7071067811865475? 0], [0 1]
865 sage: J = JordanSpinEJA(2)
867 Finite family {0: e0, 1: e1}
868 sage: J.matrix_basis()
874 return self
._matrix
_basis
877 def matrix_space(self
):
879 Return the matrix space in which this algebra's elements live, if
880 we think of them as matrices (including column vectors of the
883 "By default" this will be an `n`-by-`1` column-matrix space,
884 except when the algebra is trivial. There it's `n`-by-`n`
885 (where `n` is zero), to ensure that two elements of the matrix
886 space (empty matrices) can be multiplied. For algebras of
887 matrices, this returns the space in which their
888 real embeddings live.
892 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
894 ....: QuaternionHermitianEJA,
899 By default, the matrix representation is just a column-matrix
900 equivalent to the vector representation::
902 sage: J = JordanSpinEJA(3)
903 sage: J.matrix_space()
904 Full MatrixSpace of 3 by 1 dense matrices over Algebraic
907 The matrix representation in the trivial algebra is
908 zero-by-zero instead of the usual `n`-by-one::
910 sage: J = TrivialEJA()
911 sage: J.matrix_space()
912 Full MatrixSpace of 0 by 0 dense matrices over Algebraic
915 The matrix space for complex/quaternion Hermitian matrix EJA
916 is the space in which their real-embeddings live, not the
917 original complex/quaternion matrix space::
919 sage: J = ComplexHermitianEJA(2,field=QQ,orthonormalize=False)
920 sage: J.matrix_space()
921 Full MatrixSpace of 4 by 4 dense matrices over Rational Field
922 sage: J = QuaternionHermitianEJA(1,field=QQ,orthonormalize=False)
923 sage: J.matrix_space()
924 Full MatrixSpace of 4 by 4 dense matrices over Rational Field
927 if self
.is_trivial():
928 return MatrixSpace(self
.base_ring(), 0)
930 return self
.matrix_basis()[0].parent()
936 Return the unit element of this algebra.
940 sage: from mjo.eja.eja_algebra import (HadamardEJA,
945 We can compute unit element in the Hadamard EJA::
947 sage: J = HadamardEJA(5)
949 e0 + e1 + e2 + e3 + e4
951 The unit element in the Hadamard EJA is inherited in the
952 subalgebras generated by its elements::
954 sage: J = HadamardEJA(5)
956 e0 + e1 + e2 + e3 + e4
957 sage: x = sum(J.gens())
958 sage: A = x.subalgebra_generated_by(orthonormalize=False)
961 sage: A.one().superalgebra_element()
962 e0 + e1 + e2 + e3 + e4
966 The identity element acts like the identity, regardless of
967 whether or not we orthonormalize::
969 sage: set_random_seed()
970 sage: J = random_eja()
971 sage: x = J.random_element()
972 sage: J.one()*x == x and x*J.one() == x
974 sage: A = x.subalgebra_generated_by()
975 sage: y = A.random_element()
976 sage: A.one()*y == y and y*A.one() == y
981 sage: set_random_seed()
982 sage: J = random_eja(field=QQ, orthonormalize=False)
983 sage: x = J.random_element()
984 sage: J.one()*x == x and x*J.one() == x
986 sage: A = x.subalgebra_generated_by(orthonormalize=False)
987 sage: y = A.random_element()
988 sage: A.one()*y == y and y*A.one() == y
991 The matrix of the unit element's operator is the identity,
992 regardless of the base field and whether or not we
995 sage: set_random_seed()
996 sage: J = random_eja()
997 sage: actual = J.one().operator().matrix()
998 sage: expected = matrix.identity(J.base_ring(), J.dimension())
999 sage: actual == expected
1001 sage: x = J.random_element()
1002 sage: A = x.subalgebra_generated_by()
1003 sage: actual = A.one().operator().matrix()
1004 sage: expected = matrix.identity(A.base_ring(), A.dimension())
1005 sage: actual == expected
1010 sage: set_random_seed()
1011 sage: J = random_eja(field=QQ, orthonormalize=False)
1012 sage: actual = J.one().operator().matrix()
1013 sage: expected = matrix.identity(J.base_ring(), J.dimension())
1014 sage: actual == expected
1016 sage: x = J.random_element()
1017 sage: A = x.subalgebra_generated_by(orthonormalize=False)
1018 sage: actual = A.one().operator().matrix()
1019 sage: expected = matrix.identity(A.base_ring(), A.dimension())
1020 sage: actual == expected
1023 Ensure that the cached unit element (often precomputed by
1024 hand) agrees with the computed one::
1026 sage: set_random_seed()
1027 sage: J = random_eja()
1028 sage: cached = J.one()
1029 sage: J.one.clear_cache()
1030 sage: J.one() == cached
1035 sage: set_random_seed()
1036 sage: J = random_eja(field=QQ, orthonormalize=False)
1037 sage: cached = J.one()
1038 sage: J.one.clear_cache()
1039 sage: J.one() == cached
1043 # We can brute-force compute the matrices of the operators
1044 # that correspond to the basis elements of this algebra.
1045 # If some linear combination of those basis elements is the
1046 # algebra identity, then the same linear combination of
1047 # their matrices has to be the identity matrix.
1049 # Of course, matrices aren't vectors in sage, so we have to
1050 # appeal to the "long vectors" isometry.
1051 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
1053 # Now we use basic linear algebra to find the coefficients,
1054 # of the matrices-as-vectors-linear-combination, which should
1055 # work for the original algebra basis too.
1056 A
= matrix(self
.base_ring(), oper_vecs
)
1058 # We used the isometry on the left-hand side already, but we
1059 # still need to do it for the right-hand side. Recall that we
1060 # wanted something that summed to the identity matrix.
1061 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
1063 # Now if there's an identity element in the algebra, this
1064 # should work. We solve on the left to avoid having to
1065 # transpose the matrix "A".
1066 return self
.from_vector(A
.solve_left(b
))
1069 def peirce_decomposition(self
, c
):
1071 The Peirce decomposition of this algebra relative to the
1074 In the future, this can be extended to a complete system of
1075 orthogonal idempotents.
1079 - ``c`` -- an idempotent of this algebra.
1083 A triple (J0, J5, J1) containing two subalgebras and one subspace
1086 - ``J0`` -- the algebra on the eigenspace of ``c.operator()``
1087 corresponding to the eigenvalue zero.
1089 - ``J5`` -- the eigenspace (NOT a subalgebra) of ``c.operator()``
1090 corresponding to the eigenvalue one-half.
1092 - ``J1`` -- the algebra on the eigenspace of ``c.operator()``
1093 corresponding to the eigenvalue one.
1095 These are the only possible eigenspaces for that operator, and this
1096 algebra is a direct sum of them. The spaces ``J0`` and ``J1`` are
1097 orthogonal, and are subalgebras of this algebra with the appropriate
1102 sage: from mjo.eja.eja_algebra import random_eja, RealSymmetricEJA
1106 The canonical example comes from the symmetric matrices, which
1107 decompose into diagonal and off-diagonal parts::
1109 sage: J = RealSymmetricEJA(3)
1110 sage: C = matrix(QQ, [ [1,0,0],
1114 sage: J0,J5,J1 = J.peirce_decomposition(c)
1116 Euclidean Jordan algebra of dimension 1...
1118 Vector space of degree 6 and dimension 2...
1120 Euclidean Jordan algebra of dimension 3...
1121 sage: J0.one().to_matrix()
1125 sage: orig_df = AA.options.display_format
1126 sage: AA.options.display_format = 'radical'
1127 sage: J.from_vector(J5.basis()[0]).to_matrix()
1131 sage: J.from_vector(J5.basis()[1]).to_matrix()
1135 sage: AA.options.display_format = orig_df
1136 sage: J1.one().to_matrix()
1143 Every algebra decomposes trivially with respect to its identity
1146 sage: set_random_seed()
1147 sage: J = random_eja()
1148 sage: J0,J5,J1 = J.peirce_decomposition(J.one())
1149 sage: J0.dimension() == 0 and J5.dimension() == 0
1151 sage: J1.superalgebra() == J and J1.dimension() == J.dimension()
1154 The decomposition is into eigenspaces, and its components are
1155 therefore necessarily orthogonal. Moreover, the identity
1156 elements in the two subalgebras are the projections onto their
1157 respective subspaces of the superalgebra's identity element::
1159 sage: set_random_seed()
1160 sage: J = random_eja()
1161 sage: x = J.random_element()
1162 sage: if not J.is_trivial():
1163 ....: while x.is_nilpotent():
1164 ....: x = J.random_element()
1165 sage: c = x.subalgebra_idempotent()
1166 sage: J0,J5,J1 = J.peirce_decomposition(c)
1168 sage: for (w,y,z) in zip(J0.basis(), J5.basis(), J1.basis()):
1169 ....: w = w.superalgebra_element()
1170 ....: y = J.from_vector(y)
1171 ....: z = z.superalgebra_element()
1172 ....: ipsum += w.inner_product(y).abs()
1173 ....: ipsum += w.inner_product(z).abs()
1174 ....: ipsum += y.inner_product(z).abs()
1177 sage: J1(c) == J1.one()
1179 sage: J0(J.one() - c) == J0.one()
1183 if not c
.is_idempotent():
1184 raise ValueError("element is not idempotent: %s" % c
)
1186 # Default these to what they should be if they turn out to be
1187 # trivial, because eigenspaces_left() won't return eigenvalues
1188 # corresponding to trivial spaces (e.g. it returns only the
1189 # eigenspace corresponding to lambda=1 if you take the
1190 # decomposition relative to the identity element).
1191 trivial
= self
.subalgebra(())
1192 J0
= trivial
# eigenvalue zero
1193 J5
= VectorSpace(self
.base_ring(), 0) # eigenvalue one-half
1194 J1
= trivial
# eigenvalue one
1196 for (eigval
, eigspace
) in c
.operator().matrix().right_eigenspaces():
1197 if eigval
== ~
(self
.base_ring()(2)):
1200 gens
= tuple( self
.from_vector(b
) for b
in eigspace
.basis() )
1201 subalg
= self
.subalgebra(gens
, check_axioms
=False)
1207 raise ValueError("unexpected eigenvalue: %s" % eigval
)
1212 def random_element(self
, thorough
=False):
1214 Return a random element of this algebra.
1216 Our algebra superclass method only returns a linear
1217 combination of at most two basis elements. We instead
1218 want the vector space "random element" method that
1219 returns a more diverse selection.
1223 - ``thorough`` -- (boolean; default False) whether or not we
1224 should generate irrational coefficients for the random
1225 element when our base ring is irrational; this slows the
1226 algebra operations to a crawl, but any truly random method
1230 # For a general base ring... maybe we can trust this to do the
1231 # right thing? Unlikely, but.
1232 V
= self
.vector_space()
1233 v
= V
.random_element()
1235 if self
.base_ring() is AA
:
1236 # The "random element" method of the algebraic reals is
1237 # stupid at the moment, and only returns integers between
1238 # -2 and 2, inclusive:
1240 # https://trac.sagemath.org/ticket/30875
1242 # Instead, we implement our own "random vector" method,
1243 # and then coerce that into the algebra. We use the vector
1244 # space degree here instead of the dimension because a
1245 # subalgebra could (for example) be spanned by only two
1246 # vectors, each with five coordinates. We need to
1247 # generate all five coordinates.
1249 v
*= QQbar
.random_element().real()
1251 v
*= QQ
.random_element()
1253 return self
.from_vector(V
.coordinate_vector(v
))
1255 def random_elements(self
, count
, thorough
=False):
1257 Return ``count`` random elements as a tuple.
1261 - ``thorough`` -- (boolean; default False) whether or not we
1262 should generate irrational coefficients for the random
1263 elements when our base ring is irrational; this slows the
1264 algebra operations to a crawl, but any truly random method
1269 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1273 sage: J = JordanSpinEJA(3)
1274 sage: x,y,z = J.random_elements(3)
1275 sage: all( [ x in J, y in J, z in J ])
1277 sage: len( J.random_elements(10) ) == 10
1281 return tuple( self
.random_element(thorough
)
1282 for idx
in range(count
) )
1286 def _charpoly_coefficients(self
):
1288 The `r` polynomial coefficients of the "characteristic polynomial
1293 sage: from mjo.eja.eja_algebra import random_eja
1297 The theory shows that these are all homogeneous polynomials of
1300 sage: set_random_seed()
1301 sage: J = random_eja()
1302 sage: all(p.is_homogeneous() for p in J._charpoly_coefficients())
1306 n
= self
.dimension()
1307 R
= self
.coordinate_polynomial_ring()
1309 F
= R
.fraction_field()
1312 # From a result in my book, these are the entries of the
1313 # basis representation of L_x.
1314 return sum( vars[k
]*self
.gens()[k
].operator().matrix()[i
,j
]
1317 L_x
= matrix(F
, n
, n
, L_x_i_j
)
1320 if self
.rank
.is_in_cache():
1322 # There's no need to pad the system with redundant
1323 # columns if we *know* they'll be redundant.
1326 # Compute an extra power in case the rank is equal to
1327 # the dimension (otherwise, we would stop at x^(r-1)).
1328 x_powers
= [ (L_x
**k
)*self
.one().to_vector()
1329 for k
in range(n
+1) ]
1330 A
= matrix
.column(F
, x_powers
[:n
])
1331 AE
= A
.extended_echelon_form()
1338 # The theory says that only the first "r" coefficients are
1339 # nonzero, and they actually live in the original polynomial
1340 # ring and not the fraction field. We negate them because in
1341 # the actual characteristic polynomial, they get moved to the
1342 # other side where x^r lives. We don't bother to trim A_rref
1343 # down to a square matrix and solve the resulting system,
1344 # because the upper-left r-by-r portion of A_rref is
1345 # guaranteed to be the identity matrix, so e.g.
1347 # A_rref.solve_right(Y)
1349 # would just be returning Y.
1350 return (-E
*b
)[:r
].change_ring(R
)
1355 Return the rank of this EJA.
1357 This is a cached method because we know the rank a priori for
1358 all of the algebras we can construct. Thus we can avoid the
1359 expensive ``_charpoly_coefficients()`` call unless we truly
1360 need to compute the whole characteristic polynomial.
1364 sage: from mjo.eja.eja_algebra import (HadamardEJA,
1365 ....: JordanSpinEJA,
1366 ....: RealSymmetricEJA,
1367 ....: ComplexHermitianEJA,
1368 ....: QuaternionHermitianEJA,
1373 The rank of the Jordan spin algebra is always two::
1375 sage: JordanSpinEJA(2).rank()
1377 sage: JordanSpinEJA(3).rank()
1379 sage: JordanSpinEJA(4).rank()
1382 The rank of the `n`-by-`n` Hermitian real, complex, or
1383 quaternion matrices is `n`::
1385 sage: RealSymmetricEJA(4).rank()
1387 sage: ComplexHermitianEJA(3).rank()
1389 sage: QuaternionHermitianEJA(2).rank()
1394 Ensure that every EJA that we know how to construct has a
1395 positive integer rank, unless the algebra is trivial in
1396 which case its rank will be zero::
1398 sage: set_random_seed()
1399 sage: J = random_eja()
1403 sage: r > 0 or (r == 0 and J.is_trivial())
1406 Ensure that computing the rank actually works, since the ranks
1407 of all simple algebras are known and will be cached by default::
1409 sage: set_random_seed() # long time
1410 sage: J = random_eja() # long time
1411 sage: cached = J.rank() # long time
1412 sage: J.rank.clear_cache() # long time
1413 sage: J.rank() == cached # long time
1417 return len(self
._charpoly
_coefficients
())
1420 def subalgebra(self
, basis
, **kwargs
):
1422 Create a subalgebra of this algebra from the given basis.
1424 from mjo
.eja
.eja_subalgebra
import FiniteDimensionalEJASubalgebra
1425 return FiniteDimensionalEJASubalgebra(self
, basis
, **kwargs
)
1428 def vector_space(self
):
1430 Return the vector space that underlies this algebra.
1434 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1438 sage: J = RealSymmetricEJA(2)
1439 sage: J.vector_space()
1440 Vector space of dimension 3 over...
1443 return self
.zero().to_vector().parent().ambient_vector_space()
1447 class RationalBasisEJA(FiniteDimensionalEJA
):
1449 New class for algebras whose supplied basis elements have all rational entries.
1453 sage: from mjo.eja.eja_algebra import BilinearFormEJA
1457 The supplied basis is orthonormalized by default::
1459 sage: B = matrix(QQ, [[1, 0, 0], [0, 25, -32], [0, -32, 41]])
1460 sage: J = BilinearFormEJA(B)
1461 sage: J.matrix_basis()
1478 # Abuse the check_field parameter to check that the entries of
1479 # out basis (in ambient coordinates) are in the field QQ.
1480 if not all( all(b_i
in QQ
for b_i
in b
.list()) for b
in basis
):
1481 raise TypeError("basis not rational")
1483 self
._rational
_algebra
= None
1485 # There's no point in constructing the extra algebra if this
1486 # one is already rational.
1488 # Note: the same Jordan and inner-products work here,
1489 # because they are necessarily defined with respect to
1490 # ambient coordinates and not any particular basis.
1491 self
._rational
_algebra
= FiniteDimensionalEJA(
1496 orthonormalize
=False,
1500 super().__init
__(basis
,
1504 check_field
=check_field
,
1508 def _charpoly_coefficients(self
):
1512 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
1513 ....: JordanSpinEJA)
1517 The base ring of the resulting polynomial coefficients is what
1518 it should be, and not the rationals (unless the algebra was
1519 already over the rationals)::
1521 sage: J = JordanSpinEJA(3)
1522 sage: J._charpoly_coefficients()
1523 (X1^2 - X2^2 - X3^2, -2*X1)
1524 sage: a0 = J._charpoly_coefficients()[0]
1526 Algebraic Real Field
1527 sage: a0.base_ring()
1528 Algebraic Real Field
1531 if self
._rational
_algebra
is None:
1532 # There's no need to construct *another* algebra over the
1533 # rationals if this one is already over the
1534 # rationals. Likewise, if we never orthonormalized our
1535 # basis, we might as well just use the given one.
1536 return super()._charpoly
_coefficients
()
1538 # Do the computation over the rationals. The answer will be
1539 # the same, because all we've done is a change of basis.
1540 # Then, change back from QQ to our real base ring
1541 a
= ( a_i
.change_ring(self
.base_ring())
1542 for a_i
in self
._rational
_algebra
._charpoly
_coefficients
() )
1544 if self
._deortho
_matrix
is None:
1545 # This can happen if our base ring was, say, AA and we
1546 # chose not to (or didn't need to) orthonormalize. It's
1547 # still faster to do the computations over QQ even if
1548 # the numbers in the boxes stay the same.
1551 # Otherwise, convert the coordinate variables back to the
1552 # deorthonormalized ones.
1553 R
= self
.coordinate_polynomial_ring()
1554 from sage
.modules
.free_module_element
import vector
1555 X
= vector(R
, R
.gens())
1556 BX
= self
._deortho
_matrix
*X
1558 subs_dict
= { X[i]: BX[i] for i in range(len(X)) }
1559 return tuple( a_i
.subs(subs_dict
) for a_i
in a
)
1561 class ConcreteEJA(RationalBasisEJA
):
1563 A class for the Euclidean Jordan algebras that we know by name.
1565 These are the Jordan algebras whose basis, multiplication table,
1566 rank, and so on are known a priori. More to the point, they are
1567 the Euclidean Jordan algebras for which we are able to conjure up
1568 a "random instance."
1572 sage: from mjo.eja.eja_algebra import ConcreteEJA
1576 Our basis is normalized with respect to the algebra's inner
1577 product, unless we specify otherwise::
1579 sage: set_random_seed()
1580 sage: J = ConcreteEJA.random_instance()
1581 sage: all( b.norm() == 1 for b in J.gens() )
1584 Since our basis is orthonormal with respect to the algebra's inner
1585 product, and since we know that this algebra is an EJA, any
1586 left-multiplication operator's matrix will be symmetric because
1587 natural->EJA basis representation is an isometry and within the
1588 EJA the operator is self-adjoint by the Jordan axiom::
1590 sage: set_random_seed()
1591 sage: J = ConcreteEJA.random_instance()
1592 sage: x = J.random_element()
1593 sage: x.operator().is_self_adjoint()
1598 def _max_random_instance_size():
1600 Return an integer "size" that is an upper bound on the size of
1601 this algebra when it is used in a random test
1602 case. Unfortunately, the term "size" is ambiguous -- when
1603 dealing with `R^n` under either the Hadamard or Jordan spin
1604 product, the "size" refers to the dimension `n`. When dealing
1605 with a matrix algebra (real symmetric or complex/quaternion
1606 Hermitian), it refers to the size of the matrix, which is far
1607 less than the dimension of the underlying vector space.
1609 This method must be implemented in each subclass.
1611 raise NotImplementedError
1614 def random_instance(cls
, *args
, **kwargs
):
1616 Return a random instance of this type of algebra.
1618 This method should be implemented in each subclass.
1620 from sage
.misc
.prandom
import choice
1621 eja_class
= choice(cls
.__subclasses
__())
1623 # These all bubble up to the RationalBasisEJA superclass
1624 # constructor, so any (kw)args valid there are also valid
1626 return eja_class
.random_instance(*args
, **kwargs
)
1631 def dimension_over_reals():
1633 The dimension of this matrix's base ring over the reals.
1635 The reals are dimension one over themselves, obviously; that's
1636 just `\mathbb{R}^{1}`. Likewise, the complex numbers `a + bi`
1637 have dimension two. Finally, the quaternions have dimension
1638 four over the reals.
1640 This is used to determine the size of the matrix returned from
1641 :meth:`real_embed`, among other things.
1643 raise NotImplementedError
1646 def real_embed(cls
,M
):
1648 Embed the matrix ``M`` into a space of real matrices.
1650 The matrix ``M`` can have entries in any field at the moment:
1651 the real numbers, complex numbers, or quaternions. And although
1652 they are not a field, we can probably support octonions at some
1653 point, too. This function returns a real matrix that "acts like"
1654 the original with respect to matrix multiplication; i.e.
1656 real_embed(M*N) = real_embed(M)*real_embed(N)
1659 if M
.ncols() != M
.nrows():
1660 raise ValueError("the matrix 'M' must be square")
1665 def real_unembed(cls
,M
):
1667 The inverse of :meth:`real_embed`.
1669 if M
.ncols() != M
.nrows():
1670 raise ValueError("the matrix 'M' must be square")
1671 if not ZZ(M
.nrows()).mod(cls
.dimension_over_reals()).is_zero():
1672 raise ValueError("the matrix 'M' must be a real embedding")
1676 def jordan_product(X
,Y
):
1677 return (X
*Y
+ Y
*X
)/2
1680 def trace_inner_product(cls
,X
,Y
):
1682 Compute the trace inner-product of two real-embeddings.
1686 sage: from mjo.eja.eja_algebra import (RealSymmetricEJA,
1687 ....: ComplexHermitianEJA,
1688 ....: QuaternionHermitianEJA)
1692 This gives the same answer as it would if we computed the trace
1693 from the unembedded (original) matrices::
1695 sage: set_random_seed()
1696 sage: J = RealSymmetricEJA.random_instance()
1697 sage: x,y = J.random_elements(2)
1698 sage: Xe = x.to_matrix()
1699 sage: Ye = y.to_matrix()
1700 sage: X = J.real_unembed(Xe)
1701 sage: Y = J.real_unembed(Ye)
1702 sage: expected = (X*Y).trace()
1703 sage: actual = J.trace_inner_product(Xe,Ye)
1704 sage: actual == expected
1709 sage: set_random_seed()
1710 sage: J = ComplexHermitianEJA.random_instance()
1711 sage: x,y = J.random_elements(2)
1712 sage: Xe = x.to_matrix()
1713 sage: Ye = y.to_matrix()
1714 sage: X = J.real_unembed(Xe)
1715 sage: Y = J.real_unembed(Ye)
1716 sage: expected = (X*Y).trace().real()
1717 sage: actual = J.trace_inner_product(Xe,Ye)
1718 sage: actual == expected
1723 sage: set_random_seed()
1724 sage: J = QuaternionHermitianEJA.random_instance()
1725 sage: x,y = J.random_elements(2)
1726 sage: Xe = x.to_matrix()
1727 sage: Ye = y.to_matrix()
1728 sage: X = J.real_unembed(Xe)
1729 sage: Y = J.real_unembed(Ye)
1730 sage: expected = (X*Y).trace().coefficient_tuple()[0]
1731 sage: actual = J.trace_inner_product(Xe,Ye)
1732 sage: actual == expected
1736 Xu
= cls
.real_unembed(X
)
1737 Yu
= cls
.real_unembed(Y
)
1738 tr
= (Xu
*Yu
).trace()
1741 # Works in QQ, AA, RDF, et cetera.
1743 except AttributeError:
1744 # A quaternion doesn't have a real() method, but does
1745 # have coefficient_tuple() method that returns the
1746 # coefficients of 1, i, j, and k -- in that order.
1747 return tr
.coefficient_tuple()[0]
1750 class RealMatrixEJA(MatrixEJA
):
1752 def dimension_over_reals():
1756 class RealSymmetricEJA(ConcreteEJA
, RealMatrixEJA
):
1758 The rank-n simple EJA consisting of real symmetric n-by-n
1759 matrices, the usual symmetric Jordan product, and the trace inner
1760 product. It has dimension `(n^2 + n)/2` over the reals.
1764 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1768 sage: J = RealSymmetricEJA(2)
1769 sage: e0, e1, e2 = J.gens()
1777 In theory, our "field" can be any subfield of the reals::
1779 sage: RealSymmetricEJA(2, field=RDF)
1780 Euclidean Jordan algebra of dimension 3 over Real Double Field
1781 sage: RealSymmetricEJA(2, field=RR)
1782 Euclidean Jordan algebra of dimension 3 over Real Field with
1783 53 bits of precision
1787 The dimension of this algebra is `(n^2 + n) / 2`::
1789 sage: set_random_seed()
1790 sage: n_max = RealSymmetricEJA._max_random_instance_size()
1791 sage: n = ZZ.random_element(1, n_max)
1792 sage: J = RealSymmetricEJA(n)
1793 sage: J.dimension() == (n^2 + n)/2
1796 The Jordan multiplication is what we think it is::
1798 sage: set_random_seed()
1799 sage: J = RealSymmetricEJA.random_instance()
1800 sage: x,y = J.random_elements(2)
1801 sage: actual = (x*y).to_matrix()
1802 sage: X = x.to_matrix()
1803 sage: Y = y.to_matrix()
1804 sage: expected = (X*Y + Y*X)/2
1805 sage: actual == expected
1807 sage: J(expected) == x*y
1810 We can change the generator prefix::
1812 sage: RealSymmetricEJA(3, prefix='q').gens()
1813 (q0, q1, q2, q3, q4, q5)
1815 We can construct the (trivial) algebra of rank zero::
1817 sage: RealSymmetricEJA(0)
1818 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1822 def _denormalized_basis(cls
, n
):
1824 Return a basis for the space of real symmetric n-by-n matrices.
1828 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1832 sage: set_random_seed()
1833 sage: n = ZZ.random_element(1,5)
1834 sage: B = RealSymmetricEJA._denormalized_basis(n)
1835 sage: all( M.is_symmetric() for M in B)
1839 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1843 for j
in range(i
+1):
1844 Eij
= matrix(ZZ
, n
, lambda k
,l
: k
==i
and l
==j
)
1848 Sij
= Eij
+ Eij
.transpose()
1854 def _max_random_instance_size():
1855 return 4 # Dimension 10
1858 def random_instance(cls
, **kwargs
):
1860 Return a random instance of this type of algebra.
1862 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1863 return cls(n
, **kwargs
)
1865 def __init__(self
, n
, **kwargs
):
1866 # We know this is a valid EJA, but will double-check
1867 # if the user passes check_axioms=True.
1868 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
1870 super(RealSymmetricEJA
, self
).__init
__(self
._denormalized
_basis
(n
),
1871 self
.jordan_product
,
1872 self
.trace_inner_product
,
1875 # TODO: this could be factored out somehow, but is left here
1876 # because the MatrixEJA is not presently a subclass of the
1877 # FDEJA class that defines rank() and one().
1878 self
.rank
.set_cache(n
)
1879 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
1880 self
.one
.set_cache(self(idV
))
1884 class ComplexMatrixEJA(MatrixEJA
):
1885 # A manual dictionary-cache for the complex_extension() method,
1886 # since apparently @classmethods can't also be @cached_methods.
1887 _complex_extension
= {}
1890 def complex_extension(cls
,field
):
1892 The complex field that we embed/unembed, as an extension
1893 of the given ``field``.
1895 if field
in cls
._complex
_extension
:
1896 return cls
._complex
_extension
[field
]
1898 # Sage doesn't know how to adjoin the complex "i" (the root of
1899 # x^2 + 1) to a field in a general way. Here, we just enumerate
1900 # all of the cases that I have cared to support so far.
1902 # Sage doesn't know how to embed AA into QQbar, i.e. how
1903 # to adjoin sqrt(-1) to AA.
1905 elif not field
.is_exact():
1907 F
= field
.complex_field()
1909 # Works for QQ and... maybe some other fields.
1910 R
= PolynomialRing(field
, 'z')
1912 F
= field
.extension(z
**2 + 1, 'I', embedding
=CLF(-1).sqrt())
1914 cls
._complex
_extension
[field
] = F
1918 def dimension_over_reals():
1922 def real_embed(cls
,M
):
1924 Embed the n-by-n complex matrix ``M`` into the space of real
1925 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1926 bi` to the block matrix ``[[a,b],[-b,a]]``.
1930 sage: from mjo.eja.eja_algebra import ComplexMatrixEJA
1934 sage: F = QuadraticField(-1, 'I')
1935 sage: x1 = F(4 - 2*i)
1936 sage: x2 = F(1 + 2*i)
1939 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1940 sage: ComplexMatrixEJA.real_embed(M)
1949 Embedding is a homomorphism (isomorphism, in fact)::
1951 sage: set_random_seed()
1952 sage: n = ZZ.random_element(3)
1953 sage: F = QuadraticField(-1, 'I')
1954 sage: X = random_matrix(F, n)
1955 sage: Y = random_matrix(F, n)
1956 sage: Xe = ComplexMatrixEJA.real_embed(X)
1957 sage: Ye = ComplexMatrixEJA.real_embed(Y)
1958 sage: XYe = ComplexMatrixEJA.real_embed(X*Y)
1963 super(ComplexMatrixEJA
,cls
).real_embed(M
)
1966 # We don't need any adjoined elements...
1967 field
= M
.base_ring().base_ring()
1973 blocks
.append(matrix(field
, 2, [ [ a
, b
],
1976 return matrix
.block(field
, n
, blocks
)
1980 def real_unembed(cls
,M
):
1982 The inverse of _embed_complex_matrix().
1986 sage: from mjo.eja.eja_algebra import ComplexMatrixEJA
1990 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1991 ....: [-2, 1, -4, 3],
1992 ....: [ 9, 10, 11, 12],
1993 ....: [-10, 9, -12, 11] ])
1994 sage: ComplexMatrixEJA.real_unembed(A)
1996 [ 10*I + 9 12*I + 11]
2000 Unembedding is the inverse of embedding::
2002 sage: set_random_seed()
2003 sage: F = QuadraticField(-1, 'I')
2004 sage: M = random_matrix(F, 3)
2005 sage: Me = ComplexMatrixEJA.real_embed(M)
2006 sage: ComplexMatrixEJA.real_unembed(Me) == M
2010 super(ComplexMatrixEJA
,cls
).real_unembed(M
)
2012 d
= cls
.dimension_over_reals()
2013 F
= cls
.complex_extension(M
.base_ring())
2016 # Go top-left to bottom-right (reading order), converting every
2017 # 2-by-2 block we see to a single complex element.
2019 for k
in range(n
/d
):
2020 for j
in range(n
/d
):
2021 submat
= M
[d
*k
:d
*k
+d
,d
*j
:d
*j
+d
]
2022 if submat
[0,0] != submat
[1,1]:
2023 raise ValueError('bad on-diagonal submatrix')
2024 if submat
[0,1] != -submat
[1,0]:
2025 raise ValueError('bad off-diagonal submatrix')
2026 z
= submat
[0,0] + submat
[0,1]*i
2029 return matrix(F
, n
/d
, elements
)
2032 class ComplexHermitianEJA(ConcreteEJA
, ComplexMatrixEJA
):
2034 The rank-n simple EJA consisting of complex Hermitian n-by-n
2035 matrices over the real numbers, the usual symmetric Jordan product,
2036 and the real-part-of-trace inner product. It has dimension `n^2` over
2041 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
2045 In theory, our "field" can be any subfield of the reals::
2047 sage: ComplexHermitianEJA(2, field=RDF)
2048 Euclidean Jordan algebra of dimension 4 over Real Double Field
2049 sage: ComplexHermitianEJA(2, field=RR)
2050 Euclidean Jordan algebra of dimension 4 over Real Field with
2051 53 bits of precision
2055 The dimension of this algebra is `n^2`::
2057 sage: set_random_seed()
2058 sage: n_max = ComplexHermitianEJA._max_random_instance_size()
2059 sage: n = ZZ.random_element(1, n_max)
2060 sage: J = ComplexHermitianEJA(n)
2061 sage: J.dimension() == n^2
2064 The Jordan multiplication is what we think it is::
2066 sage: set_random_seed()
2067 sage: J = ComplexHermitianEJA.random_instance()
2068 sage: x,y = J.random_elements(2)
2069 sage: actual = (x*y).to_matrix()
2070 sage: X = x.to_matrix()
2071 sage: Y = y.to_matrix()
2072 sage: expected = (X*Y + Y*X)/2
2073 sage: actual == expected
2075 sage: J(expected) == x*y
2078 We can change the generator prefix::
2080 sage: ComplexHermitianEJA(2, prefix='z').gens()
2083 We can construct the (trivial) algebra of rank zero::
2085 sage: ComplexHermitianEJA(0)
2086 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2091 def _denormalized_basis(cls
, n
):
2093 Returns a basis for the space of complex Hermitian n-by-n matrices.
2095 Why do we embed these? Basically, because all of numerical linear
2096 algebra assumes that you're working with vectors consisting of `n`
2097 entries from a field and scalars from the same field. There's no way
2098 to tell SageMath that (for example) the vectors contain complex
2099 numbers, while the scalar field is real.
2103 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
2107 sage: set_random_seed()
2108 sage: n = ZZ.random_element(1,5)
2109 sage: B = ComplexHermitianEJA._denormalized_basis(n)
2110 sage: all( M.is_symmetric() for M in B)
2115 R
= PolynomialRing(field
, 'z')
2117 F
= field
.extension(z
**2 + 1, 'I')
2120 # This is like the symmetric case, but we need to be careful:
2122 # * We want conjugate-symmetry, not just symmetry.
2123 # * The diagonal will (as a result) be real.
2126 Eij
= matrix
.zero(F
,n
)
2128 for j
in range(i
+1):
2132 Sij
= cls
.real_embed(Eij
)
2135 # The second one has a minus because it's conjugated.
2136 Eij
[j
,i
] = 1 # Eij = Eij + Eij.transpose()
2137 Sij_real
= cls
.real_embed(Eij
)
2139 # Eij = I*Eij - I*Eij.transpose()
2142 Sij_imag
= cls
.real_embed(Eij
)
2148 # Since we embedded these, we can drop back to the "field" that we
2149 # started with instead of the complex extension "F".
2150 return tuple( s
.change_ring(field
) for s
in S
)
2153 def __init__(self
, n
, **kwargs
):
2154 # We know this is a valid EJA, but will double-check
2155 # if the user passes check_axioms=True.
2156 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2158 super(ComplexHermitianEJA
, self
).__init
__(self
._denormalized
_basis
(n
),
2159 self
.jordan_product
,
2160 self
.trace_inner_product
,
2162 # TODO: this could be factored out somehow, but is left here
2163 # because the MatrixEJA is not presently a subclass of the
2164 # FDEJA class that defines rank() and one().
2165 self
.rank
.set_cache(n
)
2166 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
2167 self
.one
.set_cache(self(idV
))
2170 def _max_random_instance_size():
2171 return 3 # Dimension 9
2174 def random_instance(cls
, **kwargs
):
2176 Return a random instance of this type of algebra.
2178 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2179 return cls(n
, **kwargs
)
2181 class QuaternionMatrixEJA(MatrixEJA
):
2183 # A manual dictionary-cache for the quaternion_extension() method,
2184 # since apparently @classmethods can't also be @cached_methods.
2185 _quaternion_extension
= {}
2188 def quaternion_extension(cls
,field
):
2190 The quaternion field that we embed/unembed, as an extension
2191 of the given ``field``.
2193 if field
in cls
._quaternion
_extension
:
2194 return cls
._quaternion
_extension
[field
]
2196 Q
= QuaternionAlgebra(field
,-1,-1)
2198 cls
._quaternion
_extension
[field
] = Q
2202 def dimension_over_reals():
2206 def real_embed(cls
,M
):
2208 Embed the n-by-n quaternion matrix ``M`` into the space of real
2209 matrices of size 4n-by-4n by first sending each quaternion entry `z
2210 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
2211 c+di],[-c + di, a-bi]]`, and then embedding those into a real
2216 sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA
2220 sage: Q = QuaternionAlgebra(QQ,-1,-1)
2221 sage: i,j,k = Q.gens()
2222 sage: x = 1 + 2*i + 3*j + 4*k
2223 sage: M = matrix(Q, 1, [[x]])
2224 sage: QuaternionMatrixEJA.real_embed(M)
2230 Embedding is a homomorphism (isomorphism, in fact)::
2232 sage: set_random_seed()
2233 sage: n = ZZ.random_element(2)
2234 sage: Q = QuaternionAlgebra(QQ,-1,-1)
2235 sage: X = random_matrix(Q, n)
2236 sage: Y = random_matrix(Q, n)
2237 sage: Xe = QuaternionMatrixEJA.real_embed(X)
2238 sage: Ye = QuaternionMatrixEJA.real_embed(Y)
2239 sage: XYe = QuaternionMatrixEJA.real_embed(X*Y)
2244 super(QuaternionMatrixEJA
,cls
).real_embed(M
)
2245 quaternions
= M
.base_ring()
2248 F
= QuadraticField(-1, 'I')
2253 t
= z
.coefficient_tuple()
2258 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
2259 [-c
+ d
*i
, a
- b
*i
]])
2260 realM
= ComplexMatrixEJA
.real_embed(cplxM
)
2261 blocks
.append(realM
)
2263 # We should have real entries by now, so use the realest field
2264 # we've got for the return value.
2265 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
2270 def real_unembed(cls
,M
):
2272 The inverse of _embed_quaternion_matrix().
2276 sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA
2280 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
2281 ....: [-2, 1, -4, 3],
2282 ....: [-3, 4, 1, -2],
2283 ....: [-4, -3, 2, 1]])
2284 sage: QuaternionMatrixEJA.real_unembed(M)
2285 [1 + 2*i + 3*j + 4*k]
2289 Unembedding is the inverse of embedding::
2291 sage: set_random_seed()
2292 sage: Q = QuaternionAlgebra(QQ, -1, -1)
2293 sage: M = random_matrix(Q, 3)
2294 sage: Me = QuaternionMatrixEJA.real_embed(M)
2295 sage: QuaternionMatrixEJA.real_unembed(Me) == M
2299 super(QuaternionMatrixEJA
,cls
).real_unembed(M
)
2301 d
= cls
.dimension_over_reals()
2303 # Use the base ring of the matrix to ensure that its entries can be
2304 # multiplied by elements of the quaternion algebra.
2305 Q
= cls
.quaternion_extension(M
.base_ring())
2308 # Go top-left to bottom-right (reading order), converting every
2309 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
2312 for l
in range(n
/d
):
2313 for m
in range(n
/d
):
2314 submat
= ComplexMatrixEJA
.real_unembed(
2315 M
[d
*l
:d
*l
+d
,d
*m
:d
*m
+d
] )
2316 if submat
[0,0] != submat
[1,1].conjugate():
2317 raise ValueError('bad on-diagonal submatrix')
2318 if submat
[0,1] != -submat
[1,0].conjugate():
2319 raise ValueError('bad off-diagonal submatrix')
2320 z
= submat
[0,0].real()
2321 z
+= submat
[0,0].imag()*i
2322 z
+= submat
[0,1].real()*j
2323 z
+= submat
[0,1].imag()*k
2326 return matrix(Q
, n
/d
, elements
)
2329 class QuaternionHermitianEJA(ConcreteEJA
, QuaternionMatrixEJA
):
2331 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
2332 matrices, the usual symmetric Jordan product, and the
2333 real-part-of-trace inner product. It has dimension `2n^2 - n` over
2338 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2342 In theory, our "field" can be any subfield of the reals::
2344 sage: QuaternionHermitianEJA(2, field=RDF)
2345 Euclidean Jordan algebra of dimension 6 over Real Double Field
2346 sage: QuaternionHermitianEJA(2, field=RR)
2347 Euclidean Jordan algebra of dimension 6 over Real Field with
2348 53 bits of precision
2352 The dimension of this algebra is `2*n^2 - n`::
2354 sage: set_random_seed()
2355 sage: n_max = QuaternionHermitianEJA._max_random_instance_size()
2356 sage: n = ZZ.random_element(1, n_max)
2357 sage: J = QuaternionHermitianEJA(n)
2358 sage: J.dimension() == 2*(n^2) - n
2361 The Jordan multiplication is what we think it is::
2363 sage: set_random_seed()
2364 sage: J = QuaternionHermitianEJA.random_instance()
2365 sage: x,y = J.random_elements(2)
2366 sage: actual = (x*y).to_matrix()
2367 sage: X = x.to_matrix()
2368 sage: Y = y.to_matrix()
2369 sage: expected = (X*Y + Y*X)/2
2370 sage: actual == expected
2372 sage: J(expected) == x*y
2375 We can change the generator prefix::
2377 sage: QuaternionHermitianEJA(2, prefix='a').gens()
2378 (a0, a1, a2, a3, a4, a5)
2380 We can construct the (trivial) algebra of rank zero::
2382 sage: QuaternionHermitianEJA(0)
2383 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2387 def _denormalized_basis(cls
, n
):
2389 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
2391 Why do we embed these? Basically, because all of numerical
2392 linear algebra assumes that you're working with vectors consisting
2393 of `n` entries from a field and scalars from the same field. There's
2394 no way to tell SageMath that (for example) the vectors contain
2395 complex numbers, while the scalar field is real.
2399 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2403 sage: set_random_seed()
2404 sage: n = ZZ.random_element(1,5)
2405 sage: B = QuaternionHermitianEJA._denormalized_basis(n)
2406 sage: all( M.is_symmetric() for M in B )
2411 Q
= QuaternionAlgebra(QQ
,-1,-1)
2414 # This is like the symmetric case, but we need to be careful:
2416 # * We want conjugate-symmetry, not just symmetry.
2417 # * The diagonal will (as a result) be real.
2420 Eij
= matrix
.zero(Q
,n
)
2422 for j
in range(i
+1):
2426 Sij
= cls
.real_embed(Eij
)
2429 # The second, third, and fourth ones have a minus
2430 # because they're conjugated.
2431 # Eij = Eij + Eij.transpose()
2433 Sij_real
= cls
.real_embed(Eij
)
2435 # Eij = I*(Eij - Eij.transpose())
2438 Sij_I
= cls
.real_embed(Eij
)
2440 # Eij = J*(Eij - Eij.transpose())
2443 Sij_J
= cls
.real_embed(Eij
)
2445 # Eij = K*(Eij - Eij.transpose())
2448 Sij_K
= cls
.real_embed(Eij
)
2454 # Since we embedded these, we can drop back to the "field" that we
2455 # started with instead of the quaternion algebra "Q".
2456 return tuple( s
.change_ring(field
) for s
in S
)
2459 def __init__(self
, n
, **kwargs
):
2460 # We know this is a valid EJA, but will double-check
2461 # if the user passes check_axioms=True.
2462 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2464 super(QuaternionHermitianEJA
, self
).__init
__(self
._denormalized
_basis
(n
),
2465 self
.jordan_product
,
2466 self
.trace_inner_product
,
2468 # TODO: this could be factored out somehow, but is left here
2469 # because the MatrixEJA is not presently a subclass of the
2470 # FDEJA class that defines rank() and one().
2471 self
.rank
.set_cache(n
)
2472 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
2473 self
.one
.set_cache(self(idV
))
2477 def _max_random_instance_size():
2479 The maximum rank of a random QuaternionHermitianEJA.
2481 return 2 # Dimension 6
2484 def random_instance(cls
, **kwargs
):
2486 Return a random instance of this type of algebra.
2488 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2489 return cls(n
, **kwargs
)
2492 class HadamardEJA(ConcreteEJA
):
2494 Return the Euclidean Jordan Algebra corresponding to the set
2495 `R^n` under the Hadamard product.
2497 Note: this is nothing more than the Cartesian product of ``n``
2498 copies of the spin algebra. Once Cartesian product algebras
2499 are implemented, this can go.
2503 sage: from mjo.eja.eja_algebra import HadamardEJA
2507 This multiplication table can be verified by hand::
2509 sage: J = HadamardEJA(3)
2510 sage: e0,e1,e2 = J.gens()
2526 We can change the generator prefix::
2528 sage: HadamardEJA(3, prefix='r').gens()
2532 def __init__(self
, n
, **kwargs
):
2534 jordan_product
= lambda x
,y
: x
2535 inner_product
= lambda x
,y
: x
2537 def jordan_product(x
,y
):
2539 return P( xi
*yi
for (xi
,yi
) in zip(x
,y
) )
2541 def inner_product(x
,y
):
2544 # New defaults for keyword arguments. Don't orthonormalize
2545 # because our basis is already orthonormal with respect to our
2546 # inner-product. Don't check the axioms, because we know this
2547 # is a valid EJA... but do double-check if the user passes
2548 # check_axioms=True. Note: we DON'T override the "check_field"
2549 # default here, because the user can pass in a field!
2550 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2551 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2553 column_basis
= tuple( b
.column() for b
in FreeModule(ZZ
, n
).basis() )
2554 super().__init
__(column_basis
,
2559 self
.rank
.set_cache(n
)
2562 self
.one
.set_cache( self
.zero() )
2564 self
.one
.set_cache( sum(self
.gens()) )
2567 def _max_random_instance_size():
2569 The maximum dimension of a random HadamardEJA.
2574 def random_instance(cls
, **kwargs
):
2576 Return a random instance of this type of algebra.
2578 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2579 return cls(n
, **kwargs
)
2582 class BilinearFormEJA(ConcreteEJA
):
2584 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2585 with the half-trace inner product and jordan product ``x*y =
2586 (<Bx,y>,y_bar>, x0*y_bar + y0*x_bar)`` where `B = 1 \times B22` is
2587 a symmetric positive-definite "bilinear form" matrix. Its
2588 dimension is the size of `B`, and it has rank two in dimensions
2589 larger than two. It reduces to the ``JordanSpinEJA`` when `B` is
2590 the identity matrix of order ``n``.
2592 We insist that the one-by-one upper-left identity block of `B` be
2593 passed in as well so that we can be passed a matrix of size zero
2594 to construct a trivial algebra.
2598 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
2599 ....: JordanSpinEJA)
2603 When no bilinear form is specified, the identity matrix is used,
2604 and the resulting algebra is the Jordan spin algebra::
2606 sage: B = matrix.identity(AA,3)
2607 sage: J0 = BilinearFormEJA(B)
2608 sage: J1 = JordanSpinEJA(3)
2609 sage: J0.multiplication_table() == J0.multiplication_table()
2612 An error is raised if the matrix `B` does not correspond to a
2613 positive-definite bilinear form::
2615 sage: B = matrix.random(QQ,2,3)
2616 sage: J = BilinearFormEJA(B)
2617 Traceback (most recent call last):
2619 ValueError: bilinear form is not positive-definite
2620 sage: B = matrix.zero(QQ,3)
2621 sage: J = BilinearFormEJA(B)
2622 Traceback (most recent call last):
2624 ValueError: bilinear form is not positive-definite
2628 We can create a zero-dimensional algebra::
2630 sage: B = matrix.identity(AA,0)
2631 sage: J = BilinearFormEJA(B)
2635 We can check the multiplication condition given in the Jordan, von
2636 Neumann, and Wigner paper (and also discussed on my "On the
2637 symmetry..." paper). Note that this relies heavily on the standard
2638 choice of basis, as does anything utilizing the bilinear form
2639 matrix. We opt not to orthonormalize the basis, because if we
2640 did, we would have to normalize the `s_{i}` in a similar manner::
2642 sage: set_random_seed()
2643 sage: n = ZZ.random_element(5)
2644 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2645 sage: B11 = matrix.identity(QQ,1)
2646 sage: B22 = M.transpose()*M
2647 sage: B = block_matrix(2,2,[ [B11,0 ],
2649 sage: J = BilinearFormEJA(B, orthonormalize=False)
2650 sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
2651 sage: V = J.vector_space()
2652 sage: sis = [ J( V([0] + (M.inverse()*ei).list()).column() )
2653 ....: for ei in eis ]
2654 sage: actual = [ sis[i]*sis[j]
2655 ....: for i in range(n-1)
2656 ....: for j in range(n-1) ]
2657 sage: expected = [ J.one() if i == j else J.zero()
2658 ....: for i in range(n-1)
2659 ....: for j in range(n-1) ]
2660 sage: actual == expected
2664 def __init__(self
, B
, **kwargs
):
2665 # The matrix "B" is supplied by the user in most cases,
2666 # so it makes sense to check whether or not its positive-
2667 # definite unless we are specifically asked not to...
2668 if ("check_axioms" not in kwargs
) or kwargs
["check_axioms"]:
2669 if not B
.is_positive_definite():
2670 raise ValueError("bilinear form is not positive-definite")
2672 # However, all of the other data for this EJA is computed
2673 # by us in manner that guarantees the axioms are
2674 # satisfied. So, again, unless we are specifically asked to
2675 # verify things, we'll skip the rest of the checks.
2676 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2678 def inner_product(x
,y
):
2679 return (y
.T
*B
*x
)[0,0]
2681 def jordan_product(x
,y
):
2687 z0
= inner_product(y
,x
)
2688 zbar
= y0
*xbar
+ x0
*ybar
2689 return P([z0
] + zbar
.list())
2692 column_basis
= tuple( b
.column() for b
in FreeModule(ZZ
, n
).basis() )
2693 super(BilinearFormEJA
, self
).__init
__(column_basis
,
2698 # The rank of this algebra is two, unless we're in a
2699 # one-dimensional ambient space (because the rank is bounded
2700 # by the ambient dimension).
2701 self
.rank
.set_cache(min(n
,2))
2704 self
.one
.set_cache( self
.zero() )
2706 self
.one
.set_cache( self
.monomial(0) )
2709 def _max_random_instance_size():
2711 The maximum dimension of a random BilinearFormEJA.
2716 def random_instance(cls
, **kwargs
):
2718 Return a random instance of this algebra.
2720 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2722 B
= matrix
.identity(ZZ
, n
)
2723 return cls(B
, **kwargs
)
2725 B11
= matrix
.identity(ZZ
, 1)
2726 M
= matrix
.random(ZZ
, n
-1)
2727 I
= matrix
.identity(ZZ
, n
-1)
2729 while alpha
.is_zero():
2730 alpha
= ZZ
.random_element().abs()
2731 B22
= M
.transpose()*M
+ alpha
*I
2733 from sage
.matrix
.special
import block_matrix
2734 B
= block_matrix(2,2, [ [B11
, ZZ(0) ],
2737 return cls(B
, **kwargs
)
2740 class JordanSpinEJA(BilinearFormEJA
):
2742 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2743 with the usual inner product and jordan product ``x*y =
2744 (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2749 sage: from mjo.eja.eja_algebra import JordanSpinEJA
2753 This multiplication table can be verified by hand::
2755 sage: J = JordanSpinEJA(4)
2756 sage: e0,e1,e2,e3 = J.gens()
2772 We can change the generator prefix::
2774 sage: JordanSpinEJA(2, prefix='B').gens()
2779 Ensure that we have the usual inner product on `R^n`::
2781 sage: set_random_seed()
2782 sage: J = JordanSpinEJA.random_instance()
2783 sage: x,y = J.random_elements(2)
2784 sage: actual = x.inner_product(y)
2785 sage: expected = x.to_vector().inner_product(y.to_vector())
2786 sage: actual == expected
2790 def __init__(self
, n
, **kwargs
):
2791 # This is a special case of the BilinearFormEJA with the
2792 # identity matrix as its bilinear form.
2793 B
= matrix
.identity(ZZ
, n
)
2795 # Don't orthonormalize because our basis is already
2796 # orthonormal with respect to our inner-product.
2797 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2799 # But also don't pass check_field=False here, because the user
2800 # can pass in a field!
2801 super(JordanSpinEJA
, self
).__init
__(B
, **kwargs
)
2804 def _max_random_instance_size():
2806 The maximum dimension of a random JordanSpinEJA.
2811 def random_instance(cls
, **kwargs
):
2813 Return a random instance of this type of algebra.
2815 Needed here to override the implementation for ``BilinearFormEJA``.
2817 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2818 return cls(n
, **kwargs
)
2821 class TrivialEJA(ConcreteEJA
):
2823 The trivial Euclidean Jordan algebra consisting of only a zero element.
2827 sage: from mjo.eja.eja_algebra import TrivialEJA
2831 sage: J = TrivialEJA()
2838 sage: 7*J.one()*12*J.one()
2840 sage: J.one().inner_product(J.one())
2842 sage: J.one().norm()
2844 sage: J.one().subalgebra_generated_by()
2845 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2850 def __init__(self
, **kwargs
):
2851 jordan_product
= lambda x
,y
: x
2852 inner_product
= lambda x
,y
: 0
2855 # New defaults for keyword arguments
2856 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2857 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2859 super(TrivialEJA
, self
).__init
__(basis
,
2863 # The rank is zero using my definition, namely the dimension of the
2864 # largest subalgebra generated by any element.
2865 self
.rank
.set_cache(0)
2866 self
.one
.set_cache( self
.zero() )
2869 def random_instance(cls
, **kwargs
):
2870 # We don't take a "size" argument so the superclass method is
2871 # inappropriate for us.
2872 return cls(**kwargs
)
2875 class CartesianProductEJA(CombinatorialFreeModule_CartesianProduct
,
2876 FiniteDimensionalEJA
):
2878 The external (orthogonal) direct sum of two or more Euclidean
2879 Jordan algebras. Every Euclidean Jordan algebra decomposes into an
2880 orthogonal direct sum of simple Euclidean Jordan algebras which is
2881 then isometric to a Cartesian product, so no generality is lost by
2882 providing only this construction.
2886 sage: from mjo.eja.eja_algebra import (random_eja,
2887 ....: CartesianProductEJA,
2889 ....: JordanSpinEJA,
2890 ....: RealSymmetricEJA)
2894 The Jordan product is inherited from our factors and implemented by
2895 our CombinatorialFreeModule Cartesian product superclass::
2897 sage: set_random_seed()
2898 sage: J1 = HadamardEJA(2)
2899 sage: J2 = RealSymmetricEJA(2)
2900 sage: J = cartesian_product([J1,J2])
2901 sage: x,y = J.random_elements(2)
2905 The ability to retrieve the original factors is implemented by our
2906 CombinatorialFreeModule Cartesian product superclass::
2908 sage: J1 = HadamardEJA(2, field=QQ)
2909 sage: J2 = JordanSpinEJA(3, field=QQ)
2910 sage: J = cartesian_product([J1,J2])
2911 sage: J.cartesian_factors()
2912 (Euclidean Jordan algebra of dimension 2 over Rational Field,
2913 Euclidean Jordan algebra of dimension 3 over Rational Field)
2915 You can provide more than two factors::
2917 sage: J1 = HadamardEJA(2)
2918 sage: J2 = JordanSpinEJA(3)
2919 sage: J3 = RealSymmetricEJA(3)
2920 sage: cartesian_product([J1,J2,J3])
2921 Euclidean Jordan algebra of dimension 2 over Algebraic Real
2922 Field (+) Euclidean Jordan algebra of dimension 3 over Algebraic
2923 Real Field (+) Euclidean Jordan algebra of dimension 6 over
2924 Algebraic Real Field
2926 Rank is additive on a Cartesian product::
2928 sage: J1 = HadamardEJA(1)
2929 sage: J2 = RealSymmetricEJA(2)
2930 sage: J = cartesian_product([J1,J2])
2931 sage: J1.rank.clear_cache()
2932 sage: J2.rank.clear_cache()
2933 sage: J.rank.clear_cache()
2936 sage: J.rank() == J1.rank() + J2.rank()
2939 The same rank computation works over the rationals, with whatever
2942 sage: J1 = HadamardEJA(1, field=QQ, orthonormalize=False)
2943 sage: J2 = RealSymmetricEJA(2, field=QQ, orthonormalize=False)
2944 sage: J = cartesian_product([J1,J2])
2945 sage: J1.rank.clear_cache()
2946 sage: J2.rank.clear_cache()
2947 sage: J.rank.clear_cache()
2950 sage: J.rank() == J1.rank() + J2.rank()
2953 The product algebra will be associative if and only if all of its
2954 components are associative::
2956 sage: J1 = HadamardEJA(2)
2957 sage: J1.is_associative()
2959 sage: J2 = HadamardEJA(3)
2960 sage: J2.is_associative()
2962 sage: J3 = RealSymmetricEJA(3)
2963 sage: J3.is_associative()
2965 sage: CP1 = cartesian_product([J1,J2])
2966 sage: CP1.is_associative()
2968 sage: CP2 = cartesian_product([J1,J3])
2969 sage: CP2.is_associative()
2974 All factors must share the same base field::
2976 sage: J1 = HadamardEJA(2, field=QQ)
2977 sage: J2 = RealSymmetricEJA(2)
2978 sage: CartesianProductEJA((J1,J2))
2979 Traceback (most recent call last):
2981 ValueError: all factors must share the same base field
2983 The cached unit element is the same one that would be computed::
2985 sage: set_random_seed() # long time
2986 sage: J1 = random_eja() # long time
2987 sage: J2 = random_eja() # long time
2988 sage: J = cartesian_product([J1,J2]) # long time
2989 sage: actual = J.one() # long time
2990 sage: J.one.clear_cache() # long time
2991 sage: expected = J.one() # long time
2992 sage: actual == expected # long time
2996 Element
= FiniteDimensionalEJAElement
2999 def __init__(self
, algebras
, **kwargs
):
3000 CombinatorialFreeModule_CartesianProduct
.__init
__(self
,
3003 field
= algebras
[0].base_ring()
3004 if not all( J
.base_ring() == field
for J
in algebras
):
3005 raise ValueError("all factors must share the same base field")
3007 associative
= all( m
.is_associative() for m
in algebras
)
3009 # The definition of matrix_space() and self.basis() relies
3010 # only on the stuff in the CFM_CartesianProduct class, which
3011 # we've already initialized.
3012 Js
= self
.cartesian_factors()
3014 MS
= self
.matrix_space()
3016 MS(tuple( self
.cartesian_projection(i
)(b
).to_matrix()
3017 for i
in range(m
) ))
3018 for b
in self
.basis()
3021 # Define jordan/inner products that operate on that matrix_basis.
3022 def jordan_product(x
,y
):
3024 (Js
[i
](x
[i
])*Js
[i
](y
[i
])).to_matrix() for i
in range(m
)
3027 def inner_product(x
, y
):
3029 Js
[i
](x
[i
]).inner_product(Js
[i
](y
[i
])) for i
in range(m
)
3032 # There's no need to check the field since it already came
3033 # from an EJA. Likewise the axioms are guaranteed to be
3034 # satisfied, unless the guy writing this class sucks.
3036 # If you want the basis to be orthonormalized, orthonormalize
3038 FiniteDimensionalEJA
.__init
__(self
,
3043 orthonormalize
=False,
3044 associative
=associative
,
3045 cartesian_product
=True,
3049 ones
= tuple(J
.one() for J
in algebras
)
3050 self
.one
.set_cache(self
._cartesian
_product
_of
_elements
(ones
))
3051 self
.rank
.set_cache(sum(J
.rank() for J
in algebras
))
3053 def matrix_space(self
):
3055 Return the space that our matrix basis lives in as a Cartesian
3060 sage: from mjo.eja.eja_algebra import (HadamardEJA,
3061 ....: RealSymmetricEJA)
3065 sage: J1 = HadamardEJA(1)
3066 sage: J2 = RealSymmetricEJA(2)
3067 sage: J = cartesian_product([J1,J2])
3068 sage: J.matrix_space()
3069 The Cartesian product of (Full MatrixSpace of 1 by 1 dense
3070 matrices over Algebraic Real Field, Full MatrixSpace of 2
3071 by 2 dense matrices over Algebraic Real Field)
3074 from sage
.categories
.cartesian_product
import cartesian_product
3075 return cartesian_product( [J
.matrix_space()
3076 for J
in self
.cartesian_factors()] )
3079 def cartesian_projection(self
, i
):
3083 sage: from mjo.eja.eja_algebra import (random_eja,
3084 ....: JordanSpinEJA,
3086 ....: RealSymmetricEJA,
3087 ....: ComplexHermitianEJA)
3091 The projection morphisms are Euclidean Jordan algebra
3094 sage: J1 = HadamardEJA(2)
3095 sage: J2 = RealSymmetricEJA(2)
3096 sage: J = cartesian_product([J1,J2])
3097 sage: J.cartesian_projection(0)
3098 Linear operator between finite-dimensional Euclidean Jordan
3099 algebras represented by the matrix:
3102 Domain: Euclidean Jordan algebra of dimension 2 over Algebraic
3103 Real Field (+) Euclidean Jordan algebra of dimension 3 over
3104 Algebraic Real Field
3105 Codomain: Euclidean Jordan algebra of dimension 2 over Algebraic
3107 sage: J.cartesian_projection(1)
3108 Linear operator between finite-dimensional Euclidean Jordan
3109 algebras represented by the matrix:
3113 Domain: Euclidean Jordan algebra of dimension 2 over Algebraic
3114 Real Field (+) Euclidean Jordan algebra of dimension 3 over
3115 Algebraic Real Field
3116 Codomain: Euclidean Jordan algebra of dimension 3 over Algebraic
3119 The projections work the way you'd expect on the vector
3120 representation of an element::
3122 sage: J1 = JordanSpinEJA(2)
3123 sage: J2 = ComplexHermitianEJA(2)
3124 sage: J = cartesian_product([J1,J2])
3125 sage: pi_left = J.cartesian_projection(0)
3126 sage: pi_right = J.cartesian_projection(1)
3127 sage: pi_left(J.one()).to_vector()
3129 sage: pi_right(J.one()).to_vector()
3131 sage: J.one().to_vector()
3136 The answer never changes::
3138 sage: set_random_seed()
3139 sage: J1 = random_eja()
3140 sage: J2 = random_eja()
3141 sage: J = cartesian_product([J1,J2])
3142 sage: P0 = J.cartesian_projection(0)
3143 sage: P1 = J.cartesian_projection(0)
3148 Ji
= self
.cartesian_factors()[i
]
3149 # Requires the fix on Trac 31421/31422 to work!
3150 Pi
= super().cartesian_projection(i
)
3151 return FiniteDimensionalEJAOperator(self
,Ji
,Pi
.matrix())
3154 def cartesian_embedding(self
, i
):
3158 sage: from mjo.eja.eja_algebra import (random_eja,
3159 ....: JordanSpinEJA,
3161 ....: RealSymmetricEJA)
3165 The embedding morphisms are Euclidean Jordan algebra
3168 sage: J1 = HadamardEJA(2)
3169 sage: J2 = RealSymmetricEJA(2)
3170 sage: J = cartesian_product([J1,J2])
3171 sage: J.cartesian_embedding(0)
3172 Linear operator between finite-dimensional Euclidean Jordan
3173 algebras represented by the matrix:
3179 Domain: Euclidean Jordan algebra of dimension 2 over
3180 Algebraic Real Field
3181 Codomain: Euclidean Jordan algebra of dimension 2 over
3182 Algebraic Real Field (+) Euclidean Jordan algebra of
3183 dimension 3 over Algebraic Real Field
3184 sage: J.cartesian_embedding(1)
3185 Linear operator between finite-dimensional Euclidean Jordan
3186 algebras represented by the matrix:
3192 Domain: Euclidean Jordan algebra of dimension 3 over
3193 Algebraic Real Field
3194 Codomain: Euclidean Jordan algebra of dimension 2 over
3195 Algebraic Real Field (+) Euclidean Jordan algebra of
3196 dimension 3 over Algebraic Real Field
3198 The embeddings work the way you'd expect on the vector
3199 representation of an element::
3201 sage: J1 = JordanSpinEJA(3)
3202 sage: J2 = RealSymmetricEJA(2)
3203 sage: J = cartesian_product([J1,J2])
3204 sage: iota_left = J.cartesian_embedding(0)
3205 sage: iota_right = J.cartesian_embedding(1)
3206 sage: iota_left(J1.zero()) == J.zero()
3208 sage: iota_right(J2.zero()) == J.zero()
3210 sage: J1.one().to_vector()
3212 sage: iota_left(J1.one()).to_vector()
3214 sage: J2.one().to_vector()
3216 sage: iota_right(J2.one()).to_vector()
3218 sage: J.one().to_vector()
3223 The answer never changes::
3225 sage: set_random_seed()
3226 sage: J1 = random_eja()
3227 sage: J2 = random_eja()
3228 sage: J = cartesian_product([J1,J2])
3229 sage: E0 = J.cartesian_embedding(0)
3230 sage: E1 = J.cartesian_embedding(0)
3234 Composing a projection with the corresponding inclusion should
3235 produce the identity map, and mismatching them should produce
3238 sage: set_random_seed()
3239 sage: J1 = random_eja()
3240 sage: J2 = random_eja()
3241 sage: J = cartesian_product([J1,J2])
3242 sage: iota_left = J.cartesian_embedding(0)
3243 sage: iota_right = J.cartesian_embedding(1)
3244 sage: pi_left = J.cartesian_projection(0)
3245 sage: pi_right = J.cartesian_projection(1)
3246 sage: pi_left*iota_left == J1.one().operator()
3248 sage: pi_right*iota_right == J2.one().operator()
3250 sage: (pi_left*iota_right).is_zero()
3252 sage: (pi_right*iota_left).is_zero()
3256 Ji
= self
.cartesian_factors()[i
]
3257 # Requires the fix on Trac 31421/31422 to work!
3258 Ei
= super().cartesian_embedding(i
)
3259 return FiniteDimensionalEJAOperator(Ji
,self
,Ei
.matrix())
3263 FiniteDimensionalEJA
.CartesianProduct
= CartesianProductEJA
3265 class RationalBasisCartesianProductEJA(CartesianProductEJA
,
3268 A separate class for products of algebras for which we know a
3273 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
3274 ....: RealSymmetricEJA)
3278 This gives us fast characteristic polynomial computations in
3279 product algebras, too::
3282 sage: J1 = JordanSpinEJA(2)
3283 sage: J2 = RealSymmetricEJA(3)
3284 sage: J = cartesian_product([J1,J2])
3285 sage: J.characteristic_polynomial_of().degree()
3291 def __init__(self
, algebras
, **kwargs
):
3292 CartesianProductEJA
.__init
__(self
, algebras
, **kwargs
)
3294 self
._rational
_algebra
= None
3295 if self
.vector_space().base_field() is not QQ
:
3296 self
._rational
_algebra
= cartesian_product([
3297 r
._rational
_algebra
for r
in algebras
3301 RationalBasisEJA
.CartesianProduct
= RationalBasisCartesianProductEJA
3303 random_eja
= ConcreteEJA
.random_instance