2 Euclidean Jordan Algebras. These are formally-real Jordan Algebras;
3 specifically those where u^2 + v^2 = 0 implies that u = v = 0. They
4 are used in optimization, and have some additional nice methods beyond
5 what can be supported in a general Jordan Algebra.
10 sage: from mjo.eja.eja_algebra import random_eja
15 Euclidean Jordan algebra of dimension...
19 from itertools
import repeat
21 from sage
.algebras
.quatalg
.quaternion_algebra
import QuaternionAlgebra
22 from sage
.categories
.magmatic_algebras
import MagmaticAlgebras
23 from sage
.categories
.sets_cat
import cartesian_product
24 from sage
.combinat
.free_module
import (CombinatorialFreeModule
,
25 CombinatorialFreeModule_CartesianProduct
)
26 from sage
.matrix
.constructor
import matrix
27 from sage
.matrix
.matrix_space
import MatrixSpace
28 from sage
.misc
.cachefunc
import cached_method
29 from sage
.misc
.table
import table
30 from sage
.modules
.free_module
import FreeModule
, VectorSpace
31 from sage
.rings
.all
import (ZZ
, QQ
, AA
, QQbar
, RR
, RLF
, CLF
,
34 from mjo
.eja
.eja_element
import FiniteDimensionalEJAElement
35 from mjo
.eja
.eja_operator
import FiniteDimensionalEJAOperator
36 from mjo
.eja
.eja_utils
import _all2list
, _mat2vec
38 class FiniteDimensionalEJA(CombinatorialFreeModule
):
40 A finite-dimensional Euclidean Jordan algebra.
44 - basis -- a tuple of basis elements in "matrix form," which
45 must be the same form as the arguments to ``jordan_product``
46 and ``inner_product``. In reality, "matrix form" can be either
47 vectors, matrices, or a Cartesian product (ordered tuple)
48 of vectors or matrices. All of these would ideally be vector
49 spaces in sage with no special-casing needed; but in reality
50 we turn vectors into column-matrices and Cartesian products
51 `(a,b)` into column matrices `(a,b)^{T}` after converting
52 `a` and `b` themselves.
54 - jordan_product -- function of two ``basis`` elements (in
55 matrix form) that returns their jordan product, also in matrix
56 form; this will be applied to ``basis`` to compute a
57 multiplication table for the algebra.
59 - inner_product -- function of two ``basis`` elements (in matrix
60 form) that returns their inner product. This will be applied
61 to ``basis`` to compute an inner-product table (basically a
62 matrix) for this algebra.
64 Element
= FiniteDimensionalEJAElement
73 cartesian_product
=False,
78 # Keep track of whether or not the matrix basis consists of
79 # tuples, since we need special cases for them damned near
80 # everywhere. This is INDEPENDENT of whether or not the
81 # algebra is a cartesian product, since a subalgebra of a
82 # cartesian product will have a basis of tuples, but will not
83 # in general itself be a cartesian product algebra.
84 self
._matrix
_basis
_is
_cartesian
= False
87 if hasattr(basis
[0], 'cartesian_factors'):
88 self
._matrix
_basis
_is
_cartesian
= True
91 if not field
.is_subring(RR
):
92 # Note: this does return true for the real algebraic
93 # field, the rationals, and any quadratic field where
94 # we've specified a real embedding.
95 raise ValueError("scalar field is not real")
97 # If the basis given to us wasn't over the field that it's
98 # supposed to be over, fix that. Or, you know, crash.
99 if not cartesian_product
:
100 # The field for a cartesian product algebra comes from one
101 # of its factors and is the same for all factors, so
102 # there's no need to "reapply" it on product algebras.
103 if self
._matrix
_basis
_is
_cartesian
:
104 # OK since if n == 0, the basis does not consist of tuples.
105 P
= basis
[0].parent()
106 basis
= tuple( P(tuple(b_i
.change_ring(field
) for b_i
in b
))
109 basis
= tuple( b
.change_ring(field
) for b
in basis
)
113 # Check commutativity of the Jordan and inner-products.
114 # This has to be done before we build the multiplication
115 # and inner-product tables/matrices, because we take
116 # advantage of symmetry in the process.
117 if not all( jordan_product(bi
,bj
) == jordan_product(bj
,bi
)
120 raise ValueError("Jordan product is not commutative")
122 if not all( inner_product(bi
,bj
) == inner_product(bj
,bi
)
125 raise ValueError("inner-product is not commutative")
128 category
= MagmaticAlgebras(field
).FiniteDimensional()
129 category
= category
.WithBasis().Unital()
131 # Element subalgebras can take advantage of this.
132 category
= category
.Associative()
133 if cartesian_product
:
134 category
= category
.CartesianProducts()
136 # Call the superclass constructor so that we can use its from_vector()
137 # method to build our multiplication table.
138 CombinatorialFreeModule
.__init
__(self
,
145 # Now comes all of the hard work. We'll be constructing an
146 # ambient vector space V that our (vectorized) basis lives in,
147 # as well as a subspace W of V spanned by those (vectorized)
148 # basis elements. The W-coordinates are the coefficients that
149 # we see in things like x = 1*e1 + 2*e2.
154 degree
= len(_all2list(basis
[0]))
156 # Build an ambient space that fits our matrix basis when
157 # written out as "long vectors."
158 V
= VectorSpace(field
, degree
)
160 # The matrix that will hole the orthonormal -> unorthonormal
161 # coordinate transformation.
162 self
._deortho
_matrix
= None
165 # Save a copy of the un-orthonormalized basis for later.
166 # Convert it to ambient V (vector) coordinates while we're
167 # at it, because we'd have to do it later anyway.
168 deortho_vector_basis
= tuple( V(_all2list(b
)) for b
in basis
)
170 from mjo
.eja
.eja_utils
import gram_schmidt
171 basis
= tuple(gram_schmidt(basis
, inner_product
))
173 # Save the (possibly orthonormalized) matrix basis for
175 self
._matrix
_basis
= basis
177 # Now create the vector space for the algebra, which will have
178 # its own set of non-ambient coordinates (in terms of the
180 vector_basis
= tuple( V(_all2list(b
)) for b
in basis
)
181 W
= V
.span_of_basis( vector_basis
, check
=check_axioms
)
184 # Now "W" is the vector space of our algebra coordinates. The
185 # variables "X1", "X2",... refer to the entries of vectors in
186 # W. Thus to convert back and forth between the orthonormal
187 # coordinates and the given ones, we need to stick the original
189 U
= V
.span_of_basis( deortho_vector_basis
, check
=check_axioms
)
190 self
._deortho
_matrix
= matrix( U
.coordinate_vector(q
)
191 for q
in vector_basis
)
194 # Now we actually compute the multiplication and inner-product
195 # tables/matrices using the possibly-orthonormalized basis.
196 self
._inner
_product
_matrix
= matrix
.identity(field
, n
)
197 self
._multiplication
_table
= [ [0 for j
in range(i
+1)]
200 # Note: the Jordan and inner-products are defined in terms
201 # of the ambient basis. It's important that their arguments
202 # are in ambient coordinates as well.
205 # ortho basis w.r.t. ambient coords
209 # The jordan product returns a matrixy answer, so we
210 # have to convert it to the algebra coordinates.
211 elt
= jordan_product(q_i
, q_j
)
212 elt
= W
.coordinate_vector(V(_all2list(elt
)))
213 self
._multiplication
_table
[i
][j
] = self
.from_vector(elt
)
215 if not orthonormalize
:
216 # If we're orthonormalizing the basis with respect
217 # to an inner-product, then the inner-product
218 # matrix with respect to the resulting basis is
219 # just going to be the identity.
220 ip
= inner_product(q_i
, q_j
)
221 self
._inner
_product
_matrix
[i
,j
] = ip
222 self
._inner
_product
_matrix
[j
,i
] = ip
224 self
._inner
_product
_matrix
._cache
= {'hermitian': True}
225 self
._inner
_product
_matrix
.set_immutable()
228 if not self
._is
_jordanian
():
229 raise ValueError("Jordan identity does not hold")
230 if not self
._inner
_product
_is
_associative
():
231 raise ValueError("inner product is not associative")
234 def _coerce_map_from_base_ring(self
):
236 Disable the map from the base ring into the algebra.
238 Performing a nonsense conversion like this automatically
239 is counterpedagogical. The fallback is to try the usual
240 element constructor, which should also fail.
244 sage: from mjo.eja.eja_algebra import random_eja
248 sage: set_random_seed()
249 sage: J = random_eja()
251 Traceback (most recent call last):
253 ValueError: not an element of this algebra
259 def product_on_basis(self
, i
, j
):
261 Returns the Jordan product of the `i` and `j`th basis elements.
263 This completely defines the Jordan product on the algebra, and
264 is used direclty by our superclass machinery to implement
269 sage: from mjo.eja.eja_algebra import random_eja
273 sage: set_random_seed()
274 sage: J = random_eja()
275 sage: n = J.dimension()
278 sage: ei_ej = J.zero()*J.zero()
280 ....: i = ZZ.random_element(n)
281 ....: j = ZZ.random_element(n)
282 ....: ei = J.gens()[i]
283 ....: ej = J.gens()[j]
284 ....: ei_ej = J.product_on_basis(i,j)
289 # We only stored the lower-triangular portion of the
290 # multiplication table.
292 return self
._multiplication
_table
[i
][j
]
294 return self
._multiplication
_table
[j
][i
]
296 def inner_product(self
, x
, y
):
298 The inner product associated with this Euclidean Jordan algebra.
300 Defaults to the trace inner product, but can be overridden by
301 subclasses if they are sure that the necessary properties are
306 sage: from mjo.eja.eja_algebra import (random_eja,
308 ....: BilinearFormEJA)
312 Our inner product is "associative," which means the following for
313 a symmetric bilinear form::
315 sage: set_random_seed()
316 sage: J = random_eja()
317 sage: x,y,z = J.random_elements(3)
318 sage: (x*y).inner_product(z) == y.inner_product(x*z)
323 Ensure that this is the usual inner product for the algebras
326 sage: set_random_seed()
327 sage: J = HadamardEJA.random_instance()
328 sage: x,y = J.random_elements(2)
329 sage: actual = x.inner_product(y)
330 sage: expected = x.to_vector().inner_product(y.to_vector())
331 sage: actual == expected
334 Ensure that this is one-half of the trace inner-product in a
335 BilinearFormEJA that isn't just the reals (when ``n`` isn't
336 one). This is in Faraut and Koranyi, and also my "On the
339 sage: set_random_seed()
340 sage: J = BilinearFormEJA.random_instance()
341 sage: n = J.dimension()
342 sage: x = J.random_element()
343 sage: y = J.random_element()
344 sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
348 B
= self
._inner
_product
_matrix
349 return (B
*x
.to_vector()).inner_product(y
.to_vector())
352 def is_associative(self
):
354 Return whether or not this algebra's Jordan product is associative.
358 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
362 sage: J = ComplexHermitianEJA(3, field=QQ, orthonormalize=False)
363 sage: J.is_associative()
365 sage: x = sum(J.gens())
366 sage: A = x.subalgebra_generated_by(orthonormalize=False)
367 sage: A.is_associative()
371 return "Associative" in self
.category().axioms()
373 def _is_jordanian(self
):
375 Whether or not this algebra's multiplication table respects the
376 Jordan identity `(x^{2})(xy) = x(x^{2}y)`.
378 We only check one arrangement of `x` and `y`, so for a
379 ``True`` result to be truly true, you should also check
380 :meth:`is_commutative`. This method should of course always
381 return ``True``, unless this algebra was constructed with
382 ``check_axioms=False`` and passed an invalid multiplication table.
384 return all( (self
.gens()[i
]**2)*(self
.gens()[i
]*self
.gens()[j
])
386 (self
.gens()[i
])*((self
.gens()[i
]**2)*self
.gens()[j
])
387 for i
in range(self
.dimension())
388 for j
in range(self
.dimension()) )
390 def _inner_product_is_associative(self
):
392 Return whether or not this algebra's inner product `B` is
393 associative; that is, whether or not `B(xy,z) = B(x,yz)`.
395 This method should of course always return ``True``, unless
396 this algebra was constructed with ``check_axioms=False`` and
397 passed an invalid Jordan or inner-product.
400 # Used to check whether or not something is zero in an inexact
401 # ring. This number is sufficient to allow the construction of
402 # QuaternionHermitianEJA(2, field=RDF) with check_axioms=True.
405 for i
in range(self
.dimension()):
406 for j
in range(self
.dimension()):
407 for k
in range(self
.dimension()):
411 diff
= (x
*y
).inner_product(z
) - x
.inner_product(y
*z
)
413 if self
.base_ring().is_exact():
417 if diff
.abs() > epsilon
:
422 def _element_constructor_(self
, elt
):
424 Construct an element of this algebra from its vector or matrix
427 This gets called only after the parent element _call_ method
428 fails to find a coercion for the argument.
432 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
434 ....: RealSymmetricEJA)
438 The identity in `S^n` is converted to the identity in the EJA::
440 sage: J = RealSymmetricEJA(3)
441 sage: I = matrix.identity(QQ,3)
442 sage: J(I) == J.one()
445 This skew-symmetric matrix can't be represented in the EJA::
447 sage: J = RealSymmetricEJA(3)
448 sage: A = matrix(QQ,3, lambda i,j: i-j)
450 Traceback (most recent call last):
452 ValueError: not an element of this algebra
454 Tuples work as well, provided that the matrix basis for the
455 algebra consists of them::
457 sage: J1 = HadamardEJA(3)
458 sage: J2 = RealSymmetricEJA(2)
459 sage: J = cartesian_product([J1,J2])
460 sage: J( (J1.matrix_basis()[1], J2.matrix_basis()[2]) )
465 Ensure that we can convert any element of the two non-matrix
466 simple algebras (whose matrix representations are columns)
467 back and forth faithfully::
469 sage: set_random_seed()
470 sage: J = HadamardEJA.random_instance()
471 sage: x = J.random_element()
472 sage: J(x.to_vector().column()) == x
474 sage: J = JordanSpinEJA.random_instance()
475 sage: x = J.random_element()
476 sage: J(x.to_vector().column()) == x
479 We cannot coerce elements between algebras just because their
480 matrix representations are compatible::
482 sage: J1 = HadamardEJA(3)
483 sage: J2 = JordanSpinEJA(3)
485 Traceback (most recent call last):
487 ValueError: not an element of this algebra
489 Traceback (most recent call last):
491 ValueError: not an element of this algebra
494 msg
= "not an element of this algebra"
495 if elt
in self
.base_ring():
496 # Ensure that no base ring -> algebra coercion is performed
497 # by this method. There's some stupidity in sage that would
498 # otherwise propagate to this method; for example, sage thinks
499 # that the integer 3 belongs to the space of 2-by-2 matrices.
500 raise ValueError(msg
)
503 # Try to convert a vector into a column-matrix...
505 except (AttributeError, TypeError):
506 # and ignore failure, because we weren't really expecting
507 # a vector as an argument anyway.
510 if elt
not in self
.matrix_space():
511 raise ValueError(msg
)
513 # Thanks for nothing! Matrix spaces aren't vector spaces in
514 # Sage, so we have to figure out its matrix-basis coordinates
515 # ourselves. We use the basis space's ring instead of the
516 # element's ring because the basis space might be an algebraic
517 # closure whereas the base ring of the 3-by-3 identity matrix
518 # could be QQ instead of QQbar.
520 # And, we also have to handle Cartesian product bases (when
521 # the matrix basis consists of tuples) here. The "good news"
522 # is that we're already converting everything to long vectors,
523 # and that strategy works for tuples as well.
525 # We pass check=False because the matrix basis is "guaranteed"
526 # to be linearly independent... right? Ha ha.
528 V
= VectorSpace(self
.base_ring(), len(elt
))
529 W
= V
.span_of_basis( (V(_all2list(s
)) for s
in self
.matrix_basis()),
533 coords
= W
.coordinate_vector(V(elt
))
534 except ArithmeticError: # vector is not in free module
535 raise ValueError(msg
)
537 return self
.from_vector(coords
)
541 Return a string representation of ``self``.
545 sage: from mjo.eja.eja_algebra import JordanSpinEJA
549 Ensure that it says what we think it says::
551 sage: JordanSpinEJA(2, field=AA)
552 Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
553 sage: JordanSpinEJA(3, field=RDF)
554 Euclidean Jordan algebra of dimension 3 over Real Double Field
557 fmt
= "Euclidean Jordan algebra of dimension {} over {}"
558 return fmt
.format(self
.dimension(), self
.base_ring())
562 def characteristic_polynomial_of(self
):
564 Return the algebra's "characteristic polynomial of" function,
565 which is itself a multivariate polynomial that, when evaluated
566 at the coordinates of some algebra element, returns that
567 element's characteristic polynomial.
569 The resulting polynomial has `n+1` variables, where `n` is the
570 dimension of this algebra. The first `n` variables correspond to
571 the coordinates of an algebra element: when evaluated at the
572 coordinates of an algebra element with respect to a certain
573 basis, the result is a univariate polynomial (in the one
574 remaining variable ``t``), namely the characteristic polynomial
579 sage: from mjo.eja.eja_algebra import JordanSpinEJA, TrivialEJA
583 The characteristic polynomial in the spin algebra is given in
584 Alizadeh, Example 11.11::
586 sage: J = JordanSpinEJA(3)
587 sage: p = J.characteristic_polynomial_of(); p
588 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
589 sage: xvec = J.one().to_vector()
593 By definition, the characteristic polynomial is a monic
594 degree-zero polynomial in a rank-zero algebra. Note that
595 Cayley-Hamilton is indeed satisfied since the polynomial
596 ``1`` evaluates to the identity element of the algebra on
599 sage: J = TrivialEJA()
600 sage: J.characteristic_polynomial_of()
607 # The list of coefficient polynomials a_0, a_1, a_2, ..., a_(r-1).
608 a
= self
._charpoly
_coefficients
()
610 # We go to a bit of trouble here to reorder the
611 # indeterminates, so that it's easier to evaluate the
612 # characteristic polynomial at x's coordinates and get back
613 # something in terms of t, which is what we want.
614 S
= PolynomialRing(self
.base_ring(),'t')
618 S
= PolynomialRing(S
, R
.variable_names())
621 return (t
**r
+ sum( a
[k
]*(t
**k
) for k
in range(r
) ))
623 def coordinate_polynomial_ring(self
):
625 The multivariate polynomial ring in which this algebra's
626 :meth:`characteristic_polynomial_of` lives.
630 sage: from mjo.eja.eja_algebra import (HadamardEJA,
631 ....: RealSymmetricEJA)
635 sage: J = HadamardEJA(2)
636 sage: J.coordinate_polynomial_ring()
637 Multivariate Polynomial Ring in X1, X2...
638 sage: J = RealSymmetricEJA(3,field=QQ,orthonormalize=False)
639 sage: J.coordinate_polynomial_ring()
640 Multivariate Polynomial Ring in X1, X2, X3, X4, X5, X6...
643 var_names
= tuple( "X%d" % z
for z
in range(1, self
.dimension()+1) )
644 return PolynomialRing(self
.base_ring(), var_names
)
646 def inner_product(self
, x
, y
):
648 The inner product associated with this Euclidean Jordan algebra.
650 Defaults to the trace inner product, but can be overridden by
651 subclasses if they are sure that the necessary properties are
656 sage: from mjo.eja.eja_algebra import (random_eja,
658 ....: BilinearFormEJA)
662 Our inner product is "associative," which means the following for
663 a symmetric bilinear form::
665 sage: set_random_seed()
666 sage: J = random_eja()
667 sage: x,y,z = J.random_elements(3)
668 sage: (x*y).inner_product(z) == y.inner_product(x*z)
673 Ensure that this is the usual inner product for the algebras
676 sage: set_random_seed()
677 sage: J = HadamardEJA.random_instance()
678 sage: x,y = J.random_elements(2)
679 sage: actual = x.inner_product(y)
680 sage: expected = x.to_vector().inner_product(y.to_vector())
681 sage: actual == expected
684 Ensure that this is one-half of the trace inner-product in a
685 BilinearFormEJA that isn't just the reals (when ``n`` isn't
686 one). This is in Faraut and Koranyi, and also my "On the
689 sage: set_random_seed()
690 sage: J = BilinearFormEJA.random_instance()
691 sage: n = J.dimension()
692 sage: x = J.random_element()
693 sage: y = J.random_element()
694 sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
697 B
= self
._inner
_product
_matrix
698 return (B
*x
.to_vector()).inner_product(y
.to_vector())
701 def is_trivial(self
):
703 Return whether or not this algebra is trivial.
705 A trivial algebra contains only the zero element.
709 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
714 sage: J = ComplexHermitianEJA(3)
720 sage: J = TrivialEJA()
725 return self
.dimension() == 0
728 def multiplication_table(self
):
730 Return a visual representation of this algebra's multiplication
731 table (on basis elements).
735 sage: from mjo.eja.eja_algebra import JordanSpinEJA
739 sage: J = JordanSpinEJA(4)
740 sage: J.multiplication_table()
741 +----++----+----+----+----+
742 | * || e0 | e1 | e2 | e3 |
743 +====++====+====+====+====+
744 | e0 || e0 | e1 | e2 | e3 |
745 +----++----+----+----+----+
746 | e1 || e1 | e0 | 0 | 0 |
747 +----++----+----+----+----+
748 | e2 || e2 | 0 | e0 | 0 |
749 +----++----+----+----+----+
750 | e3 || e3 | 0 | 0 | e0 |
751 +----++----+----+----+----+
755 # Prepend the header row.
756 M
= [["*"] + list(self
.gens())]
758 # And to each subsequent row, prepend an entry that belongs to
759 # the left-side "header column."
760 M
+= [ [self
.gens()[i
]] + [ self
.product_on_basis(i
,j
)
764 return table(M
, header_row
=True, header_column
=True, frame
=True)
767 def matrix_basis(self
):
769 Return an (often more natural) representation of this algebras
770 basis as an ordered tuple of matrices.
772 Every finite-dimensional Euclidean Jordan Algebra is a, up to
773 Jordan isomorphism, a direct sum of five simple
774 algebras---four of which comprise Hermitian matrices. And the
775 last type of algebra can of course be thought of as `n`-by-`1`
776 column matrices (ambiguusly called column vectors) to avoid
777 special cases. As a result, matrices (and column vectors) are
778 a natural representation format for Euclidean Jordan algebra
781 But, when we construct an algebra from a basis of matrices,
782 those matrix representations are lost in favor of coordinate
783 vectors *with respect to* that basis. We could eventually
784 convert back if we tried hard enough, but having the original
785 representations handy is valuable enough that we simply store
786 them and return them from this method.
788 Why implement this for non-matrix algebras? Avoiding special
789 cases for the :class:`BilinearFormEJA` pays with simplicity in
790 its own right. But mainly, we would like to be able to assume
791 that elements of a :class:`CartesianProductEJA` can be displayed
792 nicely, without having to have special classes for direct sums
793 one of whose components was a matrix algebra.
797 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
798 ....: RealSymmetricEJA)
802 sage: J = RealSymmetricEJA(2)
804 Finite family {0: e0, 1: e1, 2: e2}
805 sage: J.matrix_basis()
807 [1 0] [ 0 0.7071067811865475?] [0 0]
808 [0 0], [0.7071067811865475? 0], [0 1]
813 sage: J = JordanSpinEJA(2)
815 Finite family {0: e0, 1: e1}
816 sage: J.matrix_basis()
822 return self
._matrix
_basis
825 def matrix_space(self
):
827 Return the matrix space in which this algebra's elements live, if
828 we think of them as matrices (including column vectors of the
831 "By default" this will be an `n`-by-`1` column-matrix space,
832 except when the algebra is trivial. There it's `n`-by-`n`
833 (where `n` is zero), to ensure that two elements of the matrix
834 space (empty matrices) can be multiplied. For algebras of
835 matrices, this returns the space in which their
836 real embeddings live.
840 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
842 ....: QuaternionHermitianEJA,
847 By default, the matrix representation is just a column-matrix
848 equivalent to the vector representation::
850 sage: J = JordanSpinEJA(3)
851 sage: J.matrix_space()
852 Full MatrixSpace of 3 by 1 dense matrices over Algebraic
855 The matrix representation in the trivial algebra is
856 zero-by-zero instead of the usual `n`-by-one::
858 sage: J = TrivialEJA()
859 sage: J.matrix_space()
860 Full MatrixSpace of 0 by 0 dense matrices over Algebraic
863 The matrix space for complex/quaternion Hermitian matrix EJA
864 is the space in which their real-embeddings live, not the
865 original complex/quaternion matrix space::
867 sage: J = ComplexHermitianEJA(2,field=QQ,orthonormalize=False)
868 sage: J.matrix_space()
869 Full MatrixSpace of 4 by 4 dense matrices over Rational Field
870 sage: J = QuaternionHermitianEJA(1,field=QQ,orthonormalize=False)
871 sage: J.matrix_space()
872 Full MatrixSpace of 4 by 4 dense matrices over Rational Field
875 if self
.is_trivial():
876 return MatrixSpace(self
.base_ring(), 0)
878 return self
.matrix_basis()[0].parent()
884 Return the unit element of this algebra.
888 sage: from mjo.eja.eja_algebra import (HadamardEJA,
893 We can compute unit element in the Hadamard EJA::
895 sage: J = HadamardEJA(5)
897 e0 + e1 + e2 + e3 + e4
899 The unit element in the Hadamard EJA is inherited in the
900 subalgebras generated by its elements::
902 sage: J = HadamardEJA(5)
904 e0 + e1 + e2 + e3 + e4
905 sage: x = sum(J.gens())
906 sage: A = x.subalgebra_generated_by(orthonormalize=False)
909 sage: A.one().superalgebra_element()
910 e0 + e1 + e2 + e3 + e4
914 The identity element acts like the identity, regardless of
915 whether or not we orthonormalize::
917 sage: set_random_seed()
918 sage: J = random_eja()
919 sage: x = J.random_element()
920 sage: J.one()*x == x and x*J.one() == x
922 sage: A = x.subalgebra_generated_by()
923 sage: y = A.random_element()
924 sage: A.one()*y == y and y*A.one() == y
929 sage: set_random_seed()
930 sage: J = random_eja(field=QQ, orthonormalize=False)
931 sage: x = J.random_element()
932 sage: J.one()*x == x and x*J.one() == x
934 sage: A = x.subalgebra_generated_by(orthonormalize=False)
935 sage: y = A.random_element()
936 sage: A.one()*y == y and y*A.one() == y
939 The matrix of the unit element's operator is the identity,
940 regardless of the base field and whether or not we
943 sage: set_random_seed()
944 sage: J = random_eja()
945 sage: actual = J.one().operator().matrix()
946 sage: expected = matrix.identity(J.base_ring(), J.dimension())
947 sage: actual == expected
949 sage: x = J.random_element()
950 sage: A = x.subalgebra_generated_by()
951 sage: actual = A.one().operator().matrix()
952 sage: expected = matrix.identity(A.base_ring(), A.dimension())
953 sage: actual == expected
958 sage: set_random_seed()
959 sage: J = random_eja(field=QQ, orthonormalize=False)
960 sage: actual = J.one().operator().matrix()
961 sage: expected = matrix.identity(J.base_ring(), J.dimension())
962 sage: actual == expected
964 sage: x = J.random_element()
965 sage: A = x.subalgebra_generated_by(orthonormalize=False)
966 sage: actual = A.one().operator().matrix()
967 sage: expected = matrix.identity(A.base_ring(), A.dimension())
968 sage: actual == expected
971 Ensure that the cached unit element (often precomputed by
972 hand) agrees with the computed one::
974 sage: set_random_seed()
975 sage: J = random_eja()
976 sage: cached = J.one()
977 sage: J.one.clear_cache()
978 sage: J.one() == cached
983 sage: set_random_seed()
984 sage: J = random_eja(field=QQ, orthonormalize=False)
985 sage: cached = J.one()
986 sage: J.one.clear_cache()
987 sage: J.one() == cached
991 # We can brute-force compute the matrices of the operators
992 # that correspond to the basis elements of this algebra.
993 # If some linear combination of those basis elements is the
994 # algebra identity, then the same linear combination of
995 # their matrices has to be the identity matrix.
997 # Of course, matrices aren't vectors in sage, so we have to
998 # appeal to the "long vectors" isometry.
999 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
1001 # Now we use basic linear algebra to find the coefficients,
1002 # of the matrices-as-vectors-linear-combination, which should
1003 # work for the original algebra basis too.
1004 A
= matrix(self
.base_ring(), oper_vecs
)
1006 # We used the isometry on the left-hand side already, but we
1007 # still need to do it for the right-hand side. Recall that we
1008 # wanted something that summed to the identity matrix.
1009 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
1011 # Now if there's an identity element in the algebra, this
1012 # should work. We solve on the left to avoid having to
1013 # transpose the matrix "A".
1014 return self
.from_vector(A
.solve_left(b
))
1017 def peirce_decomposition(self
, c
):
1019 The Peirce decomposition of this algebra relative to the
1022 In the future, this can be extended to a complete system of
1023 orthogonal idempotents.
1027 - ``c`` -- an idempotent of this algebra.
1031 A triple (J0, J5, J1) containing two subalgebras and one subspace
1034 - ``J0`` -- the algebra on the eigenspace of ``c.operator()``
1035 corresponding to the eigenvalue zero.
1037 - ``J5`` -- the eigenspace (NOT a subalgebra) of ``c.operator()``
1038 corresponding to the eigenvalue one-half.
1040 - ``J1`` -- the algebra on the eigenspace of ``c.operator()``
1041 corresponding to the eigenvalue one.
1043 These are the only possible eigenspaces for that operator, and this
1044 algebra is a direct sum of them. The spaces ``J0`` and ``J1`` are
1045 orthogonal, and are subalgebras of this algebra with the appropriate
1050 sage: from mjo.eja.eja_algebra import random_eja, RealSymmetricEJA
1054 The canonical example comes from the symmetric matrices, which
1055 decompose into diagonal and off-diagonal parts::
1057 sage: J = RealSymmetricEJA(3)
1058 sage: C = matrix(QQ, [ [1,0,0],
1062 sage: J0,J5,J1 = J.peirce_decomposition(c)
1064 Euclidean Jordan algebra of dimension 1...
1066 Vector space of degree 6 and dimension 2...
1068 Euclidean Jordan algebra of dimension 3...
1069 sage: J0.one().to_matrix()
1073 sage: orig_df = AA.options.display_format
1074 sage: AA.options.display_format = 'radical'
1075 sage: J.from_vector(J5.basis()[0]).to_matrix()
1079 sage: J.from_vector(J5.basis()[1]).to_matrix()
1083 sage: AA.options.display_format = orig_df
1084 sage: J1.one().to_matrix()
1091 Every algebra decomposes trivially with respect to its identity
1094 sage: set_random_seed()
1095 sage: J = random_eja()
1096 sage: J0,J5,J1 = J.peirce_decomposition(J.one())
1097 sage: J0.dimension() == 0 and J5.dimension() == 0
1099 sage: J1.superalgebra() == J and J1.dimension() == J.dimension()
1102 The decomposition is into eigenspaces, and its components are
1103 therefore necessarily orthogonal. Moreover, the identity
1104 elements in the two subalgebras are the projections onto their
1105 respective subspaces of the superalgebra's identity element::
1107 sage: set_random_seed()
1108 sage: J = random_eja()
1109 sage: x = J.random_element()
1110 sage: if not J.is_trivial():
1111 ....: while x.is_nilpotent():
1112 ....: x = J.random_element()
1113 sage: c = x.subalgebra_idempotent()
1114 sage: J0,J5,J1 = J.peirce_decomposition(c)
1116 sage: for (w,y,z) in zip(J0.basis(), J5.basis(), J1.basis()):
1117 ....: w = w.superalgebra_element()
1118 ....: y = J.from_vector(y)
1119 ....: z = z.superalgebra_element()
1120 ....: ipsum += w.inner_product(y).abs()
1121 ....: ipsum += w.inner_product(z).abs()
1122 ....: ipsum += y.inner_product(z).abs()
1125 sage: J1(c) == J1.one()
1127 sage: J0(J.one() - c) == J0.one()
1131 if not c
.is_idempotent():
1132 raise ValueError("element is not idempotent: %s" % c
)
1134 # Default these to what they should be if they turn out to be
1135 # trivial, because eigenspaces_left() won't return eigenvalues
1136 # corresponding to trivial spaces (e.g. it returns only the
1137 # eigenspace corresponding to lambda=1 if you take the
1138 # decomposition relative to the identity element).
1139 trivial
= self
.subalgebra(())
1140 J0
= trivial
# eigenvalue zero
1141 J5
= VectorSpace(self
.base_ring(), 0) # eigenvalue one-half
1142 J1
= trivial
# eigenvalue one
1144 for (eigval
, eigspace
) in c
.operator().matrix().right_eigenspaces():
1145 if eigval
== ~
(self
.base_ring()(2)):
1148 gens
= tuple( self
.from_vector(b
) for b
in eigspace
.basis() )
1149 subalg
= self
.subalgebra(gens
, check_axioms
=False)
1155 raise ValueError("unexpected eigenvalue: %s" % eigval
)
1160 def random_element(self
, thorough
=False):
1162 Return a random element of this algebra.
1164 Our algebra superclass method only returns a linear
1165 combination of at most two basis elements. We instead
1166 want the vector space "random element" method that
1167 returns a more diverse selection.
1171 - ``thorough`` -- (boolean; default False) whether or not we
1172 should generate irrational coefficients for the random
1173 element when our base ring is irrational; this slows the
1174 algebra operations to a crawl, but any truly random method
1178 # For a general base ring... maybe we can trust this to do the
1179 # right thing? Unlikely, but.
1180 V
= self
.vector_space()
1181 v
= V
.random_element()
1183 if self
.base_ring() is AA
:
1184 # The "random element" method of the algebraic reals is
1185 # stupid at the moment, and only returns integers between
1186 # -2 and 2, inclusive:
1188 # https://trac.sagemath.org/ticket/30875
1190 # Instead, we implement our own "random vector" method,
1191 # and then coerce that into the algebra. We use the vector
1192 # space degree here instead of the dimension because a
1193 # subalgebra could (for example) be spanned by only two
1194 # vectors, each with five coordinates. We need to
1195 # generate all five coordinates.
1197 v
*= QQbar
.random_element().real()
1199 v
*= QQ
.random_element()
1201 return self
.from_vector(V
.coordinate_vector(v
))
1203 def random_elements(self
, count
, thorough
=False):
1205 Return ``count`` random elements as a tuple.
1209 - ``thorough`` -- (boolean; default False) whether or not we
1210 should generate irrational coefficients for the random
1211 elements when our base ring is irrational; this slows the
1212 algebra operations to a crawl, but any truly random method
1217 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1221 sage: J = JordanSpinEJA(3)
1222 sage: x,y,z = J.random_elements(3)
1223 sage: all( [ x in J, y in J, z in J ])
1225 sage: len( J.random_elements(10) ) == 10
1229 return tuple( self
.random_element(thorough
)
1230 for idx
in range(count
) )
1234 def _charpoly_coefficients(self
):
1236 The `r` polynomial coefficients of the "characteristic polynomial
1241 sage: from mjo.eja.eja_algebra import random_eja
1245 The theory shows that these are all homogeneous polynomials of
1248 sage: set_random_seed()
1249 sage: J = random_eja()
1250 sage: all(p.is_homogeneous() for p in J._charpoly_coefficients())
1254 n
= self
.dimension()
1255 R
= self
.coordinate_polynomial_ring()
1257 F
= R
.fraction_field()
1260 # From a result in my book, these are the entries of the
1261 # basis representation of L_x.
1262 return sum( vars[k
]*self
.gens()[k
].operator().matrix()[i
,j
]
1265 L_x
= matrix(F
, n
, n
, L_x_i_j
)
1268 if self
.rank
.is_in_cache():
1270 # There's no need to pad the system with redundant
1271 # columns if we *know* they'll be redundant.
1274 # Compute an extra power in case the rank is equal to
1275 # the dimension (otherwise, we would stop at x^(r-1)).
1276 x_powers
= [ (L_x
**k
)*self
.one().to_vector()
1277 for k
in range(n
+1) ]
1278 A
= matrix
.column(F
, x_powers
[:n
])
1279 AE
= A
.extended_echelon_form()
1286 # The theory says that only the first "r" coefficients are
1287 # nonzero, and they actually live in the original polynomial
1288 # ring and not the fraction field. We negate them because in
1289 # the actual characteristic polynomial, they get moved to the
1290 # other side where x^r lives. We don't bother to trim A_rref
1291 # down to a square matrix and solve the resulting system,
1292 # because the upper-left r-by-r portion of A_rref is
1293 # guaranteed to be the identity matrix, so e.g.
1295 # A_rref.solve_right(Y)
1297 # would just be returning Y.
1298 return (-E
*b
)[:r
].change_ring(R
)
1303 Return the rank of this EJA.
1305 This is a cached method because we know the rank a priori for
1306 all of the algebras we can construct. Thus we can avoid the
1307 expensive ``_charpoly_coefficients()`` call unless we truly
1308 need to compute the whole characteristic polynomial.
1312 sage: from mjo.eja.eja_algebra import (HadamardEJA,
1313 ....: JordanSpinEJA,
1314 ....: RealSymmetricEJA,
1315 ....: ComplexHermitianEJA,
1316 ....: QuaternionHermitianEJA,
1321 The rank of the Jordan spin algebra is always two::
1323 sage: JordanSpinEJA(2).rank()
1325 sage: JordanSpinEJA(3).rank()
1327 sage: JordanSpinEJA(4).rank()
1330 The rank of the `n`-by-`n` Hermitian real, complex, or
1331 quaternion matrices is `n`::
1333 sage: RealSymmetricEJA(4).rank()
1335 sage: ComplexHermitianEJA(3).rank()
1337 sage: QuaternionHermitianEJA(2).rank()
1342 Ensure that every EJA that we know how to construct has a
1343 positive integer rank, unless the algebra is trivial in
1344 which case its rank will be zero::
1346 sage: set_random_seed()
1347 sage: J = random_eja()
1351 sage: r > 0 or (r == 0 and J.is_trivial())
1354 Ensure that computing the rank actually works, since the ranks
1355 of all simple algebras are known and will be cached by default::
1357 sage: set_random_seed() # long time
1358 sage: J = random_eja() # long time
1359 sage: cached = J.rank() # long time
1360 sage: J.rank.clear_cache() # long time
1361 sage: J.rank() == cached # long time
1365 return len(self
._charpoly
_coefficients
())
1368 def subalgebra(self
, basis
, **kwargs
):
1370 Create a subalgebra of this algebra from the given basis.
1372 from mjo
.eja
.eja_subalgebra
import FiniteDimensionalEJASubalgebra
1373 return FiniteDimensionalEJASubalgebra(self
, basis
, **kwargs
)
1376 def vector_space(self
):
1378 Return the vector space that underlies this algebra.
1382 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1386 sage: J = RealSymmetricEJA(2)
1387 sage: J.vector_space()
1388 Vector space of dimension 3 over...
1391 return self
.zero().to_vector().parent().ambient_vector_space()
1395 class RationalBasisEJA(FiniteDimensionalEJA
):
1397 New class for algebras whose supplied basis elements have all rational entries.
1401 sage: from mjo.eja.eja_algebra import BilinearFormEJA
1405 The supplied basis is orthonormalized by default::
1407 sage: B = matrix(QQ, [[1, 0, 0], [0, 25, -32], [0, -32, 41]])
1408 sage: J = BilinearFormEJA(B)
1409 sage: J.matrix_basis()
1426 # Abuse the check_field parameter to check that the entries of
1427 # out basis (in ambient coordinates) are in the field QQ.
1428 if not all( all(b_i
in QQ
for b_i
in b
.list()) for b
in basis
):
1429 raise TypeError("basis not rational")
1431 self
._rational
_algebra
= None
1433 # There's no point in constructing the extra algebra if this
1434 # one is already rational.
1436 # Note: the same Jordan and inner-products work here,
1437 # because they are necessarily defined with respect to
1438 # ambient coordinates and not any particular basis.
1439 self
._rational
_algebra
= FiniteDimensionalEJA(
1444 orthonormalize
=False,
1448 super().__init
__(basis
,
1452 check_field
=check_field
,
1456 def _charpoly_coefficients(self
):
1460 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
1461 ....: JordanSpinEJA)
1465 The base ring of the resulting polynomial coefficients is what
1466 it should be, and not the rationals (unless the algebra was
1467 already over the rationals)::
1469 sage: J = JordanSpinEJA(3)
1470 sage: J._charpoly_coefficients()
1471 (X1^2 - X2^2 - X3^2, -2*X1)
1472 sage: a0 = J._charpoly_coefficients()[0]
1474 Algebraic Real Field
1475 sage: a0.base_ring()
1476 Algebraic Real Field
1479 if self
._rational
_algebra
is None:
1480 # There's no need to construct *another* algebra over the
1481 # rationals if this one is already over the
1482 # rationals. Likewise, if we never orthonormalized our
1483 # basis, we might as well just use the given one.
1484 return super()._charpoly
_coefficients
()
1486 # Do the computation over the rationals. The answer will be
1487 # the same, because all we've done is a change of basis.
1488 # Then, change back from QQ to our real base ring
1489 a
= ( a_i
.change_ring(self
.base_ring())
1490 for a_i
in self
._rational
_algebra
._charpoly
_coefficients
() )
1492 if self
._deortho
_matrix
is None:
1493 # This can happen if our base ring was, say, AA and we
1494 # chose not to (or didn't need to) orthonormalize. It's
1495 # still faster to do the computations over QQ even if
1496 # the numbers in the boxes stay the same.
1499 # Otherwise, convert the coordinate variables back to the
1500 # deorthonormalized ones.
1501 R
= self
.coordinate_polynomial_ring()
1502 from sage
.modules
.free_module_element
import vector
1503 X
= vector(R
, R
.gens())
1504 BX
= self
._deortho
_matrix
*X
1506 subs_dict
= { X[i]: BX[i] for i in range(len(X)) }
1507 return tuple( a_i
.subs(subs_dict
) for a_i
in a
)
1509 class ConcreteEJA(RationalBasisEJA
):
1511 A class for the Euclidean Jordan algebras that we know by name.
1513 These are the Jordan algebras whose basis, multiplication table,
1514 rank, and so on are known a priori. More to the point, they are
1515 the Euclidean Jordan algebras for which we are able to conjure up
1516 a "random instance."
1520 sage: from mjo.eja.eja_algebra import ConcreteEJA
1524 Our basis is normalized with respect to the algebra's inner
1525 product, unless we specify otherwise::
1527 sage: set_random_seed()
1528 sage: J = ConcreteEJA.random_instance()
1529 sage: all( b.norm() == 1 for b in J.gens() )
1532 Since our basis is orthonormal with respect to the algebra's inner
1533 product, and since we know that this algebra is an EJA, any
1534 left-multiplication operator's matrix will be symmetric because
1535 natural->EJA basis representation is an isometry and within the
1536 EJA the operator is self-adjoint by the Jordan axiom::
1538 sage: set_random_seed()
1539 sage: J = ConcreteEJA.random_instance()
1540 sage: x = J.random_element()
1541 sage: x.operator().is_self_adjoint()
1546 def _max_random_instance_size():
1548 Return an integer "size" that is an upper bound on the size of
1549 this algebra when it is used in a random test
1550 case. Unfortunately, the term "size" is ambiguous -- when
1551 dealing with `R^n` under either the Hadamard or Jordan spin
1552 product, the "size" refers to the dimension `n`. When dealing
1553 with a matrix algebra (real symmetric or complex/quaternion
1554 Hermitian), it refers to the size of the matrix, which is far
1555 less than the dimension of the underlying vector space.
1557 This method must be implemented in each subclass.
1559 raise NotImplementedError
1562 def random_instance(cls
, *args
, **kwargs
):
1564 Return a random instance of this type of algebra.
1566 This method should be implemented in each subclass.
1568 from sage
.misc
.prandom
import choice
1569 eja_class
= choice(cls
.__subclasses
__())
1571 # These all bubble up to the RationalBasisEJA superclass
1572 # constructor, so any (kw)args valid there are also valid
1574 return eja_class
.random_instance(*args
, **kwargs
)
1579 def dimension_over_reals():
1581 The dimension of this matrix's base ring over the reals.
1583 The reals are dimension one over themselves, obviously; that's
1584 just `\mathbb{R}^{1}`. Likewise, the complex numbers `a + bi`
1585 have dimension two. Finally, the quaternions have dimension
1586 four over the reals.
1588 This is used to determine the size of the matrix returned from
1589 :meth:`real_embed`, among other things.
1591 raise NotImplementedError
1594 def real_embed(cls
,M
):
1596 Embed the matrix ``M`` into a space of real matrices.
1598 The matrix ``M`` can have entries in any field at the moment:
1599 the real numbers, complex numbers, or quaternions. And although
1600 they are not a field, we can probably support octonions at some
1601 point, too. This function returns a real matrix that "acts like"
1602 the original with respect to matrix multiplication; i.e.
1604 real_embed(M*N) = real_embed(M)*real_embed(N)
1607 if M
.ncols() != M
.nrows():
1608 raise ValueError("the matrix 'M' must be square")
1613 def real_unembed(cls
,M
):
1615 The inverse of :meth:`real_embed`.
1617 if M
.ncols() != M
.nrows():
1618 raise ValueError("the matrix 'M' must be square")
1619 if not ZZ(M
.nrows()).mod(cls
.dimension_over_reals()).is_zero():
1620 raise ValueError("the matrix 'M' must be a real embedding")
1624 def jordan_product(X
,Y
):
1625 return (X
*Y
+ Y
*X
)/2
1628 def trace_inner_product(cls
,X
,Y
):
1630 Compute the trace inner-product of two real-embeddings.
1634 sage: from mjo.eja.eja_algebra import (RealSymmetricEJA,
1635 ....: ComplexHermitianEJA,
1636 ....: QuaternionHermitianEJA)
1640 This gives the same answer as it would if we computed the trace
1641 from the unembedded (original) matrices::
1643 sage: set_random_seed()
1644 sage: J = RealSymmetricEJA.random_instance()
1645 sage: x,y = J.random_elements(2)
1646 sage: Xe = x.to_matrix()
1647 sage: Ye = y.to_matrix()
1648 sage: X = J.real_unembed(Xe)
1649 sage: Y = J.real_unembed(Ye)
1650 sage: expected = (X*Y).trace()
1651 sage: actual = J.trace_inner_product(Xe,Ye)
1652 sage: actual == expected
1657 sage: set_random_seed()
1658 sage: J = ComplexHermitianEJA.random_instance()
1659 sage: x,y = J.random_elements(2)
1660 sage: Xe = x.to_matrix()
1661 sage: Ye = y.to_matrix()
1662 sage: X = J.real_unembed(Xe)
1663 sage: Y = J.real_unembed(Ye)
1664 sage: expected = (X*Y).trace().real()
1665 sage: actual = J.trace_inner_product(Xe,Ye)
1666 sage: actual == expected
1671 sage: set_random_seed()
1672 sage: J = QuaternionHermitianEJA.random_instance()
1673 sage: x,y = J.random_elements(2)
1674 sage: Xe = x.to_matrix()
1675 sage: Ye = y.to_matrix()
1676 sage: X = J.real_unembed(Xe)
1677 sage: Y = J.real_unembed(Ye)
1678 sage: expected = (X*Y).trace().coefficient_tuple()[0]
1679 sage: actual = J.trace_inner_product(Xe,Ye)
1680 sage: actual == expected
1684 Xu
= cls
.real_unembed(X
)
1685 Yu
= cls
.real_unembed(Y
)
1686 tr
= (Xu
*Yu
).trace()
1689 # Works in QQ, AA, RDF, et cetera.
1691 except AttributeError:
1692 # A quaternion doesn't have a real() method, but does
1693 # have coefficient_tuple() method that returns the
1694 # coefficients of 1, i, j, and k -- in that order.
1695 return tr
.coefficient_tuple()[0]
1698 class RealMatrixEJA(MatrixEJA
):
1700 def dimension_over_reals():
1704 class RealSymmetricEJA(ConcreteEJA
, RealMatrixEJA
):
1706 The rank-n simple EJA consisting of real symmetric n-by-n
1707 matrices, the usual symmetric Jordan product, and the trace inner
1708 product. It has dimension `(n^2 + n)/2` over the reals.
1712 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1716 sage: J = RealSymmetricEJA(2)
1717 sage: e0, e1, e2 = J.gens()
1725 In theory, our "field" can be any subfield of the reals::
1727 sage: RealSymmetricEJA(2, field=RDF)
1728 Euclidean Jordan algebra of dimension 3 over Real Double Field
1729 sage: RealSymmetricEJA(2, field=RR)
1730 Euclidean Jordan algebra of dimension 3 over Real Field with
1731 53 bits of precision
1735 The dimension of this algebra is `(n^2 + n) / 2`::
1737 sage: set_random_seed()
1738 sage: n_max = RealSymmetricEJA._max_random_instance_size()
1739 sage: n = ZZ.random_element(1, n_max)
1740 sage: J = RealSymmetricEJA(n)
1741 sage: J.dimension() == (n^2 + n)/2
1744 The Jordan multiplication is what we think it is::
1746 sage: set_random_seed()
1747 sage: J = RealSymmetricEJA.random_instance()
1748 sage: x,y = J.random_elements(2)
1749 sage: actual = (x*y).to_matrix()
1750 sage: X = x.to_matrix()
1751 sage: Y = y.to_matrix()
1752 sage: expected = (X*Y + Y*X)/2
1753 sage: actual == expected
1755 sage: J(expected) == x*y
1758 We can change the generator prefix::
1760 sage: RealSymmetricEJA(3, prefix='q').gens()
1761 (q0, q1, q2, q3, q4, q5)
1763 We can construct the (trivial) algebra of rank zero::
1765 sage: RealSymmetricEJA(0)
1766 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1770 def _denormalized_basis(cls
, n
):
1772 Return a basis for the space of real symmetric n-by-n matrices.
1776 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1780 sage: set_random_seed()
1781 sage: n = ZZ.random_element(1,5)
1782 sage: B = RealSymmetricEJA._denormalized_basis(n)
1783 sage: all( M.is_symmetric() for M in B)
1787 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1791 for j
in range(i
+1):
1792 Eij
= matrix(ZZ
, n
, lambda k
,l
: k
==i
and l
==j
)
1796 Sij
= Eij
+ Eij
.transpose()
1802 def _max_random_instance_size():
1803 return 4 # Dimension 10
1806 def random_instance(cls
, **kwargs
):
1808 Return a random instance of this type of algebra.
1810 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1811 return cls(n
, **kwargs
)
1813 def __init__(self
, n
, **kwargs
):
1814 # We know this is a valid EJA, but will double-check
1815 # if the user passes check_axioms=True.
1816 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
1818 super(RealSymmetricEJA
, self
).__init
__(self
._denormalized
_basis
(n
),
1819 self
.jordan_product
,
1820 self
.trace_inner_product
,
1823 # TODO: this could be factored out somehow, but is left here
1824 # because the MatrixEJA is not presently a subclass of the
1825 # FDEJA class that defines rank() and one().
1826 self
.rank
.set_cache(n
)
1827 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
1828 self
.one
.set_cache(self(idV
))
1832 class ComplexMatrixEJA(MatrixEJA
):
1833 # A manual dictionary-cache for the complex_extension() method,
1834 # since apparently @classmethods can't also be @cached_methods.
1835 _complex_extension
= {}
1838 def complex_extension(cls
,field
):
1840 The complex field that we embed/unembed, as an extension
1841 of the given ``field``.
1843 if field
in cls
._complex
_extension
:
1844 return cls
._complex
_extension
[field
]
1846 # Sage doesn't know how to adjoin the complex "i" (the root of
1847 # x^2 + 1) to a field in a general way. Here, we just enumerate
1848 # all of the cases that I have cared to support so far.
1850 # Sage doesn't know how to embed AA into QQbar, i.e. how
1851 # to adjoin sqrt(-1) to AA.
1853 elif not field
.is_exact():
1855 F
= field
.complex_field()
1857 # Works for QQ and... maybe some other fields.
1858 R
= PolynomialRing(field
, 'z')
1860 F
= field
.extension(z
**2 + 1, 'I', embedding
=CLF(-1).sqrt())
1862 cls
._complex
_extension
[field
] = F
1866 def dimension_over_reals():
1870 def real_embed(cls
,M
):
1872 Embed the n-by-n complex matrix ``M`` into the space of real
1873 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1874 bi` to the block matrix ``[[a,b],[-b,a]]``.
1878 sage: from mjo.eja.eja_algebra import ComplexMatrixEJA
1882 sage: F = QuadraticField(-1, 'I')
1883 sage: x1 = F(4 - 2*i)
1884 sage: x2 = F(1 + 2*i)
1887 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1888 sage: ComplexMatrixEJA.real_embed(M)
1897 Embedding is a homomorphism (isomorphism, in fact)::
1899 sage: set_random_seed()
1900 sage: n = ZZ.random_element(3)
1901 sage: F = QuadraticField(-1, 'I')
1902 sage: X = random_matrix(F, n)
1903 sage: Y = random_matrix(F, n)
1904 sage: Xe = ComplexMatrixEJA.real_embed(X)
1905 sage: Ye = ComplexMatrixEJA.real_embed(Y)
1906 sage: XYe = ComplexMatrixEJA.real_embed(X*Y)
1911 super(ComplexMatrixEJA
,cls
).real_embed(M
)
1914 # We don't need any adjoined elements...
1915 field
= M
.base_ring().base_ring()
1921 blocks
.append(matrix(field
, 2, [ [ a
, b
],
1924 return matrix
.block(field
, n
, blocks
)
1928 def real_unembed(cls
,M
):
1930 The inverse of _embed_complex_matrix().
1934 sage: from mjo.eja.eja_algebra import ComplexMatrixEJA
1938 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1939 ....: [-2, 1, -4, 3],
1940 ....: [ 9, 10, 11, 12],
1941 ....: [-10, 9, -12, 11] ])
1942 sage: ComplexMatrixEJA.real_unembed(A)
1944 [ 10*I + 9 12*I + 11]
1948 Unembedding is the inverse of embedding::
1950 sage: set_random_seed()
1951 sage: F = QuadraticField(-1, 'I')
1952 sage: M = random_matrix(F, 3)
1953 sage: Me = ComplexMatrixEJA.real_embed(M)
1954 sage: ComplexMatrixEJA.real_unembed(Me) == M
1958 super(ComplexMatrixEJA
,cls
).real_unembed(M
)
1960 d
= cls
.dimension_over_reals()
1961 F
= cls
.complex_extension(M
.base_ring())
1964 # Go top-left to bottom-right (reading order), converting every
1965 # 2-by-2 block we see to a single complex element.
1967 for k
in range(n
/d
):
1968 for j
in range(n
/d
):
1969 submat
= M
[d
*k
:d
*k
+d
,d
*j
:d
*j
+d
]
1970 if submat
[0,0] != submat
[1,1]:
1971 raise ValueError('bad on-diagonal submatrix')
1972 if submat
[0,1] != -submat
[1,0]:
1973 raise ValueError('bad off-diagonal submatrix')
1974 z
= submat
[0,0] + submat
[0,1]*i
1977 return matrix(F
, n
/d
, elements
)
1980 class ComplexHermitianEJA(ConcreteEJA
, ComplexMatrixEJA
):
1982 The rank-n simple EJA consisting of complex Hermitian n-by-n
1983 matrices over the real numbers, the usual symmetric Jordan product,
1984 and the real-part-of-trace inner product. It has dimension `n^2` over
1989 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1993 In theory, our "field" can be any subfield of the reals::
1995 sage: ComplexHermitianEJA(2, field=RDF)
1996 Euclidean Jordan algebra of dimension 4 over Real Double Field
1997 sage: ComplexHermitianEJA(2, field=RR)
1998 Euclidean Jordan algebra of dimension 4 over Real Field with
1999 53 bits of precision
2003 The dimension of this algebra is `n^2`::
2005 sage: set_random_seed()
2006 sage: n_max = ComplexHermitianEJA._max_random_instance_size()
2007 sage: n = ZZ.random_element(1, n_max)
2008 sage: J = ComplexHermitianEJA(n)
2009 sage: J.dimension() == n^2
2012 The Jordan multiplication is what we think it is::
2014 sage: set_random_seed()
2015 sage: J = ComplexHermitianEJA.random_instance()
2016 sage: x,y = J.random_elements(2)
2017 sage: actual = (x*y).to_matrix()
2018 sage: X = x.to_matrix()
2019 sage: Y = y.to_matrix()
2020 sage: expected = (X*Y + Y*X)/2
2021 sage: actual == expected
2023 sage: J(expected) == x*y
2026 We can change the generator prefix::
2028 sage: ComplexHermitianEJA(2, prefix='z').gens()
2031 We can construct the (trivial) algebra of rank zero::
2033 sage: ComplexHermitianEJA(0)
2034 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2039 def _denormalized_basis(cls
, n
):
2041 Returns a basis for the space of complex Hermitian n-by-n matrices.
2043 Why do we embed these? Basically, because all of numerical linear
2044 algebra assumes that you're working with vectors consisting of `n`
2045 entries from a field and scalars from the same field. There's no way
2046 to tell SageMath that (for example) the vectors contain complex
2047 numbers, while the scalar field is real.
2051 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
2055 sage: set_random_seed()
2056 sage: n = ZZ.random_element(1,5)
2057 sage: B = ComplexHermitianEJA._denormalized_basis(n)
2058 sage: all( M.is_symmetric() for M in B)
2063 R
= PolynomialRing(field
, 'z')
2065 F
= field
.extension(z
**2 + 1, 'I')
2068 # This is like the symmetric case, but we need to be careful:
2070 # * We want conjugate-symmetry, not just symmetry.
2071 # * The diagonal will (as a result) be real.
2074 Eij
= matrix
.zero(F
,n
)
2076 for j
in range(i
+1):
2080 Sij
= cls
.real_embed(Eij
)
2083 # The second one has a minus because it's conjugated.
2084 Eij
[j
,i
] = 1 # Eij = Eij + Eij.transpose()
2085 Sij_real
= cls
.real_embed(Eij
)
2087 # Eij = I*Eij - I*Eij.transpose()
2090 Sij_imag
= cls
.real_embed(Eij
)
2096 # Since we embedded these, we can drop back to the "field" that we
2097 # started with instead of the complex extension "F".
2098 return tuple( s
.change_ring(field
) for s
in S
)
2101 def __init__(self
, n
, **kwargs
):
2102 # We know this is a valid EJA, but will double-check
2103 # if the user passes check_axioms=True.
2104 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2106 super(ComplexHermitianEJA
, self
).__init
__(self
._denormalized
_basis
(n
),
2107 self
.jordan_product
,
2108 self
.trace_inner_product
,
2110 # TODO: this could be factored out somehow, but is left here
2111 # because the MatrixEJA is not presently a subclass of the
2112 # FDEJA class that defines rank() and one().
2113 self
.rank
.set_cache(n
)
2114 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
2115 self
.one
.set_cache(self(idV
))
2118 def _max_random_instance_size():
2119 return 3 # Dimension 9
2122 def random_instance(cls
, **kwargs
):
2124 Return a random instance of this type of algebra.
2126 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2127 return cls(n
, **kwargs
)
2129 class QuaternionMatrixEJA(MatrixEJA
):
2131 # A manual dictionary-cache for the quaternion_extension() method,
2132 # since apparently @classmethods can't also be @cached_methods.
2133 _quaternion_extension
= {}
2136 def quaternion_extension(cls
,field
):
2138 The quaternion field that we embed/unembed, as an extension
2139 of the given ``field``.
2141 if field
in cls
._quaternion
_extension
:
2142 return cls
._quaternion
_extension
[field
]
2144 Q
= QuaternionAlgebra(field
,-1,-1)
2146 cls
._quaternion
_extension
[field
] = Q
2150 def dimension_over_reals():
2154 def real_embed(cls
,M
):
2156 Embed the n-by-n quaternion matrix ``M`` into the space of real
2157 matrices of size 4n-by-4n by first sending each quaternion entry `z
2158 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
2159 c+di],[-c + di, a-bi]]`, and then embedding those into a real
2164 sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA
2168 sage: Q = QuaternionAlgebra(QQ,-1,-1)
2169 sage: i,j,k = Q.gens()
2170 sage: x = 1 + 2*i + 3*j + 4*k
2171 sage: M = matrix(Q, 1, [[x]])
2172 sage: QuaternionMatrixEJA.real_embed(M)
2178 Embedding is a homomorphism (isomorphism, in fact)::
2180 sage: set_random_seed()
2181 sage: n = ZZ.random_element(2)
2182 sage: Q = QuaternionAlgebra(QQ,-1,-1)
2183 sage: X = random_matrix(Q, n)
2184 sage: Y = random_matrix(Q, n)
2185 sage: Xe = QuaternionMatrixEJA.real_embed(X)
2186 sage: Ye = QuaternionMatrixEJA.real_embed(Y)
2187 sage: XYe = QuaternionMatrixEJA.real_embed(X*Y)
2192 super(QuaternionMatrixEJA
,cls
).real_embed(M
)
2193 quaternions
= M
.base_ring()
2196 F
= QuadraticField(-1, 'I')
2201 t
= z
.coefficient_tuple()
2206 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
2207 [-c
+ d
*i
, a
- b
*i
]])
2208 realM
= ComplexMatrixEJA
.real_embed(cplxM
)
2209 blocks
.append(realM
)
2211 # We should have real entries by now, so use the realest field
2212 # we've got for the return value.
2213 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
2218 def real_unembed(cls
,M
):
2220 The inverse of _embed_quaternion_matrix().
2224 sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA
2228 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
2229 ....: [-2, 1, -4, 3],
2230 ....: [-3, 4, 1, -2],
2231 ....: [-4, -3, 2, 1]])
2232 sage: QuaternionMatrixEJA.real_unembed(M)
2233 [1 + 2*i + 3*j + 4*k]
2237 Unembedding is the inverse of embedding::
2239 sage: set_random_seed()
2240 sage: Q = QuaternionAlgebra(QQ, -1, -1)
2241 sage: M = random_matrix(Q, 3)
2242 sage: Me = QuaternionMatrixEJA.real_embed(M)
2243 sage: QuaternionMatrixEJA.real_unembed(Me) == M
2247 super(QuaternionMatrixEJA
,cls
).real_unembed(M
)
2249 d
= cls
.dimension_over_reals()
2251 # Use the base ring of the matrix to ensure that its entries can be
2252 # multiplied by elements of the quaternion algebra.
2253 Q
= cls
.quaternion_extension(M
.base_ring())
2256 # Go top-left to bottom-right (reading order), converting every
2257 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
2260 for l
in range(n
/d
):
2261 for m
in range(n
/d
):
2262 submat
= ComplexMatrixEJA
.real_unembed(
2263 M
[d
*l
:d
*l
+d
,d
*m
:d
*m
+d
] )
2264 if submat
[0,0] != submat
[1,1].conjugate():
2265 raise ValueError('bad on-diagonal submatrix')
2266 if submat
[0,1] != -submat
[1,0].conjugate():
2267 raise ValueError('bad off-diagonal submatrix')
2268 z
= submat
[0,0].real()
2269 z
+= submat
[0,0].imag()*i
2270 z
+= submat
[0,1].real()*j
2271 z
+= submat
[0,1].imag()*k
2274 return matrix(Q
, n
/d
, elements
)
2277 class QuaternionHermitianEJA(ConcreteEJA
, QuaternionMatrixEJA
):
2279 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
2280 matrices, the usual symmetric Jordan product, and the
2281 real-part-of-trace inner product. It has dimension `2n^2 - n` over
2286 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2290 In theory, our "field" can be any subfield of the reals::
2292 sage: QuaternionHermitianEJA(2, field=RDF)
2293 Euclidean Jordan algebra of dimension 6 over Real Double Field
2294 sage: QuaternionHermitianEJA(2, field=RR)
2295 Euclidean Jordan algebra of dimension 6 over Real Field with
2296 53 bits of precision
2300 The dimension of this algebra is `2*n^2 - n`::
2302 sage: set_random_seed()
2303 sage: n_max = QuaternionHermitianEJA._max_random_instance_size()
2304 sage: n = ZZ.random_element(1, n_max)
2305 sage: J = QuaternionHermitianEJA(n)
2306 sage: J.dimension() == 2*(n^2) - n
2309 The Jordan multiplication is what we think it is::
2311 sage: set_random_seed()
2312 sage: J = QuaternionHermitianEJA.random_instance()
2313 sage: x,y = J.random_elements(2)
2314 sage: actual = (x*y).to_matrix()
2315 sage: X = x.to_matrix()
2316 sage: Y = y.to_matrix()
2317 sage: expected = (X*Y + Y*X)/2
2318 sage: actual == expected
2320 sage: J(expected) == x*y
2323 We can change the generator prefix::
2325 sage: QuaternionHermitianEJA(2, prefix='a').gens()
2326 (a0, a1, a2, a3, a4, a5)
2328 We can construct the (trivial) algebra of rank zero::
2330 sage: QuaternionHermitianEJA(0)
2331 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2335 def _denormalized_basis(cls
, n
):
2337 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
2339 Why do we embed these? Basically, because all of numerical
2340 linear algebra assumes that you're working with vectors consisting
2341 of `n` entries from a field and scalars from the same field. There's
2342 no way to tell SageMath that (for example) the vectors contain
2343 complex numbers, while the scalar field is real.
2347 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2351 sage: set_random_seed()
2352 sage: n = ZZ.random_element(1,5)
2353 sage: B = QuaternionHermitianEJA._denormalized_basis(n)
2354 sage: all( M.is_symmetric() for M in B )
2359 Q
= QuaternionAlgebra(QQ
,-1,-1)
2362 # This is like the symmetric case, but we need to be careful:
2364 # * We want conjugate-symmetry, not just symmetry.
2365 # * The diagonal will (as a result) be real.
2368 Eij
= matrix
.zero(Q
,n
)
2370 for j
in range(i
+1):
2374 Sij
= cls
.real_embed(Eij
)
2377 # The second, third, and fourth ones have a minus
2378 # because they're conjugated.
2379 # Eij = Eij + Eij.transpose()
2381 Sij_real
= cls
.real_embed(Eij
)
2383 # Eij = I*(Eij - Eij.transpose())
2386 Sij_I
= cls
.real_embed(Eij
)
2388 # Eij = J*(Eij - Eij.transpose())
2391 Sij_J
= cls
.real_embed(Eij
)
2393 # Eij = K*(Eij - Eij.transpose())
2396 Sij_K
= cls
.real_embed(Eij
)
2402 # Since we embedded these, we can drop back to the "field" that we
2403 # started with instead of the quaternion algebra "Q".
2404 return tuple( s
.change_ring(field
) for s
in S
)
2407 def __init__(self
, n
, **kwargs
):
2408 # We know this is a valid EJA, but will double-check
2409 # if the user passes check_axioms=True.
2410 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2412 super(QuaternionHermitianEJA
, self
).__init
__(self
._denormalized
_basis
(n
),
2413 self
.jordan_product
,
2414 self
.trace_inner_product
,
2416 # TODO: this could be factored out somehow, but is left here
2417 # because the MatrixEJA is not presently a subclass of the
2418 # FDEJA class that defines rank() and one().
2419 self
.rank
.set_cache(n
)
2420 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
2421 self
.one
.set_cache(self(idV
))
2425 def _max_random_instance_size():
2427 The maximum rank of a random QuaternionHermitianEJA.
2429 return 2 # Dimension 6
2432 def random_instance(cls
, **kwargs
):
2434 Return a random instance of this type of algebra.
2436 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2437 return cls(n
, **kwargs
)
2440 class HadamardEJA(ConcreteEJA
):
2442 Return the Euclidean Jordan Algebra corresponding to the set
2443 `R^n` under the Hadamard product.
2445 Note: this is nothing more than the Cartesian product of ``n``
2446 copies of the spin algebra. Once Cartesian product algebras
2447 are implemented, this can go.
2451 sage: from mjo.eja.eja_algebra import HadamardEJA
2455 This multiplication table can be verified by hand::
2457 sage: J = HadamardEJA(3)
2458 sage: e0,e1,e2 = J.gens()
2474 We can change the generator prefix::
2476 sage: HadamardEJA(3, prefix='r').gens()
2480 def __init__(self
, n
, **kwargs
):
2482 jordan_product
= lambda x
,y
: x
2483 inner_product
= lambda x
,y
: x
2485 def jordan_product(x
,y
):
2487 return P( xi
*yi
for (xi
,yi
) in zip(x
,y
) )
2489 def inner_product(x
,y
):
2492 # New defaults for keyword arguments. Don't orthonormalize
2493 # because our basis is already orthonormal with respect to our
2494 # inner-product. Don't check the axioms, because we know this
2495 # is a valid EJA... but do double-check if the user passes
2496 # check_axioms=True. Note: we DON'T override the "check_field"
2497 # default here, because the user can pass in a field!
2498 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2499 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2501 column_basis
= tuple( b
.column() for b
in FreeModule(ZZ
, n
).basis() )
2502 super().__init
__(column_basis
,
2507 self
.rank
.set_cache(n
)
2510 self
.one
.set_cache( self
.zero() )
2512 self
.one
.set_cache( sum(self
.gens()) )
2515 def _max_random_instance_size():
2517 The maximum dimension of a random HadamardEJA.
2522 def random_instance(cls
, **kwargs
):
2524 Return a random instance of this type of algebra.
2526 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2527 return cls(n
, **kwargs
)
2530 class BilinearFormEJA(ConcreteEJA
):
2532 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2533 with the half-trace inner product and jordan product ``x*y =
2534 (<Bx,y>,y_bar>, x0*y_bar + y0*x_bar)`` where `B = 1 \times B22` is
2535 a symmetric positive-definite "bilinear form" matrix. Its
2536 dimension is the size of `B`, and it has rank two in dimensions
2537 larger than two. It reduces to the ``JordanSpinEJA`` when `B` is
2538 the identity matrix of order ``n``.
2540 We insist that the one-by-one upper-left identity block of `B` be
2541 passed in as well so that we can be passed a matrix of size zero
2542 to construct a trivial algebra.
2546 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
2547 ....: JordanSpinEJA)
2551 When no bilinear form is specified, the identity matrix is used,
2552 and the resulting algebra is the Jordan spin algebra::
2554 sage: B = matrix.identity(AA,3)
2555 sage: J0 = BilinearFormEJA(B)
2556 sage: J1 = JordanSpinEJA(3)
2557 sage: J0.multiplication_table() == J0.multiplication_table()
2560 An error is raised if the matrix `B` does not correspond to a
2561 positive-definite bilinear form::
2563 sage: B = matrix.random(QQ,2,3)
2564 sage: J = BilinearFormEJA(B)
2565 Traceback (most recent call last):
2567 ValueError: bilinear form is not positive-definite
2568 sage: B = matrix.zero(QQ,3)
2569 sage: J = BilinearFormEJA(B)
2570 Traceback (most recent call last):
2572 ValueError: bilinear form is not positive-definite
2576 We can create a zero-dimensional algebra::
2578 sage: B = matrix.identity(AA,0)
2579 sage: J = BilinearFormEJA(B)
2583 We can check the multiplication condition given in the Jordan, von
2584 Neumann, and Wigner paper (and also discussed on my "On the
2585 symmetry..." paper). Note that this relies heavily on the standard
2586 choice of basis, as does anything utilizing the bilinear form
2587 matrix. We opt not to orthonormalize the basis, because if we
2588 did, we would have to normalize the `s_{i}` in a similar manner::
2590 sage: set_random_seed()
2591 sage: n = ZZ.random_element(5)
2592 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2593 sage: B11 = matrix.identity(QQ,1)
2594 sage: B22 = M.transpose()*M
2595 sage: B = block_matrix(2,2,[ [B11,0 ],
2597 sage: J = BilinearFormEJA(B, orthonormalize=False)
2598 sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
2599 sage: V = J.vector_space()
2600 sage: sis = [ J( V([0] + (M.inverse()*ei).list()).column() )
2601 ....: for ei in eis ]
2602 sage: actual = [ sis[i]*sis[j]
2603 ....: for i in range(n-1)
2604 ....: for j in range(n-1) ]
2605 sage: expected = [ J.one() if i == j else J.zero()
2606 ....: for i in range(n-1)
2607 ....: for j in range(n-1) ]
2608 sage: actual == expected
2612 def __init__(self
, B
, **kwargs
):
2613 # The matrix "B" is supplied by the user in most cases,
2614 # so it makes sense to check whether or not its positive-
2615 # definite unless we are specifically asked not to...
2616 if ("check_axioms" not in kwargs
) or kwargs
["check_axioms"]:
2617 if not B
.is_positive_definite():
2618 raise ValueError("bilinear form is not positive-definite")
2620 # However, all of the other data for this EJA is computed
2621 # by us in manner that guarantees the axioms are
2622 # satisfied. So, again, unless we are specifically asked to
2623 # verify things, we'll skip the rest of the checks.
2624 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2626 def inner_product(x
,y
):
2627 return (y
.T
*B
*x
)[0,0]
2629 def jordan_product(x
,y
):
2635 z0
= inner_product(y
,x
)
2636 zbar
= y0
*xbar
+ x0
*ybar
2637 return P([z0
] + zbar
.list())
2640 column_basis
= tuple( b
.column() for b
in FreeModule(ZZ
, n
).basis() )
2641 super(BilinearFormEJA
, self
).__init
__(column_basis
,
2646 # The rank of this algebra is two, unless we're in a
2647 # one-dimensional ambient space (because the rank is bounded
2648 # by the ambient dimension).
2649 self
.rank
.set_cache(min(n
,2))
2652 self
.one
.set_cache( self
.zero() )
2654 self
.one
.set_cache( self
.monomial(0) )
2657 def _max_random_instance_size():
2659 The maximum dimension of a random BilinearFormEJA.
2664 def random_instance(cls
, **kwargs
):
2666 Return a random instance of this algebra.
2668 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2670 B
= matrix
.identity(ZZ
, n
)
2671 return cls(B
, **kwargs
)
2673 B11
= matrix
.identity(ZZ
, 1)
2674 M
= matrix
.random(ZZ
, n
-1)
2675 I
= matrix
.identity(ZZ
, n
-1)
2677 while alpha
.is_zero():
2678 alpha
= ZZ
.random_element().abs()
2679 B22
= M
.transpose()*M
+ alpha
*I
2681 from sage
.matrix
.special
import block_matrix
2682 B
= block_matrix(2,2, [ [B11
, ZZ(0) ],
2685 return cls(B
, **kwargs
)
2688 class JordanSpinEJA(BilinearFormEJA
):
2690 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2691 with the usual inner product and jordan product ``x*y =
2692 (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2697 sage: from mjo.eja.eja_algebra import JordanSpinEJA
2701 This multiplication table can be verified by hand::
2703 sage: J = JordanSpinEJA(4)
2704 sage: e0,e1,e2,e3 = J.gens()
2720 We can change the generator prefix::
2722 sage: JordanSpinEJA(2, prefix='B').gens()
2727 Ensure that we have the usual inner product on `R^n`::
2729 sage: set_random_seed()
2730 sage: J = JordanSpinEJA.random_instance()
2731 sage: x,y = J.random_elements(2)
2732 sage: actual = x.inner_product(y)
2733 sage: expected = x.to_vector().inner_product(y.to_vector())
2734 sage: actual == expected
2738 def __init__(self
, n
, **kwargs
):
2739 # This is a special case of the BilinearFormEJA with the
2740 # identity matrix as its bilinear form.
2741 B
= matrix
.identity(ZZ
, n
)
2743 # Don't orthonormalize because our basis is already
2744 # orthonormal with respect to our inner-product.
2745 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2747 # But also don't pass check_field=False here, because the user
2748 # can pass in a field!
2749 super(JordanSpinEJA
, self
).__init
__(B
, **kwargs
)
2752 def _max_random_instance_size():
2754 The maximum dimension of a random JordanSpinEJA.
2759 def random_instance(cls
, **kwargs
):
2761 Return a random instance of this type of algebra.
2763 Needed here to override the implementation for ``BilinearFormEJA``.
2765 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2766 return cls(n
, **kwargs
)
2769 class TrivialEJA(ConcreteEJA
):
2771 The trivial Euclidean Jordan algebra consisting of only a zero element.
2775 sage: from mjo.eja.eja_algebra import TrivialEJA
2779 sage: J = TrivialEJA()
2786 sage: 7*J.one()*12*J.one()
2788 sage: J.one().inner_product(J.one())
2790 sage: J.one().norm()
2792 sage: J.one().subalgebra_generated_by()
2793 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2798 def __init__(self
, **kwargs
):
2799 jordan_product
= lambda x
,y
: x
2800 inner_product
= lambda x
,y
: 0
2803 # New defaults for keyword arguments
2804 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2805 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2807 super(TrivialEJA
, self
).__init
__(basis
,
2811 # The rank is zero using my definition, namely the dimension of the
2812 # largest subalgebra generated by any element.
2813 self
.rank
.set_cache(0)
2814 self
.one
.set_cache( self
.zero() )
2817 def random_instance(cls
, **kwargs
):
2818 # We don't take a "size" argument so the superclass method is
2819 # inappropriate for us.
2820 return cls(**kwargs
)
2823 class CartesianProductEJA(CombinatorialFreeModule_CartesianProduct
,
2824 FiniteDimensionalEJA
):
2826 The external (orthogonal) direct sum of two or more Euclidean
2827 Jordan algebras. Every Euclidean Jordan algebra decomposes into an
2828 orthogonal direct sum of simple Euclidean Jordan algebras which is
2829 then isometric to a Cartesian product, so no generality is lost by
2830 providing only this construction.
2834 sage: from mjo.eja.eja_algebra import (random_eja,
2835 ....: CartesianProductEJA,
2837 ....: JordanSpinEJA,
2838 ....: RealSymmetricEJA)
2842 The Jordan product is inherited from our factors and implemented by
2843 our CombinatorialFreeModule Cartesian product superclass::
2845 sage: set_random_seed()
2846 sage: J1 = HadamardEJA(2)
2847 sage: J2 = RealSymmetricEJA(2)
2848 sage: J = cartesian_product([J1,J2])
2849 sage: x,y = J.random_elements(2)
2853 The ability to retrieve the original factors is implemented by our
2854 CombinatorialFreeModule Cartesian product superclass::
2856 sage: J1 = HadamardEJA(2, field=QQ)
2857 sage: J2 = JordanSpinEJA(3, field=QQ)
2858 sage: J = cartesian_product([J1,J2])
2859 sage: J.cartesian_factors()
2860 (Euclidean Jordan algebra of dimension 2 over Rational Field,
2861 Euclidean Jordan algebra of dimension 3 over Rational Field)
2863 You can provide more than two factors::
2865 sage: J1 = HadamardEJA(2)
2866 sage: J2 = JordanSpinEJA(3)
2867 sage: J3 = RealSymmetricEJA(3)
2868 sage: cartesian_product([J1,J2,J3])
2869 Euclidean Jordan algebra of dimension 2 over Algebraic Real
2870 Field (+) Euclidean Jordan algebra of dimension 3 over Algebraic
2871 Real Field (+) Euclidean Jordan algebra of dimension 6 over
2872 Algebraic Real Field
2874 Rank is additive on a Cartesian product::
2876 sage: J1 = HadamardEJA(1)
2877 sage: J2 = RealSymmetricEJA(2)
2878 sage: J = cartesian_product([J1,J2])
2879 sage: J1.rank.clear_cache()
2880 sage: J2.rank.clear_cache()
2881 sage: J.rank.clear_cache()
2884 sage: J.rank() == J1.rank() + J2.rank()
2887 The same rank computation works over the rationals, with whatever
2890 sage: J1 = HadamardEJA(1, field=QQ, orthonormalize=False)
2891 sage: J2 = RealSymmetricEJA(2, field=QQ, orthonormalize=False)
2892 sage: J = cartesian_product([J1,J2])
2893 sage: J1.rank.clear_cache()
2894 sage: J2.rank.clear_cache()
2895 sage: J.rank.clear_cache()
2898 sage: J.rank() == J1.rank() + J2.rank()
2901 The product algebra will be associative if and only if all of its
2902 components are associative::
2904 sage: J1 = HadamardEJA(2)
2905 sage: J1.is_associative()
2907 sage: J2 = HadamardEJA(3)
2908 sage: J2.is_associative()
2910 sage: J3 = RealSymmetricEJA(3)
2911 sage: J3.is_associative()
2913 sage: CP1 = cartesian_product([J1,J2])
2914 sage: CP1.is_associative()
2916 sage: CP2 = cartesian_product([J1,J3])
2917 sage: CP2.is_associative()
2922 All factors must share the same base field::
2924 sage: J1 = HadamardEJA(2, field=QQ)
2925 sage: J2 = RealSymmetricEJA(2)
2926 sage: CartesianProductEJA((J1,J2))
2927 Traceback (most recent call last):
2929 ValueError: all factors must share the same base field
2931 The cached unit element is the same one that would be computed::
2933 sage: set_random_seed() # long time
2934 sage: J1 = random_eja() # long time
2935 sage: J2 = random_eja() # long time
2936 sage: J = cartesian_product([J1,J2]) # long time
2937 sage: actual = J.one() # long time
2938 sage: J.one.clear_cache() # long time
2939 sage: expected = J.one() # long time
2940 sage: actual == expected # long time
2944 Element
= FiniteDimensionalEJAElement
2947 def __init__(self
, algebras
, **kwargs
):
2948 CombinatorialFreeModule_CartesianProduct
.__init
__(self
,
2951 field
= algebras
[0].base_ring()
2952 if not all( J
.base_ring() == field
for J
in algebras
):
2953 raise ValueError("all factors must share the same base field")
2955 associative
= all( m
.is_associative() for m
in algebras
)
2957 # The definition of matrix_space() and self.basis() relies
2958 # only on the stuff in the CFM_CartesianProduct class, which
2959 # we've already initialized.
2960 Js
= self
.cartesian_factors()
2962 MS
= self
.matrix_space()
2964 MS(tuple( self
.cartesian_projection(i
)(b
).to_matrix()
2965 for i
in range(m
) ))
2966 for b
in self
.basis()
2969 # Define jordan/inner products that operate on that matrix_basis.
2970 def jordan_product(x
,y
):
2972 (Js
[i
](x
[i
])*Js
[i
](y
[i
])).to_matrix() for i
in range(m
)
2975 def inner_product(x
, y
):
2977 Js
[i
](x
[i
]).inner_product(Js
[i
](y
[i
])) for i
in range(m
)
2980 # There's no need to check the field since it already came
2981 # from an EJA. Likewise the axioms are guaranteed to be
2982 # satisfied, unless the guy writing this class sucks.
2984 # If you want the basis to be orthonormalized, orthonormalize
2986 FiniteDimensionalEJA
.__init
__(self
,
2991 orthonormalize
=False,
2992 associative
=associative
,
2993 cartesian_product
=True,
2997 ones
= tuple(J
.one() for J
in algebras
)
2998 self
.one
.set_cache(self
._cartesian
_product
_of
_elements
(ones
))
2999 self
.rank
.set_cache(sum(J
.rank() for J
in algebras
))
3001 def matrix_space(self
):
3003 Return the space that our matrix basis lives in as a Cartesian
3008 sage: from mjo.eja.eja_algebra import (HadamardEJA,
3009 ....: RealSymmetricEJA)
3013 sage: J1 = HadamardEJA(1)
3014 sage: J2 = RealSymmetricEJA(2)
3015 sage: J = cartesian_product([J1,J2])
3016 sage: J.matrix_space()
3017 The Cartesian product of (Full MatrixSpace of 1 by 1 dense
3018 matrices over Algebraic Real Field, Full MatrixSpace of 2
3019 by 2 dense matrices over Algebraic Real Field)
3022 from sage
.categories
.cartesian_product
import cartesian_product
3023 return cartesian_product( [J
.matrix_space()
3024 for J
in self
.cartesian_factors()] )
3027 def cartesian_projection(self
, i
):
3031 sage: from mjo.eja.eja_algebra import (random_eja,
3032 ....: JordanSpinEJA,
3034 ....: RealSymmetricEJA,
3035 ....: ComplexHermitianEJA)
3039 The projection morphisms are Euclidean Jordan algebra
3042 sage: J1 = HadamardEJA(2)
3043 sage: J2 = RealSymmetricEJA(2)
3044 sage: J = cartesian_product([J1,J2])
3045 sage: J.cartesian_projection(0)
3046 Linear operator between finite-dimensional Euclidean Jordan
3047 algebras represented by the matrix:
3050 Domain: Euclidean Jordan algebra of dimension 2 over Algebraic
3051 Real Field (+) Euclidean Jordan algebra of dimension 3 over
3052 Algebraic Real Field
3053 Codomain: Euclidean Jordan algebra of dimension 2 over Algebraic
3055 sage: J.cartesian_projection(1)
3056 Linear operator between finite-dimensional Euclidean Jordan
3057 algebras represented by the matrix:
3061 Domain: Euclidean Jordan algebra of dimension 2 over Algebraic
3062 Real Field (+) Euclidean Jordan algebra of dimension 3 over
3063 Algebraic Real Field
3064 Codomain: Euclidean Jordan algebra of dimension 3 over Algebraic
3067 The projections work the way you'd expect on the vector
3068 representation of an element::
3070 sage: J1 = JordanSpinEJA(2)
3071 sage: J2 = ComplexHermitianEJA(2)
3072 sage: J = cartesian_product([J1,J2])
3073 sage: pi_left = J.cartesian_projection(0)
3074 sage: pi_right = J.cartesian_projection(1)
3075 sage: pi_left(J.one()).to_vector()
3077 sage: pi_right(J.one()).to_vector()
3079 sage: J.one().to_vector()
3084 The answer never changes::
3086 sage: set_random_seed()
3087 sage: J1 = random_eja()
3088 sage: J2 = random_eja()
3089 sage: J = cartesian_product([J1,J2])
3090 sage: P0 = J.cartesian_projection(0)
3091 sage: P1 = J.cartesian_projection(0)
3096 Ji
= self
.cartesian_factors()[i
]
3097 # Requires the fix on Trac 31421/31422 to work!
3098 Pi
= super().cartesian_projection(i
)
3099 return FiniteDimensionalEJAOperator(self
,Ji
,Pi
.matrix())
3102 def cartesian_embedding(self
, i
):
3106 sage: from mjo.eja.eja_algebra import (random_eja,
3107 ....: JordanSpinEJA,
3109 ....: RealSymmetricEJA)
3113 The embedding morphisms are Euclidean Jordan algebra
3116 sage: J1 = HadamardEJA(2)
3117 sage: J2 = RealSymmetricEJA(2)
3118 sage: J = cartesian_product([J1,J2])
3119 sage: J.cartesian_embedding(0)
3120 Linear operator between finite-dimensional Euclidean Jordan
3121 algebras represented by the matrix:
3127 Domain: Euclidean Jordan algebra of dimension 2 over
3128 Algebraic Real Field
3129 Codomain: Euclidean Jordan algebra of dimension 2 over
3130 Algebraic Real Field (+) Euclidean Jordan algebra of
3131 dimension 3 over Algebraic Real Field
3132 sage: J.cartesian_embedding(1)
3133 Linear operator between finite-dimensional Euclidean Jordan
3134 algebras represented by the matrix:
3140 Domain: Euclidean Jordan algebra of dimension 3 over
3141 Algebraic Real Field
3142 Codomain: Euclidean Jordan algebra of dimension 2 over
3143 Algebraic Real Field (+) Euclidean Jordan algebra of
3144 dimension 3 over Algebraic Real Field
3146 The embeddings work the way you'd expect on the vector
3147 representation of an element::
3149 sage: J1 = JordanSpinEJA(3)
3150 sage: J2 = RealSymmetricEJA(2)
3151 sage: J = cartesian_product([J1,J2])
3152 sage: iota_left = J.cartesian_embedding(0)
3153 sage: iota_right = J.cartesian_embedding(1)
3154 sage: iota_left(J1.zero()) == J.zero()
3156 sage: iota_right(J2.zero()) == J.zero()
3158 sage: J1.one().to_vector()
3160 sage: iota_left(J1.one()).to_vector()
3162 sage: J2.one().to_vector()
3164 sage: iota_right(J2.one()).to_vector()
3166 sage: J.one().to_vector()
3171 The answer never changes::
3173 sage: set_random_seed()
3174 sage: J1 = random_eja()
3175 sage: J2 = random_eja()
3176 sage: J = cartesian_product([J1,J2])
3177 sage: E0 = J.cartesian_embedding(0)
3178 sage: E1 = J.cartesian_embedding(0)
3182 Composing a projection with the corresponding inclusion should
3183 produce the identity map, and mismatching them should produce
3186 sage: set_random_seed()
3187 sage: J1 = random_eja()
3188 sage: J2 = random_eja()
3189 sage: J = cartesian_product([J1,J2])
3190 sage: iota_left = J.cartesian_embedding(0)
3191 sage: iota_right = J.cartesian_embedding(1)
3192 sage: pi_left = J.cartesian_projection(0)
3193 sage: pi_right = J.cartesian_projection(1)
3194 sage: pi_left*iota_left == J1.one().operator()
3196 sage: pi_right*iota_right == J2.one().operator()
3198 sage: (pi_left*iota_right).is_zero()
3200 sage: (pi_right*iota_left).is_zero()
3204 Ji
= self
.cartesian_factors()[i
]
3205 # Requires the fix on Trac 31421/31422 to work!
3206 Ei
= super().cartesian_embedding(i
)
3207 return FiniteDimensionalEJAOperator(Ji
,self
,Ei
.matrix())
3211 FiniteDimensionalEJA
.CartesianProduct
= CartesianProductEJA
3213 class RationalBasisCartesianProductEJA(CartesianProductEJA
,
3216 A separate class for products of algebras for which we know a
3221 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
3222 ....: RealSymmetricEJA)
3226 This gives us fast characteristic polynomial computations in
3227 product algebras, too::
3230 sage: J1 = JordanSpinEJA(2)
3231 sage: J2 = RealSymmetricEJA(3)
3232 sage: J = cartesian_product([J1,J2])
3233 sage: J.characteristic_polynomial_of().degree()
3239 def __init__(self
, algebras
, **kwargs
):
3240 CartesianProductEJA
.__init
__(self
, algebras
, **kwargs
)
3242 self
._rational
_algebra
= None
3243 if self
.vector_space().base_field() is not QQ
:
3244 self
._rational
_algebra
= cartesian_product([
3245 r
._rational
_algebra
for r
in algebras
3249 RationalBasisEJA
.CartesianProduct
= RationalBasisCartesianProductEJA
3251 random_eja
= ConcreteEJA
.random_instance