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 elements (in matrix form)
55 that returns their jordan product in this algebra; this will
56 be applied to ``basis`` to compute a multiplication table for
59 - inner_product -- function of two elements (in matrix form) that
60 returns their inner product. This will be applied to ``basis`` to
61 compute an inner-product table (basically a 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 Generally this will be an `n`-by-`1` column-vector 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.
836 Matrix algebras override this with something more useful.
838 if self
.is_trivial():
839 return MatrixSpace(self
.base_ring(), 0)
841 return self
.matrix_basis()[0].parent()
847 Return the unit element of this algebra.
851 sage: from mjo.eja.eja_algebra import (HadamardEJA,
856 We can compute unit element in the Hadamard EJA::
858 sage: J = HadamardEJA(5)
860 e0 + e1 + e2 + e3 + e4
862 The unit element in the Hadamard EJA is inherited in the
863 subalgebras generated by its elements::
865 sage: J = HadamardEJA(5)
867 e0 + e1 + e2 + e3 + e4
868 sage: x = sum(J.gens())
869 sage: A = x.subalgebra_generated_by(orthonormalize=False)
872 sage: A.one().superalgebra_element()
873 e0 + e1 + e2 + e3 + e4
877 The identity element acts like the identity, regardless of
878 whether or not we orthonormalize::
880 sage: set_random_seed()
881 sage: J = random_eja()
882 sage: x = J.random_element()
883 sage: J.one()*x == x and x*J.one() == x
885 sage: A = x.subalgebra_generated_by()
886 sage: y = A.random_element()
887 sage: A.one()*y == y and y*A.one() == y
892 sage: set_random_seed()
893 sage: J = random_eja(field=QQ, orthonormalize=False)
894 sage: x = J.random_element()
895 sage: J.one()*x == x and x*J.one() == x
897 sage: A = x.subalgebra_generated_by(orthonormalize=False)
898 sage: y = A.random_element()
899 sage: A.one()*y == y and y*A.one() == y
902 The matrix of the unit element's operator is the identity,
903 regardless of the base field and whether or not we
906 sage: set_random_seed()
907 sage: J = random_eja()
908 sage: actual = J.one().operator().matrix()
909 sage: expected = matrix.identity(J.base_ring(), J.dimension())
910 sage: actual == expected
912 sage: x = J.random_element()
913 sage: A = x.subalgebra_generated_by()
914 sage: actual = A.one().operator().matrix()
915 sage: expected = matrix.identity(A.base_ring(), A.dimension())
916 sage: actual == expected
921 sage: set_random_seed()
922 sage: J = random_eja(field=QQ, orthonormalize=False)
923 sage: actual = J.one().operator().matrix()
924 sage: expected = matrix.identity(J.base_ring(), J.dimension())
925 sage: actual == expected
927 sage: x = J.random_element()
928 sage: A = x.subalgebra_generated_by(orthonormalize=False)
929 sage: actual = A.one().operator().matrix()
930 sage: expected = matrix.identity(A.base_ring(), A.dimension())
931 sage: actual == expected
934 Ensure that the cached unit element (often precomputed by
935 hand) agrees with the computed one::
937 sage: set_random_seed()
938 sage: J = random_eja()
939 sage: cached = J.one()
940 sage: J.one.clear_cache()
941 sage: J.one() == cached
946 sage: set_random_seed()
947 sage: J = random_eja(field=QQ, orthonormalize=False)
948 sage: cached = J.one()
949 sage: J.one.clear_cache()
950 sage: J.one() == cached
954 # We can brute-force compute the matrices of the operators
955 # that correspond to the basis elements of this algebra.
956 # If some linear combination of those basis elements is the
957 # algebra identity, then the same linear combination of
958 # their matrices has to be the identity matrix.
960 # Of course, matrices aren't vectors in sage, so we have to
961 # appeal to the "long vectors" isometry.
962 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
964 # Now we use basic linear algebra to find the coefficients,
965 # of the matrices-as-vectors-linear-combination, which should
966 # work for the original algebra basis too.
967 A
= matrix(self
.base_ring(), oper_vecs
)
969 # We used the isometry on the left-hand side already, but we
970 # still need to do it for the right-hand side. Recall that we
971 # wanted something that summed to the identity matrix.
972 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
974 # Now if there's an identity element in the algebra, this
975 # should work. We solve on the left to avoid having to
976 # transpose the matrix "A".
977 return self
.from_vector(A
.solve_left(b
))
980 def peirce_decomposition(self
, c
):
982 The Peirce decomposition of this algebra relative to the
985 In the future, this can be extended to a complete system of
986 orthogonal idempotents.
990 - ``c`` -- an idempotent of this algebra.
994 A triple (J0, J5, J1) containing two subalgebras and one subspace
997 - ``J0`` -- the algebra on the eigenspace of ``c.operator()``
998 corresponding to the eigenvalue zero.
1000 - ``J5`` -- the eigenspace (NOT a subalgebra) of ``c.operator()``
1001 corresponding to the eigenvalue one-half.
1003 - ``J1`` -- the algebra on the eigenspace of ``c.operator()``
1004 corresponding to the eigenvalue one.
1006 These are the only possible eigenspaces for that operator, and this
1007 algebra is a direct sum of them. The spaces ``J0`` and ``J1`` are
1008 orthogonal, and are subalgebras of this algebra with the appropriate
1013 sage: from mjo.eja.eja_algebra import random_eja, RealSymmetricEJA
1017 The canonical example comes from the symmetric matrices, which
1018 decompose into diagonal and off-diagonal parts::
1020 sage: J = RealSymmetricEJA(3)
1021 sage: C = matrix(QQ, [ [1,0,0],
1025 sage: J0,J5,J1 = J.peirce_decomposition(c)
1027 Euclidean Jordan algebra of dimension 1...
1029 Vector space of degree 6 and dimension 2...
1031 Euclidean Jordan algebra of dimension 3...
1032 sage: J0.one().to_matrix()
1036 sage: orig_df = AA.options.display_format
1037 sage: AA.options.display_format = 'radical'
1038 sage: J.from_vector(J5.basis()[0]).to_matrix()
1042 sage: J.from_vector(J5.basis()[1]).to_matrix()
1046 sage: AA.options.display_format = orig_df
1047 sage: J1.one().to_matrix()
1054 Every algebra decomposes trivially with respect to its identity
1057 sage: set_random_seed()
1058 sage: J = random_eja()
1059 sage: J0,J5,J1 = J.peirce_decomposition(J.one())
1060 sage: J0.dimension() == 0 and J5.dimension() == 0
1062 sage: J1.superalgebra() == J and J1.dimension() == J.dimension()
1065 The decomposition is into eigenspaces, and its components are
1066 therefore necessarily orthogonal. Moreover, the identity
1067 elements in the two subalgebras are the projections onto their
1068 respective subspaces of the superalgebra's identity element::
1070 sage: set_random_seed()
1071 sage: J = random_eja()
1072 sage: x = J.random_element()
1073 sage: if not J.is_trivial():
1074 ....: while x.is_nilpotent():
1075 ....: x = J.random_element()
1076 sage: c = x.subalgebra_idempotent()
1077 sage: J0,J5,J1 = J.peirce_decomposition(c)
1079 sage: for (w,y,z) in zip(J0.basis(), J5.basis(), J1.basis()):
1080 ....: w = w.superalgebra_element()
1081 ....: y = J.from_vector(y)
1082 ....: z = z.superalgebra_element()
1083 ....: ipsum += w.inner_product(y).abs()
1084 ....: ipsum += w.inner_product(z).abs()
1085 ....: ipsum += y.inner_product(z).abs()
1088 sage: J1(c) == J1.one()
1090 sage: J0(J.one() - c) == J0.one()
1094 if not c
.is_idempotent():
1095 raise ValueError("element is not idempotent: %s" % c
)
1097 # Default these to what they should be if they turn out to be
1098 # trivial, because eigenspaces_left() won't return eigenvalues
1099 # corresponding to trivial spaces (e.g. it returns only the
1100 # eigenspace corresponding to lambda=1 if you take the
1101 # decomposition relative to the identity element).
1102 trivial
= self
.subalgebra(())
1103 J0
= trivial
# eigenvalue zero
1104 J5
= VectorSpace(self
.base_ring(), 0) # eigenvalue one-half
1105 J1
= trivial
# eigenvalue one
1107 for (eigval
, eigspace
) in c
.operator().matrix().right_eigenspaces():
1108 if eigval
== ~
(self
.base_ring()(2)):
1111 gens
= tuple( self
.from_vector(b
) for b
in eigspace
.basis() )
1112 subalg
= self
.subalgebra(gens
, check_axioms
=False)
1118 raise ValueError("unexpected eigenvalue: %s" % eigval
)
1123 def random_element(self
, thorough
=False):
1125 Return a random element of this algebra.
1127 Our algebra superclass method only returns a linear
1128 combination of at most two basis elements. We instead
1129 want the vector space "random element" method that
1130 returns a more diverse selection.
1134 - ``thorough`` -- (boolean; default False) whether or not we
1135 should generate irrational coefficients for the random
1136 element when our base ring is irrational; this slows the
1137 algebra operations to a crawl, but any truly random method
1141 # For a general base ring... maybe we can trust this to do the
1142 # right thing? Unlikely, but.
1143 V
= self
.vector_space()
1144 v
= V
.random_element()
1146 if self
.base_ring() is AA
:
1147 # The "random element" method of the algebraic reals is
1148 # stupid at the moment, and only returns integers between
1149 # -2 and 2, inclusive:
1151 # https://trac.sagemath.org/ticket/30875
1153 # Instead, we implement our own "random vector" method,
1154 # and then coerce that into the algebra. We use the vector
1155 # space degree here instead of the dimension because a
1156 # subalgebra could (for example) be spanned by only two
1157 # vectors, each with five coordinates. We need to
1158 # generate all five coordinates.
1160 v
*= QQbar
.random_element().real()
1162 v
*= QQ
.random_element()
1164 return self
.from_vector(V
.coordinate_vector(v
))
1166 def random_elements(self
, count
, thorough
=False):
1168 Return ``count`` random elements as a tuple.
1172 - ``thorough`` -- (boolean; default False) whether or not we
1173 should generate irrational coefficients for the random
1174 elements when our base ring is irrational; this slows the
1175 algebra operations to a crawl, but any truly random method
1180 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1184 sage: J = JordanSpinEJA(3)
1185 sage: x,y,z = J.random_elements(3)
1186 sage: all( [ x in J, y in J, z in J ])
1188 sage: len( J.random_elements(10) ) == 10
1192 return tuple( self
.random_element(thorough
)
1193 for idx
in range(count
) )
1197 def _charpoly_coefficients(self
):
1199 The `r` polynomial coefficients of the "characteristic polynomial
1204 sage: from mjo.eja.eja_algebra import random_eja
1208 The theory shows that these are all homogeneous polynomials of
1211 sage: set_random_seed()
1212 sage: J = random_eja()
1213 sage: all(p.is_homogeneous() for p in J._charpoly_coefficients())
1217 n
= self
.dimension()
1218 R
= self
.coordinate_polynomial_ring()
1220 F
= R
.fraction_field()
1223 # From a result in my book, these are the entries of the
1224 # basis representation of L_x.
1225 return sum( vars[k
]*self
.gens()[k
].operator().matrix()[i
,j
]
1228 L_x
= matrix(F
, n
, n
, L_x_i_j
)
1231 if self
.rank
.is_in_cache():
1233 # There's no need to pad the system with redundant
1234 # columns if we *know* they'll be redundant.
1237 # Compute an extra power in case the rank is equal to
1238 # the dimension (otherwise, we would stop at x^(r-1)).
1239 x_powers
= [ (L_x
**k
)*self
.one().to_vector()
1240 for k
in range(n
+1) ]
1241 A
= matrix
.column(F
, x_powers
[:n
])
1242 AE
= A
.extended_echelon_form()
1249 # The theory says that only the first "r" coefficients are
1250 # nonzero, and they actually live in the original polynomial
1251 # ring and not the fraction field. We negate them because in
1252 # the actual characteristic polynomial, they get moved to the
1253 # other side where x^r lives. We don't bother to trim A_rref
1254 # down to a square matrix and solve the resulting system,
1255 # because the upper-left r-by-r portion of A_rref is
1256 # guaranteed to be the identity matrix, so e.g.
1258 # A_rref.solve_right(Y)
1260 # would just be returning Y.
1261 return (-E
*b
)[:r
].change_ring(R
)
1266 Return the rank of this EJA.
1268 This is a cached method because we know the rank a priori for
1269 all of the algebras we can construct. Thus we can avoid the
1270 expensive ``_charpoly_coefficients()`` call unless we truly
1271 need to compute the whole characteristic polynomial.
1275 sage: from mjo.eja.eja_algebra import (HadamardEJA,
1276 ....: JordanSpinEJA,
1277 ....: RealSymmetricEJA,
1278 ....: ComplexHermitianEJA,
1279 ....: QuaternionHermitianEJA,
1284 The rank of the Jordan spin algebra is always two::
1286 sage: JordanSpinEJA(2).rank()
1288 sage: JordanSpinEJA(3).rank()
1290 sage: JordanSpinEJA(4).rank()
1293 The rank of the `n`-by-`n` Hermitian real, complex, or
1294 quaternion matrices is `n`::
1296 sage: RealSymmetricEJA(4).rank()
1298 sage: ComplexHermitianEJA(3).rank()
1300 sage: QuaternionHermitianEJA(2).rank()
1305 Ensure that every EJA that we know how to construct has a
1306 positive integer rank, unless the algebra is trivial in
1307 which case its rank will be zero::
1309 sage: set_random_seed()
1310 sage: J = random_eja()
1314 sage: r > 0 or (r == 0 and J.is_trivial())
1317 Ensure that computing the rank actually works, since the ranks
1318 of all simple algebras are known and will be cached by default::
1320 sage: set_random_seed() # long time
1321 sage: J = random_eja() # long time
1322 sage: cached = J.rank() # long time
1323 sage: J.rank.clear_cache() # long time
1324 sage: J.rank() == cached # long time
1328 return len(self
._charpoly
_coefficients
())
1331 def subalgebra(self
, basis
, **kwargs
):
1333 Create a subalgebra of this algebra from the given basis.
1335 from mjo
.eja
.eja_subalgebra
import FiniteDimensionalEJASubalgebra
1336 return FiniteDimensionalEJASubalgebra(self
, basis
, **kwargs
)
1339 def vector_space(self
):
1341 Return the vector space that underlies this algebra.
1345 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1349 sage: J = RealSymmetricEJA(2)
1350 sage: J.vector_space()
1351 Vector space of dimension 3 over...
1354 return self
.zero().to_vector().parent().ambient_vector_space()
1358 class RationalBasisEJA(FiniteDimensionalEJA
):
1360 New class for algebras whose supplied basis elements have all rational entries.
1364 sage: from mjo.eja.eja_algebra import BilinearFormEJA
1368 The supplied basis is orthonormalized by default::
1370 sage: B = matrix(QQ, [[1, 0, 0], [0, 25, -32], [0, -32, 41]])
1371 sage: J = BilinearFormEJA(B)
1372 sage: J.matrix_basis()
1389 # Abuse the check_field parameter to check that the entries of
1390 # out basis (in ambient coordinates) are in the field QQ.
1391 if not all( all(b_i
in QQ
for b_i
in b
.list()) for b
in basis
):
1392 raise TypeError("basis not rational")
1394 self
._rational
_algebra
= None
1396 # There's no point in constructing the extra algebra if this
1397 # one is already rational.
1399 # Note: the same Jordan and inner-products work here,
1400 # because they are necessarily defined with respect to
1401 # ambient coordinates and not any particular basis.
1402 self
._rational
_algebra
= FiniteDimensionalEJA(
1407 orthonormalize
=False,
1411 super().__init
__(basis
,
1415 check_field
=check_field
,
1419 def _charpoly_coefficients(self
):
1423 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
1424 ....: JordanSpinEJA)
1428 The base ring of the resulting polynomial coefficients is what
1429 it should be, and not the rationals (unless the algebra was
1430 already over the rationals)::
1432 sage: J = JordanSpinEJA(3)
1433 sage: J._charpoly_coefficients()
1434 (X1^2 - X2^2 - X3^2, -2*X1)
1435 sage: a0 = J._charpoly_coefficients()[0]
1437 Algebraic Real Field
1438 sage: a0.base_ring()
1439 Algebraic Real Field
1442 if self
._rational
_algebra
is None:
1443 # There's no need to construct *another* algebra over the
1444 # rationals if this one is already over the
1445 # rationals. Likewise, if we never orthonormalized our
1446 # basis, we might as well just use the given one.
1447 return super()._charpoly
_coefficients
()
1449 # Do the computation over the rationals. The answer will be
1450 # the same, because all we've done is a change of basis.
1451 # Then, change back from QQ to our real base ring
1452 a
= ( a_i
.change_ring(self
.base_ring())
1453 for a_i
in self
._rational
_algebra
._charpoly
_coefficients
() )
1455 if self
._deortho
_matrix
is None:
1456 # This can happen if our base ring was, say, AA and we
1457 # chose not to (or didn't need to) orthonormalize. It's
1458 # still faster to do the computations over QQ even if
1459 # the numbers in the boxes stay the same.
1462 # Otherwise, convert the coordinate variables back to the
1463 # deorthonormalized ones.
1464 R
= self
.coordinate_polynomial_ring()
1465 from sage
.modules
.free_module_element
import vector
1466 X
= vector(R
, R
.gens())
1467 BX
= self
._deortho
_matrix
*X
1469 subs_dict
= { X[i]: BX[i] for i in range(len(X)) }
1470 return tuple( a_i
.subs(subs_dict
) for a_i
in a
)
1472 class ConcreteEJA(RationalBasisEJA
):
1474 A class for the Euclidean Jordan algebras that we know by name.
1476 These are the Jordan algebras whose basis, multiplication table,
1477 rank, and so on are known a priori. More to the point, they are
1478 the Euclidean Jordan algebras for which we are able to conjure up
1479 a "random instance."
1483 sage: from mjo.eja.eja_algebra import ConcreteEJA
1487 Our basis is normalized with respect to the algebra's inner
1488 product, unless we specify otherwise::
1490 sage: set_random_seed()
1491 sage: J = ConcreteEJA.random_instance()
1492 sage: all( b.norm() == 1 for b in J.gens() )
1495 Since our basis is orthonormal with respect to the algebra's inner
1496 product, and since we know that this algebra is an EJA, any
1497 left-multiplication operator's matrix will be symmetric because
1498 natural->EJA basis representation is an isometry and within the
1499 EJA the operator is self-adjoint by the Jordan axiom::
1501 sage: set_random_seed()
1502 sage: J = ConcreteEJA.random_instance()
1503 sage: x = J.random_element()
1504 sage: x.operator().is_self_adjoint()
1509 def _max_random_instance_size():
1511 Return an integer "size" that is an upper bound on the size of
1512 this algebra when it is used in a random test
1513 case. Unfortunately, the term "size" is ambiguous -- when
1514 dealing with `R^n` under either the Hadamard or Jordan spin
1515 product, the "size" refers to the dimension `n`. When dealing
1516 with a matrix algebra (real symmetric or complex/quaternion
1517 Hermitian), it refers to the size of the matrix, which is far
1518 less than the dimension of the underlying vector space.
1520 This method must be implemented in each subclass.
1522 raise NotImplementedError
1525 def random_instance(cls
, *args
, **kwargs
):
1527 Return a random instance of this type of algebra.
1529 This method should be implemented in each subclass.
1531 from sage
.misc
.prandom
import choice
1532 eja_class
= choice(cls
.__subclasses
__())
1534 # These all bubble up to the RationalBasisEJA superclass
1535 # constructor, so any (kw)args valid there are also valid
1537 return eja_class
.random_instance(*args
, **kwargs
)
1542 def dimension_over_reals():
1544 The dimension of this matrix's base ring over the reals.
1546 The reals are dimension one over themselves, obviously; that's
1547 just `\mathbb{R}^{1}`. Likewise, the complex numbers `a + bi`
1548 have dimension two. Finally, the quaternions have dimension
1549 four over the reals.
1551 This is used to determine the size of the matrix returned from
1552 :meth:`real_embed`, among other things.
1554 raise NotImplementedError
1557 def real_embed(cls
,M
):
1559 Embed the matrix ``M`` into a space of real matrices.
1561 The matrix ``M`` can have entries in any field at the moment:
1562 the real numbers, complex numbers, or quaternions. And although
1563 they are not a field, we can probably support octonions at some
1564 point, too. This function returns a real matrix that "acts like"
1565 the original with respect to matrix multiplication; i.e.
1567 real_embed(M*N) = real_embed(M)*real_embed(N)
1570 if M
.ncols() != M
.nrows():
1571 raise ValueError("the matrix 'M' must be square")
1576 def real_unembed(cls
,M
):
1578 The inverse of :meth:`real_embed`.
1580 if M
.ncols() != M
.nrows():
1581 raise ValueError("the matrix 'M' must be square")
1582 if not ZZ(M
.nrows()).mod(cls
.dimension_over_reals()).is_zero():
1583 raise ValueError("the matrix 'M' must be a real embedding")
1587 def jordan_product(X
,Y
):
1588 return (X
*Y
+ Y
*X
)/2
1591 def trace_inner_product(cls
,X
,Y
):
1593 Compute the trace inner-product of two real-embeddings.
1597 sage: from mjo.eja.eja_algebra import (RealSymmetricEJA,
1598 ....: ComplexHermitianEJA,
1599 ....: QuaternionHermitianEJA)
1603 This gives the same answer as it would if we computed the trace
1604 from the unembedded (original) matrices::
1606 sage: set_random_seed()
1607 sage: J = RealSymmetricEJA.random_instance()
1608 sage: x,y = J.random_elements(2)
1609 sage: Xe = x.to_matrix()
1610 sage: Ye = y.to_matrix()
1611 sage: X = J.real_unembed(Xe)
1612 sage: Y = J.real_unembed(Ye)
1613 sage: expected = (X*Y).trace()
1614 sage: actual = J.trace_inner_product(Xe,Ye)
1615 sage: actual == expected
1620 sage: set_random_seed()
1621 sage: J = ComplexHermitianEJA.random_instance()
1622 sage: x,y = J.random_elements(2)
1623 sage: Xe = x.to_matrix()
1624 sage: Ye = y.to_matrix()
1625 sage: X = J.real_unembed(Xe)
1626 sage: Y = J.real_unembed(Ye)
1627 sage: expected = (X*Y).trace().real()
1628 sage: actual = J.trace_inner_product(Xe,Ye)
1629 sage: actual == expected
1634 sage: set_random_seed()
1635 sage: J = QuaternionHermitianEJA.random_instance()
1636 sage: x,y = J.random_elements(2)
1637 sage: Xe = x.to_matrix()
1638 sage: Ye = y.to_matrix()
1639 sage: X = J.real_unembed(Xe)
1640 sage: Y = J.real_unembed(Ye)
1641 sage: expected = (X*Y).trace().coefficient_tuple()[0]
1642 sage: actual = J.trace_inner_product(Xe,Ye)
1643 sage: actual == expected
1647 Xu
= cls
.real_unembed(X
)
1648 Yu
= cls
.real_unembed(Y
)
1649 tr
= (Xu
*Yu
).trace()
1652 # Works in QQ, AA, RDF, et cetera.
1654 except AttributeError:
1655 # A quaternion doesn't have a real() method, but does
1656 # have coefficient_tuple() method that returns the
1657 # coefficients of 1, i, j, and k -- in that order.
1658 return tr
.coefficient_tuple()[0]
1661 class RealMatrixEJA(MatrixEJA
):
1663 def dimension_over_reals():
1667 class RealSymmetricEJA(ConcreteEJA
, RealMatrixEJA
):
1669 The rank-n simple EJA consisting of real symmetric n-by-n
1670 matrices, the usual symmetric Jordan product, and the trace inner
1671 product. It has dimension `(n^2 + n)/2` over the reals.
1675 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1679 sage: J = RealSymmetricEJA(2)
1680 sage: e0, e1, e2 = J.gens()
1688 In theory, our "field" can be any subfield of the reals::
1690 sage: RealSymmetricEJA(2, field=RDF)
1691 Euclidean Jordan algebra of dimension 3 over Real Double Field
1692 sage: RealSymmetricEJA(2, field=RR)
1693 Euclidean Jordan algebra of dimension 3 over Real Field with
1694 53 bits of precision
1698 The dimension of this algebra is `(n^2 + n) / 2`::
1700 sage: set_random_seed()
1701 sage: n_max = RealSymmetricEJA._max_random_instance_size()
1702 sage: n = ZZ.random_element(1, n_max)
1703 sage: J = RealSymmetricEJA(n)
1704 sage: J.dimension() == (n^2 + n)/2
1707 The Jordan multiplication is what we think it is::
1709 sage: set_random_seed()
1710 sage: J = RealSymmetricEJA.random_instance()
1711 sage: x,y = J.random_elements(2)
1712 sage: actual = (x*y).to_matrix()
1713 sage: X = x.to_matrix()
1714 sage: Y = y.to_matrix()
1715 sage: expected = (X*Y + Y*X)/2
1716 sage: actual == expected
1718 sage: J(expected) == x*y
1721 We can change the generator prefix::
1723 sage: RealSymmetricEJA(3, prefix='q').gens()
1724 (q0, q1, q2, q3, q4, q5)
1726 We can construct the (trivial) algebra of rank zero::
1728 sage: RealSymmetricEJA(0)
1729 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1733 def _denormalized_basis(cls
, n
):
1735 Return a basis for the space of real symmetric n-by-n matrices.
1739 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1743 sage: set_random_seed()
1744 sage: n = ZZ.random_element(1,5)
1745 sage: B = RealSymmetricEJA._denormalized_basis(n)
1746 sage: all( M.is_symmetric() for M in B)
1750 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1754 for j
in range(i
+1):
1755 Eij
= matrix(ZZ
, n
, lambda k
,l
: k
==i
and l
==j
)
1759 Sij
= Eij
+ Eij
.transpose()
1765 def _max_random_instance_size():
1766 return 4 # Dimension 10
1769 def random_instance(cls
, **kwargs
):
1771 Return a random instance of this type of algebra.
1773 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1774 return cls(n
, **kwargs
)
1776 def __init__(self
, n
, **kwargs
):
1777 # We know this is a valid EJA, but will double-check
1778 # if the user passes check_axioms=True.
1779 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
1781 super(RealSymmetricEJA
, self
).__init
__(self
._denormalized
_basis
(n
),
1782 self
.jordan_product
,
1783 self
.trace_inner_product
,
1786 # TODO: this could be factored out somehow, but is left here
1787 # because the MatrixEJA is not presently a subclass of the
1788 # FDEJA class that defines rank() and one().
1789 self
.rank
.set_cache(n
)
1790 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
1791 self
.one
.set_cache(self(idV
))
1795 class ComplexMatrixEJA(MatrixEJA
):
1796 # A manual dictionary-cache for the complex_extension() method,
1797 # since apparently @classmethods can't also be @cached_methods.
1798 _complex_extension
= {}
1801 def complex_extension(cls
,field
):
1803 The complex field that we embed/unembed, as an extension
1804 of the given ``field``.
1806 if field
in cls
._complex
_extension
:
1807 return cls
._complex
_extension
[field
]
1809 # Sage doesn't know how to adjoin the complex "i" (the root of
1810 # x^2 + 1) to a field in a general way. Here, we just enumerate
1811 # all of the cases that I have cared to support so far.
1813 # Sage doesn't know how to embed AA into QQbar, i.e. how
1814 # to adjoin sqrt(-1) to AA.
1816 elif not field
.is_exact():
1818 F
= field
.complex_field()
1820 # Works for QQ and... maybe some other fields.
1821 R
= PolynomialRing(field
, 'z')
1823 F
= field
.extension(z
**2 + 1, 'I', embedding
=CLF(-1).sqrt())
1825 cls
._complex
_extension
[field
] = F
1829 def dimension_over_reals():
1833 def real_embed(cls
,M
):
1835 Embed the n-by-n complex matrix ``M`` into the space of real
1836 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1837 bi` to the block matrix ``[[a,b],[-b,a]]``.
1841 sage: from mjo.eja.eja_algebra import ComplexMatrixEJA
1845 sage: F = QuadraticField(-1, 'I')
1846 sage: x1 = F(4 - 2*i)
1847 sage: x2 = F(1 + 2*i)
1850 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1851 sage: ComplexMatrixEJA.real_embed(M)
1860 Embedding is a homomorphism (isomorphism, in fact)::
1862 sage: set_random_seed()
1863 sage: n = ZZ.random_element(3)
1864 sage: F = QuadraticField(-1, 'I')
1865 sage: X = random_matrix(F, n)
1866 sage: Y = random_matrix(F, n)
1867 sage: Xe = ComplexMatrixEJA.real_embed(X)
1868 sage: Ye = ComplexMatrixEJA.real_embed(Y)
1869 sage: XYe = ComplexMatrixEJA.real_embed(X*Y)
1874 super(ComplexMatrixEJA
,cls
).real_embed(M
)
1877 # We don't need any adjoined elements...
1878 field
= M
.base_ring().base_ring()
1884 blocks
.append(matrix(field
, 2, [ [ a
, b
],
1887 return matrix
.block(field
, n
, blocks
)
1891 def real_unembed(cls
,M
):
1893 The inverse of _embed_complex_matrix().
1897 sage: from mjo.eja.eja_algebra import ComplexMatrixEJA
1901 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1902 ....: [-2, 1, -4, 3],
1903 ....: [ 9, 10, 11, 12],
1904 ....: [-10, 9, -12, 11] ])
1905 sage: ComplexMatrixEJA.real_unembed(A)
1907 [ 10*I + 9 12*I + 11]
1911 Unembedding is the inverse of embedding::
1913 sage: set_random_seed()
1914 sage: F = QuadraticField(-1, 'I')
1915 sage: M = random_matrix(F, 3)
1916 sage: Me = ComplexMatrixEJA.real_embed(M)
1917 sage: ComplexMatrixEJA.real_unembed(Me) == M
1921 super(ComplexMatrixEJA
,cls
).real_unembed(M
)
1923 d
= cls
.dimension_over_reals()
1924 F
= cls
.complex_extension(M
.base_ring())
1927 # Go top-left to bottom-right (reading order), converting every
1928 # 2-by-2 block we see to a single complex element.
1930 for k
in range(n
/d
):
1931 for j
in range(n
/d
):
1932 submat
= M
[d
*k
:d
*k
+d
,d
*j
:d
*j
+d
]
1933 if submat
[0,0] != submat
[1,1]:
1934 raise ValueError('bad on-diagonal submatrix')
1935 if submat
[0,1] != -submat
[1,0]:
1936 raise ValueError('bad off-diagonal submatrix')
1937 z
= submat
[0,0] + submat
[0,1]*i
1940 return matrix(F
, n
/d
, elements
)
1943 class ComplexHermitianEJA(ConcreteEJA
, ComplexMatrixEJA
):
1945 The rank-n simple EJA consisting of complex Hermitian n-by-n
1946 matrices over the real numbers, the usual symmetric Jordan product,
1947 and the real-part-of-trace inner product. It has dimension `n^2` over
1952 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1956 In theory, our "field" can be any subfield of the reals::
1958 sage: ComplexHermitianEJA(2, field=RDF)
1959 Euclidean Jordan algebra of dimension 4 over Real Double Field
1960 sage: ComplexHermitianEJA(2, field=RR)
1961 Euclidean Jordan algebra of dimension 4 over Real Field with
1962 53 bits of precision
1966 The dimension of this algebra is `n^2`::
1968 sage: set_random_seed()
1969 sage: n_max = ComplexHermitianEJA._max_random_instance_size()
1970 sage: n = ZZ.random_element(1, n_max)
1971 sage: J = ComplexHermitianEJA(n)
1972 sage: J.dimension() == n^2
1975 The Jordan multiplication is what we think it is::
1977 sage: set_random_seed()
1978 sage: J = ComplexHermitianEJA.random_instance()
1979 sage: x,y = J.random_elements(2)
1980 sage: actual = (x*y).to_matrix()
1981 sage: X = x.to_matrix()
1982 sage: Y = y.to_matrix()
1983 sage: expected = (X*Y + Y*X)/2
1984 sage: actual == expected
1986 sage: J(expected) == x*y
1989 We can change the generator prefix::
1991 sage: ComplexHermitianEJA(2, prefix='z').gens()
1994 We can construct the (trivial) algebra of rank zero::
1996 sage: ComplexHermitianEJA(0)
1997 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2002 def _denormalized_basis(cls
, n
):
2004 Returns a basis for the space of complex Hermitian n-by-n matrices.
2006 Why do we embed these? Basically, because all of numerical linear
2007 algebra assumes that you're working with vectors consisting of `n`
2008 entries from a field and scalars from the same field. There's no way
2009 to tell SageMath that (for example) the vectors contain complex
2010 numbers, while the scalar field is real.
2014 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
2018 sage: set_random_seed()
2019 sage: n = ZZ.random_element(1,5)
2020 sage: B = ComplexHermitianEJA._denormalized_basis(n)
2021 sage: all( M.is_symmetric() for M in B)
2026 R
= PolynomialRing(field
, 'z')
2028 F
= field
.extension(z
**2 + 1, 'I')
2031 # This is like the symmetric case, but we need to be careful:
2033 # * We want conjugate-symmetry, not just symmetry.
2034 # * The diagonal will (as a result) be real.
2037 Eij
= matrix
.zero(F
,n
)
2039 for j
in range(i
+1):
2043 Sij
= cls
.real_embed(Eij
)
2046 # The second one has a minus because it's conjugated.
2047 Eij
[j
,i
] = 1 # Eij = Eij + Eij.transpose()
2048 Sij_real
= cls
.real_embed(Eij
)
2050 # Eij = I*Eij - I*Eij.transpose()
2053 Sij_imag
= cls
.real_embed(Eij
)
2059 # Since we embedded these, we can drop back to the "field" that we
2060 # started with instead of the complex extension "F".
2061 return tuple( s
.change_ring(field
) for s
in S
)
2064 def __init__(self
, n
, **kwargs
):
2065 # We know this is a valid EJA, but will double-check
2066 # if the user passes check_axioms=True.
2067 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2069 super(ComplexHermitianEJA
, self
).__init
__(self
._denormalized
_basis
(n
),
2070 self
.jordan_product
,
2071 self
.trace_inner_product
,
2073 # TODO: this could be factored out somehow, but is left here
2074 # because the MatrixEJA is not presently a subclass of the
2075 # FDEJA class that defines rank() and one().
2076 self
.rank
.set_cache(n
)
2077 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
2078 self
.one
.set_cache(self(idV
))
2081 def _max_random_instance_size():
2082 return 3 # Dimension 9
2085 def random_instance(cls
, **kwargs
):
2087 Return a random instance of this type of algebra.
2089 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2090 return cls(n
, **kwargs
)
2092 class QuaternionMatrixEJA(MatrixEJA
):
2094 # A manual dictionary-cache for the quaternion_extension() method,
2095 # since apparently @classmethods can't also be @cached_methods.
2096 _quaternion_extension
= {}
2099 def quaternion_extension(cls
,field
):
2101 The quaternion field that we embed/unembed, as an extension
2102 of the given ``field``.
2104 if field
in cls
._quaternion
_extension
:
2105 return cls
._quaternion
_extension
[field
]
2107 Q
= QuaternionAlgebra(field
,-1,-1)
2109 cls
._quaternion
_extension
[field
] = Q
2113 def dimension_over_reals():
2117 def real_embed(cls
,M
):
2119 Embed the n-by-n quaternion matrix ``M`` into the space of real
2120 matrices of size 4n-by-4n by first sending each quaternion entry `z
2121 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
2122 c+di],[-c + di, a-bi]]`, and then embedding those into a real
2127 sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA
2131 sage: Q = QuaternionAlgebra(QQ,-1,-1)
2132 sage: i,j,k = Q.gens()
2133 sage: x = 1 + 2*i + 3*j + 4*k
2134 sage: M = matrix(Q, 1, [[x]])
2135 sage: QuaternionMatrixEJA.real_embed(M)
2141 Embedding is a homomorphism (isomorphism, in fact)::
2143 sage: set_random_seed()
2144 sage: n = ZZ.random_element(2)
2145 sage: Q = QuaternionAlgebra(QQ,-1,-1)
2146 sage: X = random_matrix(Q, n)
2147 sage: Y = random_matrix(Q, n)
2148 sage: Xe = QuaternionMatrixEJA.real_embed(X)
2149 sage: Ye = QuaternionMatrixEJA.real_embed(Y)
2150 sage: XYe = QuaternionMatrixEJA.real_embed(X*Y)
2155 super(QuaternionMatrixEJA
,cls
).real_embed(M
)
2156 quaternions
= M
.base_ring()
2159 F
= QuadraticField(-1, 'I')
2164 t
= z
.coefficient_tuple()
2169 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
2170 [-c
+ d
*i
, a
- b
*i
]])
2171 realM
= ComplexMatrixEJA
.real_embed(cplxM
)
2172 blocks
.append(realM
)
2174 # We should have real entries by now, so use the realest field
2175 # we've got for the return value.
2176 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
2181 def real_unembed(cls
,M
):
2183 The inverse of _embed_quaternion_matrix().
2187 sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA
2191 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
2192 ....: [-2, 1, -4, 3],
2193 ....: [-3, 4, 1, -2],
2194 ....: [-4, -3, 2, 1]])
2195 sage: QuaternionMatrixEJA.real_unembed(M)
2196 [1 + 2*i + 3*j + 4*k]
2200 Unembedding is the inverse of embedding::
2202 sage: set_random_seed()
2203 sage: Q = QuaternionAlgebra(QQ, -1, -1)
2204 sage: M = random_matrix(Q, 3)
2205 sage: Me = QuaternionMatrixEJA.real_embed(M)
2206 sage: QuaternionMatrixEJA.real_unembed(Me) == M
2210 super(QuaternionMatrixEJA
,cls
).real_unembed(M
)
2212 d
= cls
.dimension_over_reals()
2214 # Use the base ring of the matrix to ensure that its entries can be
2215 # multiplied by elements of the quaternion algebra.
2216 Q
= cls
.quaternion_extension(M
.base_ring())
2219 # Go top-left to bottom-right (reading order), converting every
2220 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
2223 for l
in range(n
/d
):
2224 for m
in range(n
/d
):
2225 submat
= ComplexMatrixEJA
.real_unembed(
2226 M
[d
*l
:d
*l
+d
,d
*m
:d
*m
+d
] )
2227 if submat
[0,0] != submat
[1,1].conjugate():
2228 raise ValueError('bad on-diagonal submatrix')
2229 if submat
[0,1] != -submat
[1,0].conjugate():
2230 raise ValueError('bad off-diagonal submatrix')
2231 z
= submat
[0,0].real()
2232 z
+= submat
[0,0].imag()*i
2233 z
+= submat
[0,1].real()*j
2234 z
+= submat
[0,1].imag()*k
2237 return matrix(Q
, n
/d
, elements
)
2240 class QuaternionHermitianEJA(ConcreteEJA
, QuaternionMatrixEJA
):
2242 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
2243 matrices, the usual symmetric Jordan product, and the
2244 real-part-of-trace inner product. It has dimension `2n^2 - n` over
2249 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2253 In theory, our "field" can be any subfield of the reals::
2255 sage: QuaternionHermitianEJA(2, field=RDF)
2256 Euclidean Jordan algebra of dimension 6 over Real Double Field
2257 sage: QuaternionHermitianEJA(2, field=RR)
2258 Euclidean Jordan algebra of dimension 6 over Real Field with
2259 53 bits of precision
2263 The dimension of this algebra is `2*n^2 - n`::
2265 sage: set_random_seed()
2266 sage: n_max = QuaternionHermitianEJA._max_random_instance_size()
2267 sage: n = ZZ.random_element(1, n_max)
2268 sage: J = QuaternionHermitianEJA(n)
2269 sage: J.dimension() == 2*(n^2) - n
2272 The Jordan multiplication is what we think it is::
2274 sage: set_random_seed()
2275 sage: J = QuaternionHermitianEJA.random_instance()
2276 sage: x,y = J.random_elements(2)
2277 sage: actual = (x*y).to_matrix()
2278 sage: X = x.to_matrix()
2279 sage: Y = y.to_matrix()
2280 sage: expected = (X*Y + Y*X)/2
2281 sage: actual == expected
2283 sage: J(expected) == x*y
2286 We can change the generator prefix::
2288 sage: QuaternionHermitianEJA(2, prefix='a').gens()
2289 (a0, a1, a2, a3, a4, a5)
2291 We can construct the (trivial) algebra of rank zero::
2293 sage: QuaternionHermitianEJA(0)
2294 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2298 def _denormalized_basis(cls
, n
):
2300 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
2302 Why do we embed these? Basically, because all of numerical
2303 linear algebra assumes that you're working with vectors consisting
2304 of `n` entries from a field and scalars from the same field. There's
2305 no way to tell SageMath that (for example) the vectors contain
2306 complex numbers, while the scalar field is real.
2310 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2314 sage: set_random_seed()
2315 sage: n = ZZ.random_element(1,5)
2316 sage: B = QuaternionHermitianEJA._denormalized_basis(n)
2317 sage: all( M.is_symmetric() for M in B )
2322 Q
= QuaternionAlgebra(QQ
,-1,-1)
2325 # This is like the symmetric case, but we need to be careful:
2327 # * We want conjugate-symmetry, not just symmetry.
2328 # * The diagonal will (as a result) be real.
2331 Eij
= matrix
.zero(Q
,n
)
2333 for j
in range(i
+1):
2337 Sij
= cls
.real_embed(Eij
)
2340 # The second, third, and fourth ones have a minus
2341 # because they're conjugated.
2342 # Eij = Eij + Eij.transpose()
2344 Sij_real
= cls
.real_embed(Eij
)
2346 # Eij = I*(Eij - Eij.transpose())
2349 Sij_I
= cls
.real_embed(Eij
)
2351 # Eij = J*(Eij - Eij.transpose())
2354 Sij_J
= cls
.real_embed(Eij
)
2356 # Eij = K*(Eij - Eij.transpose())
2359 Sij_K
= cls
.real_embed(Eij
)
2365 # Since we embedded these, we can drop back to the "field" that we
2366 # started with instead of the quaternion algebra "Q".
2367 return tuple( s
.change_ring(field
) for s
in S
)
2370 def __init__(self
, n
, **kwargs
):
2371 # We know this is a valid EJA, but will double-check
2372 # if the user passes check_axioms=True.
2373 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2375 super(QuaternionHermitianEJA
, self
).__init
__(self
._denormalized
_basis
(n
),
2376 self
.jordan_product
,
2377 self
.trace_inner_product
,
2379 # TODO: this could be factored out somehow, but is left here
2380 # because the MatrixEJA is not presently a subclass of the
2381 # FDEJA class that defines rank() and one().
2382 self
.rank
.set_cache(n
)
2383 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
2384 self
.one
.set_cache(self(idV
))
2388 def _max_random_instance_size():
2390 The maximum rank of a random QuaternionHermitianEJA.
2392 return 2 # Dimension 6
2395 def random_instance(cls
, **kwargs
):
2397 Return a random instance of this type of algebra.
2399 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2400 return cls(n
, **kwargs
)
2403 class HadamardEJA(ConcreteEJA
):
2405 Return the Euclidean Jordan Algebra corresponding to the set
2406 `R^n` under the Hadamard product.
2408 Note: this is nothing more than the Cartesian product of ``n``
2409 copies of the spin algebra. Once Cartesian product algebras
2410 are implemented, this can go.
2414 sage: from mjo.eja.eja_algebra import HadamardEJA
2418 This multiplication table can be verified by hand::
2420 sage: J = HadamardEJA(3)
2421 sage: e0,e1,e2 = J.gens()
2437 We can change the generator prefix::
2439 sage: HadamardEJA(3, prefix='r').gens()
2443 def __init__(self
, n
, **kwargs
):
2445 jordan_product
= lambda x
,y
: x
2446 inner_product
= lambda x
,y
: x
2448 def jordan_product(x
,y
):
2450 return P( xi
*yi
for (xi
,yi
) in zip(x
,y
) )
2452 def inner_product(x
,y
):
2455 # New defaults for keyword arguments. Don't orthonormalize
2456 # because our basis is already orthonormal with respect to our
2457 # inner-product. Don't check the axioms, because we know this
2458 # is a valid EJA... but do double-check if the user passes
2459 # check_axioms=True. Note: we DON'T override the "check_field"
2460 # default here, because the user can pass in a field!
2461 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2462 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2464 column_basis
= tuple( b
.column() for b
in FreeModule(ZZ
, n
).basis() )
2465 super().__init
__(column_basis
,
2470 self
.rank
.set_cache(n
)
2473 self
.one
.set_cache( self
.zero() )
2475 self
.one
.set_cache( sum(self
.gens()) )
2478 def _max_random_instance_size():
2480 The maximum dimension of a random HadamardEJA.
2485 def random_instance(cls
, **kwargs
):
2487 Return a random instance of this type of algebra.
2489 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2490 return cls(n
, **kwargs
)
2493 class BilinearFormEJA(ConcreteEJA
):
2495 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2496 with the half-trace inner product and jordan product ``x*y =
2497 (<Bx,y>,y_bar>, x0*y_bar + y0*x_bar)`` where `B = 1 \times B22` is
2498 a symmetric positive-definite "bilinear form" matrix. Its
2499 dimension is the size of `B`, and it has rank two in dimensions
2500 larger than two. It reduces to the ``JordanSpinEJA`` when `B` is
2501 the identity matrix of order ``n``.
2503 We insist that the one-by-one upper-left identity block of `B` be
2504 passed in as well so that we can be passed a matrix of size zero
2505 to construct a trivial algebra.
2509 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
2510 ....: JordanSpinEJA)
2514 When no bilinear form is specified, the identity matrix is used,
2515 and the resulting algebra is the Jordan spin algebra::
2517 sage: B = matrix.identity(AA,3)
2518 sage: J0 = BilinearFormEJA(B)
2519 sage: J1 = JordanSpinEJA(3)
2520 sage: J0.multiplication_table() == J0.multiplication_table()
2523 An error is raised if the matrix `B` does not correspond to a
2524 positive-definite bilinear form::
2526 sage: B = matrix.random(QQ,2,3)
2527 sage: J = BilinearFormEJA(B)
2528 Traceback (most recent call last):
2530 ValueError: bilinear form is not positive-definite
2531 sage: B = matrix.zero(QQ,3)
2532 sage: J = BilinearFormEJA(B)
2533 Traceback (most recent call last):
2535 ValueError: bilinear form is not positive-definite
2539 We can create a zero-dimensional algebra::
2541 sage: B = matrix.identity(AA,0)
2542 sage: J = BilinearFormEJA(B)
2546 We can check the multiplication condition given in the Jordan, von
2547 Neumann, and Wigner paper (and also discussed on my "On the
2548 symmetry..." paper). Note that this relies heavily on the standard
2549 choice of basis, as does anything utilizing the bilinear form
2550 matrix. We opt not to orthonormalize the basis, because if we
2551 did, we would have to normalize the `s_{i}` in a similar manner::
2553 sage: set_random_seed()
2554 sage: n = ZZ.random_element(5)
2555 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2556 sage: B11 = matrix.identity(QQ,1)
2557 sage: B22 = M.transpose()*M
2558 sage: B = block_matrix(2,2,[ [B11,0 ],
2560 sage: J = BilinearFormEJA(B, orthonormalize=False)
2561 sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
2562 sage: V = J.vector_space()
2563 sage: sis = [ J( V([0] + (M.inverse()*ei).list()).column() )
2564 ....: for ei in eis ]
2565 sage: actual = [ sis[i]*sis[j]
2566 ....: for i in range(n-1)
2567 ....: for j in range(n-1) ]
2568 sage: expected = [ J.one() if i == j else J.zero()
2569 ....: for i in range(n-1)
2570 ....: for j in range(n-1) ]
2571 sage: actual == expected
2575 def __init__(self
, B
, **kwargs
):
2576 # The matrix "B" is supplied by the user in most cases,
2577 # so it makes sense to check whether or not its positive-
2578 # definite unless we are specifically asked not to...
2579 if ("check_axioms" not in kwargs
) or kwargs
["check_axioms"]:
2580 if not B
.is_positive_definite():
2581 raise ValueError("bilinear form is not positive-definite")
2583 # However, all of the other data for this EJA is computed
2584 # by us in manner that guarantees the axioms are
2585 # satisfied. So, again, unless we are specifically asked to
2586 # verify things, we'll skip the rest of the checks.
2587 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2589 def inner_product(x
,y
):
2590 return (y
.T
*B
*x
)[0,0]
2592 def jordan_product(x
,y
):
2598 z0
= inner_product(y
,x
)
2599 zbar
= y0
*xbar
+ x0
*ybar
2600 return P([z0
] + zbar
.list())
2603 column_basis
= tuple( b
.column() for b
in FreeModule(ZZ
, n
).basis() )
2604 super(BilinearFormEJA
, self
).__init
__(column_basis
,
2609 # The rank of this algebra is two, unless we're in a
2610 # one-dimensional ambient space (because the rank is bounded
2611 # by the ambient dimension).
2612 self
.rank
.set_cache(min(n
,2))
2615 self
.one
.set_cache( self
.zero() )
2617 self
.one
.set_cache( self
.monomial(0) )
2620 def _max_random_instance_size():
2622 The maximum dimension of a random BilinearFormEJA.
2627 def random_instance(cls
, **kwargs
):
2629 Return a random instance of this algebra.
2631 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2633 B
= matrix
.identity(ZZ
, n
)
2634 return cls(B
, **kwargs
)
2636 B11
= matrix
.identity(ZZ
, 1)
2637 M
= matrix
.random(ZZ
, n
-1)
2638 I
= matrix
.identity(ZZ
, n
-1)
2640 while alpha
.is_zero():
2641 alpha
= ZZ
.random_element().abs()
2642 B22
= M
.transpose()*M
+ alpha
*I
2644 from sage
.matrix
.special
import block_matrix
2645 B
= block_matrix(2,2, [ [B11
, ZZ(0) ],
2648 return cls(B
, **kwargs
)
2651 class JordanSpinEJA(BilinearFormEJA
):
2653 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2654 with the usual inner product and jordan product ``x*y =
2655 (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2660 sage: from mjo.eja.eja_algebra import JordanSpinEJA
2664 This multiplication table can be verified by hand::
2666 sage: J = JordanSpinEJA(4)
2667 sage: e0,e1,e2,e3 = J.gens()
2683 We can change the generator prefix::
2685 sage: JordanSpinEJA(2, prefix='B').gens()
2690 Ensure that we have the usual inner product on `R^n`::
2692 sage: set_random_seed()
2693 sage: J = JordanSpinEJA.random_instance()
2694 sage: x,y = J.random_elements(2)
2695 sage: actual = x.inner_product(y)
2696 sage: expected = x.to_vector().inner_product(y.to_vector())
2697 sage: actual == expected
2701 def __init__(self
, n
, **kwargs
):
2702 # This is a special case of the BilinearFormEJA with the
2703 # identity matrix as its bilinear form.
2704 B
= matrix
.identity(ZZ
, n
)
2706 # Don't orthonormalize because our basis is already
2707 # orthonormal with respect to our inner-product.
2708 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2710 # But also don't pass check_field=False here, because the user
2711 # can pass in a field!
2712 super(JordanSpinEJA
, self
).__init
__(B
, **kwargs
)
2715 def _max_random_instance_size():
2717 The maximum dimension of a random JordanSpinEJA.
2722 def random_instance(cls
, **kwargs
):
2724 Return a random instance of this type of algebra.
2726 Needed here to override the implementation for ``BilinearFormEJA``.
2728 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2729 return cls(n
, **kwargs
)
2732 class TrivialEJA(ConcreteEJA
):
2734 The trivial Euclidean Jordan algebra consisting of only a zero element.
2738 sage: from mjo.eja.eja_algebra import TrivialEJA
2742 sage: J = TrivialEJA()
2749 sage: 7*J.one()*12*J.one()
2751 sage: J.one().inner_product(J.one())
2753 sage: J.one().norm()
2755 sage: J.one().subalgebra_generated_by()
2756 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2761 def __init__(self
, **kwargs
):
2762 jordan_product
= lambda x
,y
: x
2763 inner_product
= lambda x
,y
: 0
2766 # New defaults for keyword arguments
2767 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2768 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2770 super(TrivialEJA
, self
).__init
__(basis
,
2774 # The rank is zero using my definition, namely the dimension of the
2775 # largest subalgebra generated by any element.
2776 self
.rank
.set_cache(0)
2777 self
.one
.set_cache( self
.zero() )
2780 def random_instance(cls
, **kwargs
):
2781 # We don't take a "size" argument so the superclass method is
2782 # inappropriate for us.
2783 return cls(**kwargs
)
2786 class CartesianProductEJA(CombinatorialFreeModule_CartesianProduct
,
2787 FiniteDimensionalEJA
):
2789 The external (orthogonal) direct sum of two or more Euclidean
2790 Jordan algebras. Every Euclidean Jordan algebra decomposes into an
2791 orthogonal direct sum of simple Euclidean Jordan algebras which is
2792 then isometric to a Cartesian product, so no generality is lost by
2793 providing only this construction.
2797 sage: from mjo.eja.eja_algebra import (random_eja,
2798 ....: CartesianProductEJA,
2800 ....: JordanSpinEJA,
2801 ....: RealSymmetricEJA)
2805 The Jordan product is inherited from our factors and implemented by
2806 our CombinatorialFreeModule Cartesian product superclass::
2808 sage: set_random_seed()
2809 sage: J1 = HadamardEJA(2)
2810 sage: J2 = RealSymmetricEJA(2)
2811 sage: J = cartesian_product([J1,J2])
2812 sage: x,y = J.random_elements(2)
2816 The ability to retrieve the original factors is implemented by our
2817 CombinatorialFreeModule Cartesian product superclass::
2819 sage: J1 = HadamardEJA(2, field=QQ)
2820 sage: J2 = JordanSpinEJA(3, field=QQ)
2821 sage: J = cartesian_product([J1,J2])
2822 sage: J.cartesian_factors()
2823 (Euclidean Jordan algebra of dimension 2 over Rational Field,
2824 Euclidean Jordan algebra of dimension 3 over Rational Field)
2826 You can provide more than two factors::
2828 sage: J1 = HadamardEJA(2)
2829 sage: J2 = JordanSpinEJA(3)
2830 sage: J3 = RealSymmetricEJA(3)
2831 sage: cartesian_product([J1,J2,J3])
2832 Euclidean Jordan algebra of dimension 2 over Algebraic Real
2833 Field (+) Euclidean Jordan algebra of dimension 3 over Algebraic
2834 Real Field (+) Euclidean Jordan algebra of dimension 6 over
2835 Algebraic Real Field
2837 Rank is additive on a Cartesian product::
2839 sage: J1 = HadamardEJA(1)
2840 sage: J2 = RealSymmetricEJA(2)
2841 sage: J = cartesian_product([J1,J2])
2842 sage: J1.rank.clear_cache()
2843 sage: J2.rank.clear_cache()
2844 sage: J.rank.clear_cache()
2847 sage: J.rank() == J1.rank() + J2.rank()
2850 The same rank computation works over the rationals, with whatever
2853 sage: J1 = HadamardEJA(1, field=QQ, orthonormalize=False)
2854 sage: J2 = RealSymmetricEJA(2, field=QQ, orthonormalize=False)
2855 sage: J = cartesian_product([J1,J2])
2856 sage: J1.rank.clear_cache()
2857 sage: J2.rank.clear_cache()
2858 sage: J.rank.clear_cache()
2861 sage: J.rank() == J1.rank() + J2.rank()
2864 The product algebra will be associative if and only if all of its
2865 components are associative::
2867 sage: J1 = HadamardEJA(2)
2868 sage: J1.is_associative()
2870 sage: J2 = HadamardEJA(3)
2871 sage: J2.is_associative()
2873 sage: J3 = RealSymmetricEJA(3)
2874 sage: J3.is_associative()
2876 sage: CP1 = cartesian_product([J1,J2])
2877 sage: CP1.is_associative()
2879 sage: CP2 = cartesian_product([J1,J3])
2880 sage: CP2.is_associative()
2885 All factors must share the same base field::
2887 sage: J1 = HadamardEJA(2, field=QQ)
2888 sage: J2 = RealSymmetricEJA(2)
2889 sage: CartesianProductEJA((J1,J2))
2890 Traceback (most recent call last):
2892 ValueError: all factors must share the same base field
2894 The cached unit element is the same one that would be computed::
2896 sage: set_random_seed() # long time
2897 sage: J1 = random_eja() # long time
2898 sage: J2 = random_eja() # long time
2899 sage: J = cartesian_product([J1,J2]) # long time
2900 sage: actual = J.one() # long time
2901 sage: J.one.clear_cache() # long time
2902 sage: expected = J.one() # long time
2903 sage: actual == expected # long time
2907 Element
= FiniteDimensionalEJAElement
2910 def __init__(self
, algebras
, **kwargs
):
2911 CombinatorialFreeModule_CartesianProduct
.__init
__(self
,
2914 field
= algebras
[0].base_ring()
2915 if not all( J
.base_ring() == field
for J
in algebras
):
2916 raise ValueError("all factors must share the same base field")
2918 associative
= all( m
.is_associative() for m
in algebras
)
2920 # The definition of matrix_space() and self.basis() relies
2921 # only on the stuff in the CFM_CartesianProduct class, which
2922 # we've already initialized.
2923 Js
= self
.cartesian_factors()
2925 MS
= self
.matrix_space()
2927 MS(tuple( self
.cartesian_projection(i
)(b
).to_matrix()
2928 for i
in range(m
) ))
2929 for b
in self
.basis()
2932 # Define jordan/inner products that operate on that matrix_basis.
2933 def jordan_product(x
,y
):
2935 (Js
[i
](x
[i
])*Js
[i
](y
[i
])).to_matrix() for i
in range(m
)
2938 def inner_product(x
, y
):
2940 Js
[i
](x
[i
]).inner_product(Js
[i
](y
[i
])) for i
in range(m
)
2943 # There's no need to check the field since it already came
2944 # from an EJA. Likewise the axioms are guaranteed to be
2945 # satisfied, unless the guy writing this class sucks.
2947 # If you want the basis to be orthonormalized, orthonormalize
2949 FiniteDimensionalEJA
.__init
__(self
,
2954 orthonormalize
=False,
2955 associative
=associative
,
2956 cartesian_product
=True,
2960 ones
= tuple(J
.one() for J
in algebras
)
2961 self
.one
.set_cache(self
._cartesian
_product
_of
_elements
(ones
))
2962 self
.rank
.set_cache(sum(J
.rank() for J
in algebras
))
2964 def matrix_space(self
):
2966 Return the space that our matrix basis lives in as a Cartesian
2971 sage: from mjo.eja.eja_algebra import (HadamardEJA,
2972 ....: RealSymmetricEJA)
2976 sage: J1 = HadamardEJA(1)
2977 sage: J2 = RealSymmetricEJA(2)
2978 sage: J = cartesian_product([J1,J2])
2979 sage: J.matrix_space()
2980 The Cartesian product of (Full MatrixSpace of 1 by 1 dense
2981 matrices over Algebraic Real Field, Full MatrixSpace of 2
2982 by 2 dense matrices over Algebraic Real Field)
2985 from sage
.categories
.cartesian_product
import cartesian_product
2986 return cartesian_product( [J
.matrix_space()
2987 for J
in self
.cartesian_factors()] )
2990 def cartesian_projection(self
, i
):
2994 sage: from mjo.eja.eja_algebra import (random_eja,
2995 ....: JordanSpinEJA,
2997 ....: RealSymmetricEJA,
2998 ....: ComplexHermitianEJA)
3002 The projection morphisms are Euclidean Jordan algebra
3005 sage: J1 = HadamardEJA(2)
3006 sage: J2 = RealSymmetricEJA(2)
3007 sage: J = cartesian_product([J1,J2])
3008 sage: J.cartesian_projection(0)
3009 Linear operator between finite-dimensional Euclidean Jordan
3010 algebras represented by the matrix:
3013 Domain: Euclidean Jordan algebra of dimension 2 over Algebraic
3014 Real Field (+) Euclidean Jordan algebra of dimension 3 over
3015 Algebraic Real Field
3016 Codomain: Euclidean Jordan algebra of dimension 2 over Algebraic
3018 sage: J.cartesian_projection(1)
3019 Linear operator between finite-dimensional Euclidean Jordan
3020 algebras represented by the matrix:
3024 Domain: Euclidean Jordan algebra of dimension 2 over Algebraic
3025 Real Field (+) Euclidean Jordan algebra of dimension 3 over
3026 Algebraic Real Field
3027 Codomain: Euclidean Jordan algebra of dimension 3 over Algebraic
3030 The projections work the way you'd expect on the vector
3031 representation of an element::
3033 sage: J1 = JordanSpinEJA(2)
3034 sage: J2 = ComplexHermitianEJA(2)
3035 sage: J = cartesian_product([J1,J2])
3036 sage: pi_left = J.cartesian_projection(0)
3037 sage: pi_right = J.cartesian_projection(1)
3038 sage: pi_left(J.one()).to_vector()
3040 sage: pi_right(J.one()).to_vector()
3042 sage: J.one().to_vector()
3047 The answer never changes::
3049 sage: set_random_seed()
3050 sage: J1 = random_eja()
3051 sage: J2 = random_eja()
3052 sage: J = cartesian_product([J1,J2])
3053 sage: P0 = J.cartesian_projection(0)
3054 sage: P1 = J.cartesian_projection(0)
3059 Ji
= self
.cartesian_factors()[i
]
3060 # Requires the fix on Trac 31421/31422 to work!
3061 Pi
= super().cartesian_projection(i
)
3062 return FiniteDimensionalEJAOperator(self
,Ji
,Pi
.matrix())
3065 def cartesian_embedding(self
, i
):
3069 sage: from mjo.eja.eja_algebra import (random_eja,
3070 ....: JordanSpinEJA,
3072 ....: RealSymmetricEJA)
3076 The embedding morphisms are Euclidean Jordan algebra
3079 sage: J1 = HadamardEJA(2)
3080 sage: J2 = RealSymmetricEJA(2)
3081 sage: J = cartesian_product([J1,J2])
3082 sage: J.cartesian_embedding(0)
3083 Linear operator between finite-dimensional Euclidean Jordan
3084 algebras represented by the matrix:
3090 Domain: Euclidean Jordan algebra of dimension 2 over
3091 Algebraic Real Field
3092 Codomain: Euclidean Jordan algebra of dimension 2 over
3093 Algebraic Real Field (+) Euclidean Jordan algebra of
3094 dimension 3 over Algebraic Real Field
3095 sage: J.cartesian_embedding(1)
3096 Linear operator between finite-dimensional Euclidean Jordan
3097 algebras represented by the matrix:
3103 Domain: Euclidean Jordan algebra of dimension 3 over
3104 Algebraic Real Field
3105 Codomain: Euclidean Jordan algebra of dimension 2 over
3106 Algebraic Real Field (+) Euclidean Jordan algebra of
3107 dimension 3 over Algebraic Real Field
3109 The embeddings work the way you'd expect on the vector
3110 representation of an element::
3112 sage: J1 = JordanSpinEJA(3)
3113 sage: J2 = RealSymmetricEJA(2)
3114 sage: J = cartesian_product([J1,J2])
3115 sage: iota_left = J.cartesian_embedding(0)
3116 sage: iota_right = J.cartesian_embedding(1)
3117 sage: iota_left(J1.zero()) == J.zero()
3119 sage: iota_right(J2.zero()) == J.zero()
3121 sage: J1.one().to_vector()
3123 sage: iota_left(J1.one()).to_vector()
3125 sage: J2.one().to_vector()
3127 sage: iota_right(J2.one()).to_vector()
3129 sage: J.one().to_vector()
3134 The answer never changes::
3136 sage: set_random_seed()
3137 sage: J1 = random_eja()
3138 sage: J2 = random_eja()
3139 sage: J = cartesian_product([J1,J2])
3140 sage: E0 = J.cartesian_embedding(0)
3141 sage: E1 = J.cartesian_embedding(0)
3145 Composing a projection with the corresponding inclusion should
3146 produce the identity map, and mismatching them should produce
3149 sage: set_random_seed()
3150 sage: J1 = random_eja()
3151 sage: J2 = random_eja()
3152 sage: J = cartesian_product([J1,J2])
3153 sage: iota_left = J.cartesian_embedding(0)
3154 sage: iota_right = J.cartesian_embedding(1)
3155 sage: pi_left = J.cartesian_projection(0)
3156 sage: pi_right = J.cartesian_projection(1)
3157 sage: pi_left*iota_left == J1.one().operator()
3159 sage: pi_right*iota_right == J2.one().operator()
3161 sage: (pi_left*iota_right).is_zero()
3163 sage: (pi_right*iota_left).is_zero()
3167 Ji
= self
.cartesian_factors()[i
]
3168 # Requires the fix on Trac 31421/31422 to work!
3169 Ei
= super().cartesian_embedding(i
)
3170 return FiniteDimensionalEJAOperator(Ji
,self
,Ei
.matrix())
3174 FiniteDimensionalEJA
.CartesianProduct
= CartesianProductEJA
3176 random_eja
= ConcreteEJA
.random_instance