2 Euclidean Jordan Algebras. These are formally-real Jordan Algebras;
3 specifically those where u^2 + v^2 = 0 implies that u = v = 0. They
4 are used in optimization, and have some additional nice methods beyond
5 what can be supported in a general Jordan Algebra.
10 sage: from mjo.eja.eja_algebra import random_eja
15 Euclidean Jordan algebra of dimension...
19 from itertools
import repeat
21 from sage
.algebras
.quatalg
.quaternion_algebra
import QuaternionAlgebra
22 from sage
.categories
.magmatic_algebras
import MagmaticAlgebras
23 from sage
.categories
.sets_cat
import cartesian_product
24 from sage
.combinat
.free_module
import (CombinatorialFreeModule
,
25 CombinatorialFreeModule_CartesianProduct
)
26 from sage
.matrix
.constructor
import matrix
27 from sage
.matrix
.matrix_space
import MatrixSpace
28 from sage
.misc
.cachefunc
import cached_method
29 from sage
.misc
.table
import table
30 from sage
.modules
.free_module
import FreeModule
, VectorSpace
31 from sage
.rings
.all
import (ZZ
, QQ
, AA
, QQbar
, RR
, RLF
, CLF
,
34 from mjo
.eja
.eja_element
import (CartesianProductEJAElement
,
35 FiniteDimensionalEJAElement
)
36 from mjo
.eja
.eja_operator
import FiniteDimensionalEJAOperator
37 from mjo
.eja
.eja_utils
import _mat2vec
39 class FiniteDimensionalEJA(CombinatorialFreeModule
):
41 A finite-dimensional Euclidean Jordan algebra.
45 - basis -- a tuple of basis elements in "matrix form," which
46 must be the same form as the arguments to ``jordan_product``
47 and ``inner_product``. In reality, "matrix form" can be either
48 vectors, matrices, or a Cartesian product (ordered tuple)
49 of vectors or matrices. All of these would ideally be vector
50 spaces in sage with no special-casing needed; but in reality
51 we turn vectors into column-matrices and Cartesian products
52 `(a,b)` into column matrices `(a,b)^{T}` after converting
53 `a` and `b` themselves.
55 - jordan_product -- function of two elements (in matrix form)
56 that returns their jordan product in this algebra; this will
57 be applied to ``basis`` to compute a multiplication table for
60 - inner_product -- function of two elements (in matrix form) that
61 returns their inner product. This will be applied to ``basis`` to
62 compute an inner-product table (basically a matrix) for this algebra.
65 Element
= FiniteDimensionalEJAElement
74 cartesian_product
=False,
80 if not field
.is_subring(RR
):
81 # Note: this does return true for the real algebraic
82 # field, the rationals, and any quadratic field where
83 # we've specified a real embedding.
84 raise ValueError("scalar field is not real")
86 # If the basis given to us wasn't over the field that it's
87 # supposed to be over, fix that. Or, you know, crash.
88 if not cartesian_product
:
89 # The field for a cartesian product algebra comes from one
90 # of its factors and is the same for all factors, so
91 # there's no need to "reapply" it on product algebras.
92 basis
= tuple( b
.change_ring(field
) for b
in basis
)
96 # Check commutativity of the Jordan and inner-products.
97 # This has to be done before we build the multiplication
98 # and inner-product tables/matrices, because we take
99 # advantage of symmetry in the process.
100 if not all( jordan_product(bi
,bj
) == jordan_product(bj
,bi
)
103 raise ValueError("Jordan product is not commutative")
105 if not all( inner_product(bi
,bj
) == inner_product(bj
,bi
)
108 raise ValueError("inner-product is not commutative")
111 category
= MagmaticAlgebras(field
).FiniteDimensional()
112 category
= category
.WithBasis().Unital()
114 # Element subalgebras can take advantage of this.
115 category
= category
.Associative()
116 if cartesian_product
:
117 category
= category
.CartesianProducts()
119 # Call the superclass constructor so that we can use its from_vector()
120 # method to build our multiplication table.
122 super().__init
__(field
,
128 # Now comes all of the hard work. We'll be constructing an
129 # ambient vector space V that our (vectorized) basis lives in,
130 # as well as a subspace W of V spanned by those (vectorized)
131 # basis elements. The W-coordinates are the coefficients that
132 # we see in things like x = 1*e1 + 2*e2.
136 # flatten a vector, matrix, or cartesian product of those
137 # things into a long list.
138 if cartesian_product
:
139 return sum(( b_i
.list() for b_i
in b
), [])
145 degree
= len(flatten(basis
[0]))
147 # Build an ambient space that fits our matrix basis when
148 # written out as "long vectors."
149 V
= VectorSpace(field
, degree
)
151 # The matrix that will hole the orthonormal -> unorthonormal
152 # coordinate transformation.
153 self
._deortho
_matrix
= None
156 # Save a copy of the un-orthonormalized basis for later.
157 # Convert it to ambient V (vector) coordinates while we're
158 # at it, because we'd have to do it later anyway.
159 deortho_vector_basis
= tuple( V(flatten(b
)) for b
in basis
)
161 from mjo
.eja
.eja_utils
import gram_schmidt
162 basis
= tuple(gram_schmidt(basis
, inner_product
))
164 # Save the (possibly orthonormalized) matrix basis for
166 self
._matrix
_basis
= basis
168 # Now create the vector space for the algebra, which will have
169 # its own set of non-ambient coordinates (in terms of the
171 vector_basis
= tuple( V(flatten(b
)) for b
in basis
)
172 W
= V
.span_of_basis( vector_basis
, check
=check_axioms
)
175 # Now "W" is the vector space of our algebra coordinates. The
176 # variables "X1", "X2",... refer to the entries of vectors in
177 # W. Thus to convert back and forth between the orthonormal
178 # coordinates and the given ones, we need to stick the original
180 U
= V
.span_of_basis( deortho_vector_basis
, check
=check_axioms
)
181 self
._deortho
_matrix
= matrix( U
.coordinate_vector(q
)
182 for q
in vector_basis
)
185 # Now we actually compute the multiplication and inner-product
186 # tables/matrices using the possibly-orthonormalized basis.
187 self
._inner
_product
_matrix
= matrix
.identity(field
, n
)
188 self
._multiplication
_table
= [ [0 for j
in range(i
+1)]
191 # Note: the Jordan and inner-products are defined in terms
192 # of the ambient basis. It's important that their arguments
193 # are in ambient coordinates as well.
196 # ortho basis w.r.t. ambient coords
200 # The jordan product returns a matrixy answer, so we
201 # have to convert it to the algebra coordinates.
202 elt
= jordan_product(q_i
, q_j
)
203 elt
= W
.coordinate_vector(V(flatten(elt
)))
204 self
._multiplication
_table
[i
][j
] = self
.from_vector(elt
)
206 if not orthonormalize
:
207 # If we're orthonormalizing the basis with respect
208 # to an inner-product, then the inner-product
209 # matrix with respect to the resulting basis is
210 # just going to be the identity.
211 ip
= inner_product(q_i
, q_j
)
212 self
._inner
_product
_matrix
[i
,j
] = ip
213 self
._inner
_product
_matrix
[j
,i
] = ip
215 self
._inner
_product
_matrix
._cache
= {'hermitian': True}
216 self
._inner
_product
_matrix
.set_immutable()
219 if not self
._is
_jordanian
():
220 raise ValueError("Jordan identity does not hold")
221 if not self
._inner
_product
_is
_associative
():
222 raise ValueError("inner product is not associative")
225 def _coerce_map_from_base_ring(self
):
227 Disable the map from the base ring into the algebra.
229 Performing a nonsense conversion like this automatically
230 is counterpedagogical. The fallback is to try the usual
231 element constructor, which should also fail.
235 sage: from mjo.eja.eja_algebra import random_eja
239 sage: set_random_seed()
240 sage: J = random_eja()
242 Traceback (most recent call last):
244 ValueError: not an element of this algebra
250 def product_on_basis(self
, i
, j
):
251 # We only stored the lower-triangular portion of the
252 # multiplication table.
254 return self
._multiplication
_table
[i
][j
]
256 return self
._multiplication
_table
[j
][i
]
258 def inner_product(self
, x
, y
):
260 The inner product associated with this Euclidean Jordan algebra.
262 Defaults to the trace inner product, but can be overridden by
263 subclasses if they are sure that the necessary properties are
268 sage: from mjo.eja.eja_algebra import (random_eja,
270 ....: BilinearFormEJA)
274 Our inner product is "associative," which means the following for
275 a symmetric bilinear form::
277 sage: set_random_seed()
278 sage: J = random_eja()
279 sage: x,y,z = J.random_elements(3)
280 sage: (x*y).inner_product(z) == y.inner_product(x*z)
285 Ensure that this is the usual inner product for the algebras
288 sage: set_random_seed()
289 sage: J = HadamardEJA.random_instance()
290 sage: x,y = J.random_elements(2)
291 sage: actual = x.inner_product(y)
292 sage: expected = x.to_vector().inner_product(y.to_vector())
293 sage: actual == expected
296 Ensure that this is one-half of the trace inner-product in a
297 BilinearFormEJA that isn't just the reals (when ``n`` isn't
298 one). This is in Faraut and Koranyi, and also my "On the
301 sage: set_random_seed()
302 sage: J = BilinearFormEJA.random_instance()
303 sage: n = J.dimension()
304 sage: x = J.random_element()
305 sage: y = J.random_element()
306 sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
309 B
= self
._inner
_product
_matrix
310 return (B
*x
.to_vector()).inner_product(y
.to_vector())
313 def _is_commutative(self
):
315 Whether or not this algebra's multiplication table is commutative.
317 This method should of course always return ``True``, unless
318 this algebra was constructed with ``check_axioms=False`` and
319 passed an invalid multiplication table.
321 return all( self
.product_on_basis(i
,j
) == self
.product_on_basis(i
,j
)
322 for i
in range(self
.dimension())
323 for j
in range(self
.dimension()) )
325 def _is_jordanian(self
):
327 Whether or not this algebra's multiplication table respects the
328 Jordan identity `(x^{2})(xy) = x(x^{2}y)`.
330 We only check one arrangement of `x` and `y`, so for a
331 ``True`` result to be truly true, you should also check
332 :meth:`_is_commutative`. This method should of course always
333 return ``True``, unless this algebra was constructed with
334 ``check_axioms=False`` and passed an invalid multiplication table.
336 return all( (self
.gens()[i
]**2)*(self
.gens()[i
]*self
.gens()[j
])
338 (self
.gens()[i
])*((self
.gens()[i
]**2)*self
.gens()[j
])
339 for i
in range(self
.dimension())
340 for j
in range(self
.dimension()) )
342 def _inner_product_is_associative(self
):
344 Return whether or not this algebra's inner product `B` is
345 associative; that is, whether or not `B(xy,z) = B(x,yz)`.
347 This method should of course always return ``True``, unless
348 this algebra was constructed with ``check_axioms=False`` and
349 passed an invalid Jordan or inner-product.
352 # Used to check whether or not something is zero in an inexact
353 # ring. This number is sufficient to allow the construction of
354 # QuaternionHermitianEJA(2, field=RDF) with check_axioms=True.
357 for i
in range(self
.dimension()):
358 for j
in range(self
.dimension()):
359 for k
in range(self
.dimension()):
363 diff
= (x
*y
).inner_product(z
) - x
.inner_product(y
*z
)
365 if self
.base_ring().is_exact():
369 if diff
.abs() > epsilon
:
374 def _element_constructor_(self
, elt
):
376 Construct an element of this algebra from its vector or matrix
379 This gets called only after the parent element _call_ method
380 fails to find a coercion for the argument.
384 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
386 ....: RealSymmetricEJA)
390 The identity in `S^n` is converted to the identity in the EJA::
392 sage: J = RealSymmetricEJA(3)
393 sage: I = matrix.identity(QQ,3)
394 sage: J(I) == J.one()
397 This skew-symmetric matrix can't be represented in the EJA::
399 sage: J = RealSymmetricEJA(3)
400 sage: A = matrix(QQ,3, lambda i,j: i-j)
402 Traceback (most recent call last):
404 ValueError: not an element of this algebra
408 Ensure that we can convert any element of the two non-matrix
409 simple algebras (whose matrix representations are columns)
410 back and forth faithfully::
412 sage: set_random_seed()
413 sage: J = HadamardEJA.random_instance()
414 sage: x = J.random_element()
415 sage: J(x.to_vector().column()) == x
417 sage: J = JordanSpinEJA.random_instance()
418 sage: x = J.random_element()
419 sage: J(x.to_vector().column()) == x
423 msg
= "not an element of this algebra"
425 # The superclass implementation of random_element()
426 # needs to be able to coerce "0" into the algebra.
428 elif elt
in self
.base_ring():
429 # Ensure that no base ring -> algebra coercion is performed
430 # by this method. There's some stupidity in sage that would
431 # otherwise propagate to this method; for example, sage thinks
432 # that the integer 3 belongs to the space of 2-by-2 matrices.
433 raise ValueError(msg
)
437 except (AttributeError, TypeError):
438 # Try to convert a vector into a column-matrix
441 if elt
not in self
.matrix_space():
442 raise ValueError(msg
)
444 # Thanks for nothing! Matrix spaces aren't vector spaces in
445 # Sage, so we have to figure out its matrix-basis coordinates
446 # ourselves. We use the basis space's ring instead of the
447 # element's ring because the basis space might be an algebraic
448 # closure whereas the base ring of the 3-by-3 identity matrix
449 # could be QQ instead of QQbar.
451 # We pass check=False because the matrix basis is "guaranteed"
452 # to be linearly independent... right? Ha ha.
453 V
= VectorSpace(self
.base_ring(), elt
.nrows()*elt
.ncols())
454 W
= V
.span_of_basis( (_mat2vec(s
) for s
in self
.matrix_basis()),
458 coords
= W
.coordinate_vector(_mat2vec(elt
))
459 except ArithmeticError: # vector is not in free module
460 raise ValueError(msg
)
462 return self
.from_vector(coords
)
466 Return a string representation of ``self``.
470 sage: from mjo.eja.eja_algebra import JordanSpinEJA
474 Ensure that it says what we think it says::
476 sage: JordanSpinEJA(2, field=AA)
477 Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
478 sage: JordanSpinEJA(3, field=RDF)
479 Euclidean Jordan algebra of dimension 3 over Real Double Field
482 fmt
= "Euclidean Jordan algebra of dimension {} over {}"
483 return fmt
.format(self
.dimension(), self
.base_ring())
487 def characteristic_polynomial_of(self
):
489 Return the algebra's "characteristic polynomial of" function,
490 which is itself a multivariate polynomial that, when evaluated
491 at the coordinates of some algebra element, returns that
492 element's characteristic polynomial.
494 The resulting polynomial has `n+1` variables, where `n` is the
495 dimension of this algebra. The first `n` variables correspond to
496 the coordinates of an algebra element: when evaluated at the
497 coordinates of an algebra element with respect to a certain
498 basis, the result is a univariate polynomial (in the one
499 remaining variable ``t``), namely the characteristic polynomial
504 sage: from mjo.eja.eja_algebra import JordanSpinEJA, TrivialEJA
508 The characteristic polynomial in the spin algebra is given in
509 Alizadeh, Example 11.11::
511 sage: J = JordanSpinEJA(3)
512 sage: p = J.characteristic_polynomial_of(); p
513 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
514 sage: xvec = J.one().to_vector()
518 By definition, the characteristic polynomial is a monic
519 degree-zero polynomial in a rank-zero algebra. Note that
520 Cayley-Hamilton is indeed satisfied since the polynomial
521 ``1`` evaluates to the identity element of the algebra on
524 sage: J = TrivialEJA()
525 sage: J.characteristic_polynomial_of()
532 # The list of coefficient polynomials a_0, a_1, a_2, ..., a_(r-1).
533 a
= self
._charpoly
_coefficients
()
535 # We go to a bit of trouble here to reorder the
536 # indeterminates, so that it's easier to evaluate the
537 # characteristic polynomial at x's coordinates and get back
538 # something in terms of t, which is what we want.
539 S
= PolynomialRing(self
.base_ring(),'t')
543 S
= PolynomialRing(S
, R
.variable_names())
546 return (t
**r
+ sum( a
[k
]*(t
**k
) for k
in range(r
) ))
548 def coordinate_polynomial_ring(self
):
550 The multivariate polynomial ring in which this algebra's
551 :meth:`characteristic_polynomial_of` lives.
555 sage: from mjo.eja.eja_algebra import (HadamardEJA,
556 ....: RealSymmetricEJA)
560 sage: J = HadamardEJA(2)
561 sage: J.coordinate_polynomial_ring()
562 Multivariate Polynomial Ring in X1, X2...
563 sage: J = RealSymmetricEJA(3,field=QQ,orthonormalize=False)
564 sage: J.coordinate_polynomial_ring()
565 Multivariate Polynomial Ring in X1, X2, X3, X4, X5, X6...
568 var_names
= tuple( "X%d" % z
for z
in range(1, self
.dimension()+1) )
569 return PolynomialRing(self
.base_ring(), var_names
)
571 def inner_product(self
, x
, y
):
573 The inner product associated with this Euclidean Jordan algebra.
575 Defaults to the trace inner product, but can be overridden by
576 subclasses if they are sure that the necessary properties are
581 sage: from mjo.eja.eja_algebra import (random_eja,
583 ....: BilinearFormEJA)
587 Our inner product is "associative," which means the following for
588 a symmetric bilinear form::
590 sage: set_random_seed()
591 sage: J = random_eja()
592 sage: x,y,z = J.random_elements(3)
593 sage: (x*y).inner_product(z) == y.inner_product(x*z)
598 Ensure that this is the usual inner product for the algebras
601 sage: set_random_seed()
602 sage: J = HadamardEJA.random_instance()
603 sage: x,y = J.random_elements(2)
604 sage: actual = x.inner_product(y)
605 sage: expected = x.to_vector().inner_product(y.to_vector())
606 sage: actual == expected
609 Ensure that this is one-half of the trace inner-product in a
610 BilinearFormEJA that isn't just the reals (when ``n`` isn't
611 one). This is in Faraut and Koranyi, and also my "On the
614 sage: set_random_seed()
615 sage: J = BilinearFormEJA.random_instance()
616 sage: n = J.dimension()
617 sage: x = J.random_element()
618 sage: y = J.random_element()
619 sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
622 B
= self
._inner
_product
_matrix
623 return (B
*x
.to_vector()).inner_product(y
.to_vector())
626 def is_trivial(self
):
628 Return whether or not this algebra is trivial.
630 A trivial algebra contains only the zero element.
634 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
639 sage: J = ComplexHermitianEJA(3)
645 sage: J = TrivialEJA()
650 return self
.dimension() == 0
653 def multiplication_table(self
):
655 Return a visual representation of this algebra's multiplication
656 table (on basis elements).
660 sage: from mjo.eja.eja_algebra import JordanSpinEJA
664 sage: J = JordanSpinEJA(4)
665 sage: J.multiplication_table()
666 +----++----+----+----+----+
667 | * || e0 | e1 | e2 | e3 |
668 +====++====+====+====+====+
669 | e0 || e0 | e1 | e2 | e3 |
670 +----++----+----+----+----+
671 | e1 || e1 | e0 | 0 | 0 |
672 +----++----+----+----+----+
673 | e2 || e2 | 0 | e0 | 0 |
674 +----++----+----+----+----+
675 | e3 || e3 | 0 | 0 | e0 |
676 +----++----+----+----+----+
680 # Prepend the header row.
681 M
= [["*"] + list(self
.gens())]
683 # And to each subsequent row, prepend an entry that belongs to
684 # the left-side "header column."
685 M
+= [ [self
.gens()[i
]] + [ self
.product_on_basis(i
,j
)
689 return table(M
, header_row
=True, header_column
=True, frame
=True)
692 def matrix_basis(self
):
694 Return an (often more natural) representation of this algebras
695 basis as an ordered tuple of matrices.
697 Every finite-dimensional Euclidean Jordan Algebra is a, up to
698 Jordan isomorphism, a direct sum of five simple
699 algebras---four of which comprise Hermitian matrices. And the
700 last type of algebra can of course be thought of as `n`-by-`1`
701 column matrices (ambiguusly called column vectors) to avoid
702 special cases. As a result, matrices (and column vectors) are
703 a natural representation format for Euclidean Jordan algebra
706 But, when we construct an algebra from a basis of matrices,
707 those matrix representations are lost in favor of coordinate
708 vectors *with respect to* that basis. We could eventually
709 convert back if we tried hard enough, but having the original
710 representations handy is valuable enough that we simply store
711 them and return them from this method.
713 Why implement this for non-matrix algebras? Avoiding special
714 cases for the :class:`BilinearFormEJA` pays with simplicity in
715 its own right. But mainly, we would like to be able to assume
716 that elements of a :class:`CartesianProductEJA` can be displayed
717 nicely, without having to have special classes for direct sums
718 one of whose components was a matrix algebra.
722 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
723 ....: RealSymmetricEJA)
727 sage: J = RealSymmetricEJA(2)
729 Finite family {0: e0, 1: e1, 2: e2}
730 sage: J.matrix_basis()
732 [1 0] [ 0 0.7071067811865475?] [0 0]
733 [0 0], [0.7071067811865475? 0], [0 1]
738 sage: J = JordanSpinEJA(2)
740 Finite family {0: e0, 1: e1}
741 sage: J.matrix_basis()
747 return self
._matrix
_basis
750 def matrix_space(self
):
752 Return the matrix space in which this algebra's elements live, if
753 we think of them as matrices (including column vectors of the
756 Generally this will be an `n`-by-`1` column-vector space,
757 except when the algebra is trivial. There it's `n`-by-`n`
758 (where `n` is zero), to ensure that two elements of the matrix
759 space (empty matrices) can be multiplied.
761 Matrix algebras override this with something more useful.
763 if self
.is_trivial():
764 return MatrixSpace(self
.base_ring(), 0)
766 return self
.matrix_basis()[0].parent()
772 Return the unit element of this algebra.
776 sage: from mjo.eja.eja_algebra import (HadamardEJA,
781 We can compute unit element in the Hadamard EJA::
783 sage: J = HadamardEJA(5)
785 e0 + e1 + e2 + e3 + e4
787 The unit element in the Hadamard EJA is inherited in the
788 subalgebras generated by its elements::
790 sage: J = HadamardEJA(5)
792 e0 + e1 + e2 + e3 + e4
793 sage: x = sum(J.gens())
794 sage: A = x.subalgebra_generated_by(orthonormalize=False)
797 sage: A.one().superalgebra_element()
798 e0 + e1 + e2 + e3 + e4
802 The identity element acts like the identity, regardless of
803 whether or not we orthonormalize::
805 sage: set_random_seed()
806 sage: J = random_eja()
807 sage: x = J.random_element()
808 sage: J.one()*x == x and x*J.one() == x
810 sage: A = x.subalgebra_generated_by()
811 sage: y = A.random_element()
812 sage: A.one()*y == y and y*A.one() == y
817 sage: set_random_seed()
818 sage: J = random_eja(field=QQ, orthonormalize=False)
819 sage: x = J.random_element()
820 sage: J.one()*x == x and x*J.one() == x
822 sage: A = x.subalgebra_generated_by(orthonormalize=False)
823 sage: y = A.random_element()
824 sage: A.one()*y == y and y*A.one() == y
827 The matrix of the unit element's operator is the identity,
828 regardless of the base field and whether or not we
831 sage: set_random_seed()
832 sage: J = random_eja()
833 sage: actual = J.one().operator().matrix()
834 sage: expected = matrix.identity(J.base_ring(), J.dimension())
835 sage: actual == expected
837 sage: x = J.random_element()
838 sage: A = x.subalgebra_generated_by()
839 sage: actual = A.one().operator().matrix()
840 sage: expected = matrix.identity(A.base_ring(), A.dimension())
841 sage: actual == expected
846 sage: set_random_seed()
847 sage: J = random_eja(field=QQ, orthonormalize=False)
848 sage: actual = J.one().operator().matrix()
849 sage: expected = matrix.identity(J.base_ring(), J.dimension())
850 sage: actual == expected
852 sage: x = J.random_element()
853 sage: A = x.subalgebra_generated_by(orthonormalize=False)
854 sage: actual = A.one().operator().matrix()
855 sage: expected = matrix.identity(A.base_ring(), A.dimension())
856 sage: actual == expected
859 Ensure that the cached unit element (often precomputed by
860 hand) agrees with the computed one::
862 sage: set_random_seed()
863 sage: J = random_eja()
864 sage: cached = J.one()
865 sage: J.one.clear_cache()
866 sage: J.one() == cached
871 sage: set_random_seed()
872 sage: J = random_eja(field=QQ, orthonormalize=False)
873 sage: cached = J.one()
874 sage: J.one.clear_cache()
875 sage: J.one() == cached
879 # We can brute-force compute the matrices of the operators
880 # that correspond to the basis elements of this algebra.
881 # If some linear combination of those basis elements is the
882 # algebra identity, then the same linear combination of
883 # their matrices has to be the identity matrix.
885 # Of course, matrices aren't vectors in sage, so we have to
886 # appeal to the "long vectors" isometry.
887 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
889 # Now we use basic linear algebra to find the coefficients,
890 # of the matrices-as-vectors-linear-combination, which should
891 # work for the original algebra basis too.
892 A
= matrix(self
.base_ring(), oper_vecs
)
894 # We used the isometry on the left-hand side already, but we
895 # still need to do it for the right-hand side. Recall that we
896 # wanted something that summed to the identity matrix.
897 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
899 # Now if there's an identity element in the algebra, this
900 # should work. We solve on the left to avoid having to
901 # transpose the matrix "A".
902 return self
.from_vector(A
.solve_left(b
))
905 def peirce_decomposition(self
, c
):
907 The Peirce decomposition of this algebra relative to the
910 In the future, this can be extended to a complete system of
911 orthogonal idempotents.
915 - ``c`` -- an idempotent of this algebra.
919 A triple (J0, J5, J1) containing two subalgebras and one subspace
922 - ``J0`` -- the algebra on the eigenspace of ``c.operator()``
923 corresponding to the eigenvalue zero.
925 - ``J5`` -- the eigenspace (NOT a subalgebra) of ``c.operator()``
926 corresponding to the eigenvalue one-half.
928 - ``J1`` -- the algebra on the eigenspace of ``c.operator()``
929 corresponding to the eigenvalue one.
931 These are the only possible eigenspaces for that operator, and this
932 algebra is a direct sum of them. The spaces ``J0`` and ``J1`` are
933 orthogonal, and are subalgebras of this algebra with the appropriate
938 sage: from mjo.eja.eja_algebra import random_eja, RealSymmetricEJA
942 The canonical example comes from the symmetric matrices, which
943 decompose into diagonal and off-diagonal parts::
945 sage: J = RealSymmetricEJA(3)
946 sage: C = matrix(QQ, [ [1,0,0],
950 sage: J0,J5,J1 = J.peirce_decomposition(c)
952 Euclidean Jordan algebra of dimension 1...
954 Vector space of degree 6 and dimension 2...
956 Euclidean Jordan algebra of dimension 3...
957 sage: J0.one().to_matrix()
961 sage: orig_df = AA.options.display_format
962 sage: AA.options.display_format = 'radical'
963 sage: J.from_vector(J5.basis()[0]).to_matrix()
967 sage: J.from_vector(J5.basis()[1]).to_matrix()
971 sage: AA.options.display_format = orig_df
972 sage: J1.one().to_matrix()
979 Every algebra decomposes trivially with respect to its identity
982 sage: set_random_seed()
983 sage: J = random_eja()
984 sage: J0,J5,J1 = J.peirce_decomposition(J.one())
985 sage: J0.dimension() == 0 and J5.dimension() == 0
987 sage: J1.superalgebra() == J and J1.dimension() == J.dimension()
990 The decomposition is into eigenspaces, and its components are
991 therefore necessarily orthogonal. Moreover, the identity
992 elements in the two subalgebras are the projections onto their
993 respective subspaces of the superalgebra's identity element::
995 sage: set_random_seed()
996 sage: J = random_eja()
997 sage: x = J.random_element()
998 sage: if not J.is_trivial():
999 ....: while x.is_nilpotent():
1000 ....: x = J.random_element()
1001 sage: c = x.subalgebra_idempotent()
1002 sage: J0,J5,J1 = J.peirce_decomposition(c)
1004 sage: for (w,y,z) in zip(J0.basis(), J5.basis(), J1.basis()):
1005 ....: w = w.superalgebra_element()
1006 ....: y = J.from_vector(y)
1007 ....: z = z.superalgebra_element()
1008 ....: ipsum += w.inner_product(y).abs()
1009 ....: ipsum += w.inner_product(z).abs()
1010 ....: ipsum += y.inner_product(z).abs()
1013 sage: J1(c) == J1.one()
1015 sage: J0(J.one() - c) == J0.one()
1019 if not c
.is_idempotent():
1020 raise ValueError("element is not idempotent: %s" % c
)
1022 from mjo
.eja
.eja_subalgebra
import FiniteDimensionalEJASubalgebra
1024 # Default these to what they should be if they turn out to be
1025 # trivial, because eigenspaces_left() won't return eigenvalues
1026 # corresponding to trivial spaces (e.g. it returns only the
1027 # eigenspace corresponding to lambda=1 if you take the
1028 # decomposition relative to the identity element).
1029 trivial
= FiniteDimensionalEJASubalgebra(self
, ())
1030 J0
= trivial
# eigenvalue zero
1031 J5
= VectorSpace(self
.base_ring(), 0) # eigenvalue one-half
1032 J1
= trivial
# eigenvalue one
1034 for (eigval
, eigspace
) in c
.operator().matrix().right_eigenspaces():
1035 if eigval
== ~
(self
.base_ring()(2)):
1038 gens
= tuple( self
.from_vector(b
) for b
in eigspace
.basis() )
1039 subalg
= FiniteDimensionalEJASubalgebra(self
,
1047 raise ValueError("unexpected eigenvalue: %s" % eigval
)
1052 def random_element(self
, thorough
=False):
1054 Return a random element of this algebra.
1056 Our algebra superclass method only returns a linear
1057 combination of at most two basis elements. We instead
1058 want the vector space "random element" method that
1059 returns a more diverse selection.
1063 - ``thorough`` -- (boolean; default False) whether or not we
1064 should generate irrational coefficients for the random
1065 element when our base ring is irrational; this slows the
1066 algebra operations to a crawl, but any truly random method
1070 # For a general base ring... maybe we can trust this to do the
1071 # right thing? Unlikely, but.
1072 V
= self
.vector_space()
1073 v
= V
.random_element()
1075 if self
.base_ring() is AA
:
1076 # The "random element" method of the algebraic reals is
1077 # stupid at the moment, and only returns integers between
1078 # -2 and 2, inclusive:
1080 # https://trac.sagemath.org/ticket/30875
1082 # Instead, we implement our own "random vector" method,
1083 # and then coerce that into the algebra. We use the vector
1084 # space degree here instead of the dimension because a
1085 # subalgebra could (for example) be spanned by only two
1086 # vectors, each with five coordinates. We need to
1087 # generate all five coordinates.
1089 v
*= QQbar
.random_element().real()
1091 v
*= QQ
.random_element()
1093 return self
.from_vector(V
.coordinate_vector(v
))
1095 def random_elements(self
, count
, thorough
=False):
1097 Return ``count`` random elements as a tuple.
1101 - ``thorough`` -- (boolean; default False) whether or not we
1102 should generate irrational coefficients for the random
1103 elements when our base ring is irrational; this slows the
1104 algebra operations to a crawl, but any truly random method
1109 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1113 sage: J = JordanSpinEJA(3)
1114 sage: x,y,z = J.random_elements(3)
1115 sage: all( [ x in J, y in J, z in J ])
1117 sage: len( J.random_elements(10) ) == 10
1121 return tuple( self
.random_element(thorough
)
1122 for idx
in range(count
) )
1126 def _charpoly_coefficients(self
):
1128 The `r` polynomial coefficients of the "characteristic polynomial
1133 sage: from mjo.eja.eja_algebra import random_eja
1137 The theory shows that these are all homogeneous polynomials of
1140 sage: set_random_seed()
1141 sage: J = random_eja()
1142 sage: all(p.is_homogeneous() for p in J._charpoly_coefficients())
1146 n
= self
.dimension()
1147 R
= self
.coordinate_polynomial_ring()
1149 F
= R
.fraction_field()
1152 # From a result in my book, these are the entries of the
1153 # basis representation of L_x.
1154 return sum( vars[k
]*self
.gens()[k
].operator().matrix()[i
,j
]
1157 L_x
= matrix(F
, n
, n
, L_x_i_j
)
1160 if self
.rank
.is_in_cache():
1162 # There's no need to pad the system with redundant
1163 # columns if we *know* they'll be redundant.
1166 # Compute an extra power in case the rank is equal to
1167 # the dimension (otherwise, we would stop at x^(r-1)).
1168 x_powers
= [ (L_x
**k
)*self
.one().to_vector()
1169 for k
in range(n
+1) ]
1170 A
= matrix
.column(F
, x_powers
[:n
])
1171 AE
= A
.extended_echelon_form()
1178 # The theory says that only the first "r" coefficients are
1179 # nonzero, and they actually live in the original polynomial
1180 # ring and not the fraction field. We negate them because in
1181 # the actual characteristic polynomial, they get moved to the
1182 # other side where x^r lives. We don't bother to trim A_rref
1183 # down to a square matrix and solve the resulting system,
1184 # because the upper-left r-by-r portion of A_rref is
1185 # guaranteed to be the identity matrix, so e.g.
1187 # A_rref.solve_right(Y)
1189 # would just be returning Y.
1190 return (-E
*b
)[:r
].change_ring(R
)
1195 Return the rank of this EJA.
1197 This is a cached method because we know the rank a priori for
1198 all of the algebras we can construct. Thus we can avoid the
1199 expensive ``_charpoly_coefficients()`` call unless we truly
1200 need to compute the whole characteristic polynomial.
1204 sage: from mjo.eja.eja_algebra import (HadamardEJA,
1205 ....: JordanSpinEJA,
1206 ....: RealSymmetricEJA,
1207 ....: ComplexHermitianEJA,
1208 ....: QuaternionHermitianEJA,
1213 The rank of the Jordan spin algebra is always two::
1215 sage: JordanSpinEJA(2).rank()
1217 sage: JordanSpinEJA(3).rank()
1219 sage: JordanSpinEJA(4).rank()
1222 The rank of the `n`-by-`n` Hermitian real, complex, or
1223 quaternion matrices is `n`::
1225 sage: RealSymmetricEJA(4).rank()
1227 sage: ComplexHermitianEJA(3).rank()
1229 sage: QuaternionHermitianEJA(2).rank()
1234 Ensure that every EJA that we know how to construct has a
1235 positive integer rank, unless the algebra is trivial in
1236 which case its rank will be zero::
1238 sage: set_random_seed()
1239 sage: J = random_eja()
1243 sage: r > 0 or (r == 0 and J.is_trivial())
1246 Ensure that computing the rank actually works, since the ranks
1247 of all simple algebras are known and will be cached by default::
1249 sage: set_random_seed() # long time
1250 sage: J = random_eja() # long time
1251 sage: cached = J.rank() # long time
1252 sage: J.rank.clear_cache() # long time
1253 sage: J.rank() == cached # long time
1257 return len(self
._charpoly
_coefficients
())
1260 def vector_space(self
):
1262 Return the vector space that underlies this algebra.
1266 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1270 sage: J = RealSymmetricEJA(2)
1271 sage: J.vector_space()
1272 Vector space of dimension 3 over...
1275 return self
.zero().to_vector().parent().ambient_vector_space()
1278 Element
= FiniteDimensionalEJAElement
1280 class RationalBasisEJA(FiniteDimensionalEJA
):
1282 New class for algebras whose supplied basis elements have all rational entries.
1286 sage: from mjo.eja.eja_algebra import BilinearFormEJA
1290 The supplied basis is orthonormalized by default::
1292 sage: B = matrix(QQ, [[1, 0, 0], [0, 25, -32], [0, -32, 41]])
1293 sage: J = BilinearFormEJA(B)
1294 sage: J.matrix_basis()
1311 # Abuse the check_field parameter to check that the entries of
1312 # out basis (in ambient coordinates) are in the field QQ.
1313 if not all( all(b_i
in QQ
for b_i
in b
.list()) for b
in basis
):
1314 raise TypeError("basis not rational")
1316 self
._rational
_algebra
= None
1318 # There's no point in constructing the extra algebra if this
1319 # one is already rational.
1321 # Note: the same Jordan and inner-products work here,
1322 # because they are necessarily defined with respect to
1323 # ambient coordinates and not any particular basis.
1324 self
._rational
_algebra
= FiniteDimensionalEJA(
1329 orthonormalize
=False,
1333 super().__init
__(basis
,
1337 check_field
=check_field
,
1341 def _charpoly_coefficients(self
):
1345 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
1346 ....: JordanSpinEJA)
1350 The base ring of the resulting polynomial coefficients is what
1351 it should be, and not the rationals (unless the algebra was
1352 already over the rationals)::
1354 sage: J = JordanSpinEJA(3)
1355 sage: J._charpoly_coefficients()
1356 (X1^2 - X2^2 - X3^2, -2*X1)
1357 sage: a0 = J._charpoly_coefficients()[0]
1359 Algebraic Real Field
1360 sage: a0.base_ring()
1361 Algebraic Real Field
1364 if self
._rational
_algebra
is None:
1365 # There's no need to construct *another* algebra over the
1366 # rationals if this one is already over the
1367 # rationals. Likewise, if we never orthonormalized our
1368 # basis, we might as well just use the given one.
1369 return super()._charpoly
_coefficients
()
1371 # Do the computation over the rationals. The answer will be
1372 # the same, because all we've done is a change of basis.
1373 # Then, change back from QQ to our real base ring
1374 a
= ( a_i
.change_ring(self
.base_ring())
1375 for a_i
in self
._rational
_algebra
._charpoly
_coefficients
() )
1377 if self
._deortho
_matrix
is None:
1378 # This can happen if our base ring was, say, AA and we
1379 # chose not to (or didn't need to) orthonormalize. It's
1380 # still faster to do the computations over QQ even if
1381 # the numbers in the boxes stay the same.
1384 # Otherwise, convert the coordinate variables back to the
1385 # deorthonormalized ones.
1386 R
= self
.coordinate_polynomial_ring()
1387 from sage
.modules
.free_module_element
import vector
1388 X
= vector(R
, R
.gens())
1389 BX
= self
._deortho
_matrix
*X
1391 subs_dict
= { X[i]: BX[i] for i in range(len(X)) }
1392 return tuple( a_i
.subs(subs_dict
) for a_i
in a
)
1394 class ConcreteEJA(RationalBasisEJA
):
1396 A class for the Euclidean Jordan algebras that we know by name.
1398 These are the Jordan algebras whose basis, multiplication table,
1399 rank, and so on are known a priori. More to the point, they are
1400 the Euclidean Jordan algebras for which we are able to conjure up
1401 a "random instance."
1405 sage: from mjo.eja.eja_algebra import ConcreteEJA
1409 Our basis is normalized with respect to the algebra's inner
1410 product, unless we specify otherwise::
1412 sage: set_random_seed()
1413 sage: J = ConcreteEJA.random_instance()
1414 sage: all( b.norm() == 1 for b in J.gens() )
1417 Since our basis is orthonormal with respect to the algebra's inner
1418 product, and since we know that this algebra is an EJA, any
1419 left-multiplication operator's matrix will be symmetric because
1420 natural->EJA basis representation is an isometry and within the
1421 EJA the operator is self-adjoint by the Jordan axiom::
1423 sage: set_random_seed()
1424 sage: J = ConcreteEJA.random_instance()
1425 sage: x = J.random_element()
1426 sage: x.operator().is_self_adjoint()
1431 def _max_random_instance_size():
1433 Return an integer "size" that is an upper bound on the size of
1434 this algebra when it is used in a random test
1435 case. Unfortunately, the term "size" is ambiguous -- when
1436 dealing with `R^n` under either the Hadamard or Jordan spin
1437 product, the "size" refers to the dimension `n`. When dealing
1438 with a matrix algebra (real symmetric or complex/quaternion
1439 Hermitian), it refers to the size of the matrix, which is far
1440 less than the dimension of the underlying vector space.
1442 This method must be implemented in each subclass.
1444 raise NotImplementedError
1447 def random_instance(cls
, *args
, **kwargs
):
1449 Return a random instance of this type of algebra.
1451 This method should be implemented in each subclass.
1453 from sage
.misc
.prandom
import choice
1454 eja_class
= choice(cls
.__subclasses
__())
1456 # These all bubble up to the RationalBasisEJA superclass
1457 # constructor, so any (kw)args valid there are also valid
1459 return eja_class
.random_instance(*args
, **kwargs
)
1464 def dimension_over_reals():
1466 The dimension of this matrix's base ring over the reals.
1468 The reals are dimension one over themselves, obviously; that's
1469 just `\mathbb{R}^{1}`. Likewise, the complex numbers `a + bi`
1470 have dimension two. Finally, the quaternions have dimension
1471 four over the reals.
1473 This is used to determine the size of the matrix returned from
1474 :meth:`real_embed`, among other things.
1476 raise NotImplementedError
1479 def real_embed(cls
,M
):
1481 Embed the matrix ``M`` into a space of real matrices.
1483 The matrix ``M`` can have entries in any field at the moment:
1484 the real numbers, complex numbers, or quaternions. And although
1485 they are not a field, we can probably support octonions at some
1486 point, too. This function returns a real matrix that "acts like"
1487 the original with respect to matrix multiplication; i.e.
1489 real_embed(M*N) = real_embed(M)*real_embed(N)
1492 if M
.ncols() != M
.nrows():
1493 raise ValueError("the matrix 'M' must be square")
1498 def real_unembed(cls
,M
):
1500 The inverse of :meth:`real_embed`.
1502 if M
.ncols() != M
.nrows():
1503 raise ValueError("the matrix 'M' must be square")
1504 if not ZZ(M
.nrows()).mod(cls
.dimension_over_reals()).is_zero():
1505 raise ValueError("the matrix 'M' must be a real embedding")
1509 def jordan_product(X
,Y
):
1510 return (X
*Y
+ Y
*X
)/2
1513 def trace_inner_product(cls
,X
,Y
):
1515 Compute the trace inner-product of two real-embeddings.
1519 sage: from mjo.eja.eja_algebra import (RealSymmetricEJA,
1520 ....: ComplexHermitianEJA,
1521 ....: QuaternionHermitianEJA)
1525 This gives the same answer as it would if we computed the trace
1526 from the unembedded (original) matrices::
1528 sage: set_random_seed()
1529 sage: J = RealSymmetricEJA.random_instance()
1530 sage: x,y = J.random_elements(2)
1531 sage: Xe = x.to_matrix()
1532 sage: Ye = y.to_matrix()
1533 sage: X = J.real_unembed(Xe)
1534 sage: Y = J.real_unembed(Ye)
1535 sage: expected = (X*Y).trace()
1536 sage: actual = J.trace_inner_product(Xe,Ye)
1537 sage: actual == expected
1542 sage: set_random_seed()
1543 sage: J = ComplexHermitianEJA.random_instance()
1544 sage: x,y = J.random_elements(2)
1545 sage: Xe = x.to_matrix()
1546 sage: Ye = y.to_matrix()
1547 sage: X = J.real_unembed(Xe)
1548 sage: Y = J.real_unembed(Ye)
1549 sage: expected = (X*Y).trace().real()
1550 sage: actual = J.trace_inner_product(Xe,Ye)
1551 sage: actual == expected
1556 sage: set_random_seed()
1557 sage: J = QuaternionHermitianEJA.random_instance()
1558 sage: x,y = J.random_elements(2)
1559 sage: Xe = x.to_matrix()
1560 sage: Ye = y.to_matrix()
1561 sage: X = J.real_unembed(Xe)
1562 sage: Y = J.real_unembed(Ye)
1563 sage: expected = (X*Y).trace().coefficient_tuple()[0]
1564 sage: actual = J.trace_inner_product(Xe,Ye)
1565 sage: actual == expected
1569 Xu
= cls
.real_unembed(X
)
1570 Yu
= cls
.real_unembed(Y
)
1571 tr
= (Xu
*Yu
).trace()
1574 # Works in QQ, AA, RDF, et cetera.
1576 except AttributeError:
1577 # A quaternion doesn't have a real() method, but does
1578 # have coefficient_tuple() method that returns the
1579 # coefficients of 1, i, j, and k -- in that order.
1580 return tr
.coefficient_tuple()[0]
1583 class RealMatrixEJA(MatrixEJA
):
1585 def dimension_over_reals():
1589 class RealSymmetricEJA(ConcreteEJA
, RealMatrixEJA
):
1591 The rank-n simple EJA consisting of real symmetric n-by-n
1592 matrices, the usual symmetric Jordan product, and the trace inner
1593 product. It has dimension `(n^2 + n)/2` over the reals.
1597 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1601 sage: J = RealSymmetricEJA(2)
1602 sage: e0, e1, e2 = J.gens()
1610 In theory, our "field" can be any subfield of the reals::
1612 sage: RealSymmetricEJA(2, field=RDF)
1613 Euclidean Jordan algebra of dimension 3 over Real Double Field
1614 sage: RealSymmetricEJA(2, field=RR)
1615 Euclidean Jordan algebra of dimension 3 over Real Field with
1616 53 bits of precision
1620 The dimension of this algebra is `(n^2 + n) / 2`::
1622 sage: set_random_seed()
1623 sage: n_max = RealSymmetricEJA._max_random_instance_size()
1624 sage: n = ZZ.random_element(1, n_max)
1625 sage: J = RealSymmetricEJA(n)
1626 sage: J.dimension() == (n^2 + n)/2
1629 The Jordan multiplication is what we think it is::
1631 sage: set_random_seed()
1632 sage: J = RealSymmetricEJA.random_instance()
1633 sage: x,y = J.random_elements(2)
1634 sage: actual = (x*y).to_matrix()
1635 sage: X = x.to_matrix()
1636 sage: Y = y.to_matrix()
1637 sage: expected = (X*Y + Y*X)/2
1638 sage: actual == expected
1640 sage: J(expected) == x*y
1643 We can change the generator prefix::
1645 sage: RealSymmetricEJA(3, prefix='q').gens()
1646 (q0, q1, q2, q3, q4, q5)
1648 We can construct the (trivial) algebra of rank zero::
1650 sage: RealSymmetricEJA(0)
1651 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1655 def _denormalized_basis(cls
, n
):
1657 Return a basis for the space of real symmetric n-by-n matrices.
1661 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1665 sage: set_random_seed()
1666 sage: n = ZZ.random_element(1,5)
1667 sage: B = RealSymmetricEJA._denormalized_basis(n)
1668 sage: all( M.is_symmetric() for M in B)
1672 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1676 for j
in range(i
+1):
1677 Eij
= matrix(ZZ
, n
, lambda k
,l
: k
==i
and l
==j
)
1681 Sij
= Eij
+ Eij
.transpose()
1687 def _max_random_instance_size():
1688 return 4 # Dimension 10
1691 def random_instance(cls
, **kwargs
):
1693 Return a random instance of this type of algebra.
1695 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1696 return cls(n
, **kwargs
)
1698 def __init__(self
, n
, **kwargs
):
1699 # We know this is a valid EJA, but will double-check
1700 # if the user passes check_axioms=True.
1701 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
1703 super(RealSymmetricEJA
, self
).__init
__(self
._denormalized
_basis
(n
),
1704 self
.jordan_product
,
1705 self
.trace_inner_product
,
1708 # TODO: this could be factored out somehow, but is left here
1709 # because the MatrixEJA is not presently a subclass of the
1710 # FDEJA class that defines rank() and one().
1711 self
.rank
.set_cache(n
)
1712 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
1713 self
.one
.set_cache(self(idV
))
1717 class ComplexMatrixEJA(MatrixEJA
):
1718 # A manual dictionary-cache for the complex_extension() method,
1719 # since apparently @classmethods can't also be @cached_methods.
1720 _complex_extension
= {}
1723 def complex_extension(cls
,field
):
1725 The complex field that we embed/unembed, as an extension
1726 of the given ``field``.
1728 if field
in cls
._complex
_extension
:
1729 return cls
._complex
_extension
[field
]
1731 # Sage doesn't know how to adjoin the complex "i" (the root of
1732 # x^2 + 1) to a field in a general way. Here, we just enumerate
1733 # all of the cases that I have cared to support so far.
1735 # Sage doesn't know how to embed AA into QQbar, i.e. how
1736 # to adjoin sqrt(-1) to AA.
1738 elif not field
.is_exact():
1740 F
= field
.complex_field()
1742 # Works for QQ and... maybe some other fields.
1743 R
= PolynomialRing(field
, 'z')
1745 F
= field
.extension(z
**2 + 1, 'I', embedding
=CLF(-1).sqrt())
1747 cls
._complex
_extension
[field
] = F
1751 def dimension_over_reals():
1755 def real_embed(cls
,M
):
1757 Embed the n-by-n complex matrix ``M`` into the space of real
1758 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1759 bi` to the block matrix ``[[a,b],[-b,a]]``.
1763 sage: from mjo.eja.eja_algebra import ComplexMatrixEJA
1767 sage: F = QuadraticField(-1, 'I')
1768 sage: x1 = F(4 - 2*i)
1769 sage: x2 = F(1 + 2*i)
1772 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1773 sage: ComplexMatrixEJA.real_embed(M)
1782 Embedding is a homomorphism (isomorphism, in fact)::
1784 sage: set_random_seed()
1785 sage: n = ZZ.random_element(3)
1786 sage: F = QuadraticField(-1, 'I')
1787 sage: X = random_matrix(F, n)
1788 sage: Y = random_matrix(F, n)
1789 sage: Xe = ComplexMatrixEJA.real_embed(X)
1790 sage: Ye = ComplexMatrixEJA.real_embed(Y)
1791 sage: XYe = ComplexMatrixEJA.real_embed(X*Y)
1796 super(ComplexMatrixEJA
,cls
).real_embed(M
)
1799 # We don't need any adjoined elements...
1800 field
= M
.base_ring().base_ring()
1806 blocks
.append(matrix(field
, 2, [ [ a
, b
],
1809 return matrix
.block(field
, n
, blocks
)
1813 def real_unembed(cls
,M
):
1815 The inverse of _embed_complex_matrix().
1819 sage: from mjo.eja.eja_algebra import ComplexMatrixEJA
1823 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1824 ....: [-2, 1, -4, 3],
1825 ....: [ 9, 10, 11, 12],
1826 ....: [-10, 9, -12, 11] ])
1827 sage: ComplexMatrixEJA.real_unembed(A)
1829 [ 10*I + 9 12*I + 11]
1833 Unembedding is the inverse of embedding::
1835 sage: set_random_seed()
1836 sage: F = QuadraticField(-1, 'I')
1837 sage: M = random_matrix(F, 3)
1838 sage: Me = ComplexMatrixEJA.real_embed(M)
1839 sage: ComplexMatrixEJA.real_unembed(Me) == M
1843 super(ComplexMatrixEJA
,cls
).real_unembed(M
)
1845 d
= cls
.dimension_over_reals()
1846 F
= cls
.complex_extension(M
.base_ring())
1849 # Go top-left to bottom-right (reading order), converting every
1850 # 2-by-2 block we see to a single complex element.
1852 for k
in range(n
/d
):
1853 for j
in range(n
/d
):
1854 submat
= M
[d
*k
:d
*k
+d
,d
*j
:d
*j
+d
]
1855 if submat
[0,0] != submat
[1,1]:
1856 raise ValueError('bad on-diagonal submatrix')
1857 if submat
[0,1] != -submat
[1,0]:
1858 raise ValueError('bad off-diagonal submatrix')
1859 z
= submat
[0,0] + submat
[0,1]*i
1862 return matrix(F
, n
/d
, elements
)
1865 class ComplexHermitianEJA(ConcreteEJA
, ComplexMatrixEJA
):
1867 The rank-n simple EJA consisting of complex Hermitian n-by-n
1868 matrices over the real numbers, the usual symmetric Jordan product,
1869 and the real-part-of-trace inner product. It has dimension `n^2` over
1874 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1878 In theory, our "field" can be any subfield of the reals::
1880 sage: ComplexHermitianEJA(2, field=RDF)
1881 Euclidean Jordan algebra of dimension 4 over Real Double Field
1882 sage: ComplexHermitianEJA(2, field=RR)
1883 Euclidean Jordan algebra of dimension 4 over Real Field with
1884 53 bits of precision
1888 The dimension of this algebra is `n^2`::
1890 sage: set_random_seed()
1891 sage: n_max = ComplexHermitianEJA._max_random_instance_size()
1892 sage: n = ZZ.random_element(1, n_max)
1893 sage: J = ComplexHermitianEJA(n)
1894 sage: J.dimension() == n^2
1897 The Jordan multiplication is what we think it is::
1899 sage: set_random_seed()
1900 sage: J = ComplexHermitianEJA.random_instance()
1901 sage: x,y = J.random_elements(2)
1902 sage: actual = (x*y).to_matrix()
1903 sage: X = x.to_matrix()
1904 sage: Y = y.to_matrix()
1905 sage: expected = (X*Y + Y*X)/2
1906 sage: actual == expected
1908 sage: J(expected) == x*y
1911 We can change the generator prefix::
1913 sage: ComplexHermitianEJA(2, prefix='z').gens()
1916 We can construct the (trivial) algebra of rank zero::
1918 sage: ComplexHermitianEJA(0)
1919 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1924 def _denormalized_basis(cls
, n
):
1926 Returns a basis for the space of complex Hermitian n-by-n matrices.
1928 Why do we embed these? Basically, because all of numerical linear
1929 algebra assumes that you're working with vectors consisting of `n`
1930 entries from a field and scalars from the same field. There's no way
1931 to tell SageMath that (for example) the vectors contain complex
1932 numbers, while the scalar field is real.
1936 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1940 sage: set_random_seed()
1941 sage: n = ZZ.random_element(1,5)
1942 sage: B = ComplexHermitianEJA._denormalized_basis(n)
1943 sage: all( M.is_symmetric() for M in B)
1948 R
= PolynomialRing(field
, 'z')
1950 F
= field
.extension(z
**2 + 1, 'I')
1953 # This is like the symmetric case, but we need to be careful:
1955 # * We want conjugate-symmetry, not just symmetry.
1956 # * The diagonal will (as a result) be real.
1959 Eij
= matrix
.zero(F
,n
)
1961 for j
in range(i
+1):
1965 Sij
= cls
.real_embed(Eij
)
1968 # The second one has a minus because it's conjugated.
1969 Eij
[j
,i
] = 1 # Eij = Eij + Eij.transpose()
1970 Sij_real
= cls
.real_embed(Eij
)
1972 # Eij = I*Eij - I*Eij.transpose()
1975 Sij_imag
= cls
.real_embed(Eij
)
1981 # Since we embedded these, we can drop back to the "field" that we
1982 # started with instead of the complex extension "F".
1983 return tuple( s
.change_ring(field
) for s
in S
)
1986 def __init__(self
, n
, **kwargs
):
1987 # We know this is a valid EJA, but will double-check
1988 # if the user passes check_axioms=True.
1989 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
1991 super(ComplexHermitianEJA
, self
).__init
__(self
._denormalized
_basis
(n
),
1992 self
.jordan_product
,
1993 self
.trace_inner_product
,
1995 # TODO: this could be factored out somehow, but is left here
1996 # because the MatrixEJA is not presently a subclass of the
1997 # FDEJA class that defines rank() and one().
1998 self
.rank
.set_cache(n
)
1999 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
2000 self
.one
.set_cache(self(idV
))
2003 def _max_random_instance_size():
2004 return 3 # Dimension 9
2007 def random_instance(cls
, **kwargs
):
2009 Return a random instance of this type of algebra.
2011 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2012 return cls(n
, **kwargs
)
2014 class QuaternionMatrixEJA(MatrixEJA
):
2016 # A manual dictionary-cache for the quaternion_extension() method,
2017 # since apparently @classmethods can't also be @cached_methods.
2018 _quaternion_extension
= {}
2021 def quaternion_extension(cls
,field
):
2023 The quaternion field that we embed/unembed, as an extension
2024 of the given ``field``.
2026 if field
in cls
._quaternion
_extension
:
2027 return cls
._quaternion
_extension
[field
]
2029 Q
= QuaternionAlgebra(field
,-1,-1)
2031 cls
._quaternion
_extension
[field
] = Q
2035 def dimension_over_reals():
2039 def real_embed(cls
,M
):
2041 Embed the n-by-n quaternion matrix ``M`` into the space of real
2042 matrices of size 4n-by-4n by first sending each quaternion entry `z
2043 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
2044 c+di],[-c + di, a-bi]]`, and then embedding those into a real
2049 sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA
2053 sage: Q = QuaternionAlgebra(QQ,-1,-1)
2054 sage: i,j,k = Q.gens()
2055 sage: x = 1 + 2*i + 3*j + 4*k
2056 sage: M = matrix(Q, 1, [[x]])
2057 sage: QuaternionMatrixEJA.real_embed(M)
2063 Embedding is a homomorphism (isomorphism, in fact)::
2065 sage: set_random_seed()
2066 sage: n = ZZ.random_element(2)
2067 sage: Q = QuaternionAlgebra(QQ,-1,-1)
2068 sage: X = random_matrix(Q, n)
2069 sage: Y = random_matrix(Q, n)
2070 sage: Xe = QuaternionMatrixEJA.real_embed(X)
2071 sage: Ye = QuaternionMatrixEJA.real_embed(Y)
2072 sage: XYe = QuaternionMatrixEJA.real_embed(X*Y)
2077 super(QuaternionMatrixEJA
,cls
).real_embed(M
)
2078 quaternions
= M
.base_ring()
2081 F
= QuadraticField(-1, 'I')
2086 t
= z
.coefficient_tuple()
2091 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
2092 [-c
+ d
*i
, a
- b
*i
]])
2093 realM
= ComplexMatrixEJA
.real_embed(cplxM
)
2094 blocks
.append(realM
)
2096 # We should have real entries by now, so use the realest field
2097 # we've got for the return value.
2098 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
2103 def real_unembed(cls
,M
):
2105 The inverse of _embed_quaternion_matrix().
2109 sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA
2113 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
2114 ....: [-2, 1, -4, 3],
2115 ....: [-3, 4, 1, -2],
2116 ....: [-4, -3, 2, 1]])
2117 sage: QuaternionMatrixEJA.real_unembed(M)
2118 [1 + 2*i + 3*j + 4*k]
2122 Unembedding is the inverse of embedding::
2124 sage: set_random_seed()
2125 sage: Q = QuaternionAlgebra(QQ, -1, -1)
2126 sage: M = random_matrix(Q, 3)
2127 sage: Me = QuaternionMatrixEJA.real_embed(M)
2128 sage: QuaternionMatrixEJA.real_unembed(Me) == M
2132 super(QuaternionMatrixEJA
,cls
).real_unembed(M
)
2134 d
= cls
.dimension_over_reals()
2136 # Use the base ring of the matrix to ensure that its entries can be
2137 # multiplied by elements of the quaternion algebra.
2138 Q
= cls
.quaternion_extension(M
.base_ring())
2141 # Go top-left to bottom-right (reading order), converting every
2142 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
2145 for l
in range(n
/d
):
2146 for m
in range(n
/d
):
2147 submat
= ComplexMatrixEJA
.real_unembed(
2148 M
[d
*l
:d
*l
+d
,d
*m
:d
*m
+d
] )
2149 if submat
[0,0] != submat
[1,1].conjugate():
2150 raise ValueError('bad on-diagonal submatrix')
2151 if submat
[0,1] != -submat
[1,0].conjugate():
2152 raise ValueError('bad off-diagonal submatrix')
2153 z
= submat
[0,0].real()
2154 z
+= submat
[0,0].imag()*i
2155 z
+= submat
[0,1].real()*j
2156 z
+= submat
[0,1].imag()*k
2159 return matrix(Q
, n
/d
, elements
)
2162 class QuaternionHermitianEJA(ConcreteEJA
, QuaternionMatrixEJA
):
2164 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
2165 matrices, the usual symmetric Jordan product, and the
2166 real-part-of-trace inner product. It has dimension `2n^2 - n` over
2171 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2175 In theory, our "field" can be any subfield of the reals::
2177 sage: QuaternionHermitianEJA(2, field=RDF)
2178 Euclidean Jordan algebra of dimension 6 over Real Double Field
2179 sage: QuaternionHermitianEJA(2, field=RR)
2180 Euclidean Jordan algebra of dimension 6 over Real Field with
2181 53 bits of precision
2185 The dimension of this algebra is `2*n^2 - n`::
2187 sage: set_random_seed()
2188 sage: n_max = QuaternionHermitianEJA._max_random_instance_size()
2189 sage: n = ZZ.random_element(1, n_max)
2190 sage: J = QuaternionHermitianEJA(n)
2191 sage: J.dimension() == 2*(n^2) - n
2194 The Jordan multiplication is what we think it is::
2196 sage: set_random_seed()
2197 sage: J = QuaternionHermitianEJA.random_instance()
2198 sage: x,y = J.random_elements(2)
2199 sage: actual = (x*y).to_matrix()
2200 sage: X = x.to_matrix()
2201 sage: Y = y.to_matrix()
2202 sage: expected = (X*Y + Y*X)/2
2203 sage: actual == expected
2205 sage: J(expected) == x*y
2208 We can change the generator prefix::
2210 sage: QuaternionHermitianEJA(2, prefix='a').gens()
2211 (a0, a1, a2, a3, a4, a5)
2213 We can construct the (trivial) algebra of rank zero::
2215 sage: QuaternionHermitianEJA(0)
2216 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2220 def _denormalized_basis(cls
, n
):
2222 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
2224 Why do we embed these? Basically, because all of numerical
2225 linear algebra assumes that you're working with vectors consisting
2226 of `n` entries from a field and scalars from the same field. There's
2227 no way to tell SageMath that (for example) the vectors contain
2228 complex numbers, while the scalar field is real.
2232 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2236 sage: set_random_seed()
2237 sage: n = ZZ.random_element(1,5)
2238 sage: B = QuaternionHermitianEJA._denormalized_basis(n)
2239 sage: all( M.is_symmetric() for M in B )
2244 Q
= QuaternionAlgebra(QQ
,-1,-1)
2247 # This is like the symmetric case, but we need to be careful:
2249 # * We want conjugate-symmetry, not just symmetry.
2250 # * The diagonal will (as a result) be real.
2253 Eij
= matrix
.zero(Q
,n
)
2255 for j
in range(i
+1):
2259 Sij
= cls
.real_embed(Eij
)
2262 # The second, third, and fourth ones have a minus
2263 # because they're conjugated.
2264 # Eij = Eij + Eij.transpose()
2266 Sij_real
= cls
.real_embed(Eij
)
2268 # Eij = I*(Eij - Eij.transpose())
2271 Sij_I
= cls
.real_embed(Eij
)
2273 # Eij = J*(Eij - Eij.transpose())
2276 Sij_J
= cls
.real_embed(Eij
)
2278 # Eij = K*(Eij - Eij.transpose())
2281 Sij_K
= cls
.real_embed(Eij
)
2287 # Since we embedded these, we can drop back to the "field" that we
2288 # started with instead of the quaternion algebra "Q".
2289 return tuple( s
.change_ring(field
) for s
in S
)
2292 def __init__(self
, n
, **kwargs
):
2293 # We know this is a valid EJA, but will double-check
2294 # if the user passes check_axioms=True.
2295 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2297 super(QuaternionHermitianEJA
, self
).__init
__(self
._denormalized
_basis
(n
),
2298 self
.jordan_product
,
2299 self
.trace_inner_product
,
2301 # TODO: this could be factored out somehow, but is left here
2302 # because the MatrixEJA is not presently a subclass of the
2303 # FDEJA class that defines rank() and one().
2304 self
.rank
.set_cache(n
)
2305 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
2306 self
.one
.set_cache(self(idV
))
2310 def _max_random_instance_size():
2312 The maximum rank of a random QuaternionHermitianEJA.
2314 return 2 # Dimension 6
2317 def random_instance(cls
, **kwargs
):
2319 Return a random instance of this type of algebra.
2321 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2322 return cls(n
, **kwargs
)
2325 class HadamardEJA(ConcreteEJA
):
2327 Return the Euclidean Jordan Algebra corresponding to the set
2328 `R^n` under the Hadamard product.
2330 Note: this is nothing more than the Cartesian product of ``n``
2331 copies of the spin algebra. Once Cartesian product algebras
2332 are implemented, this can go.
2336 sage: from mjo.eja.eja_algebra import HadamardEJA
2340 This multiplication table can be verified by hand::
2342 sage: J = HadamardEJA(3)
2343 sage: e0,e1,e2 = J.gens()
2359 We can change the generator prefix::
2361 sage: HadamardEJA(3, prefix='r').gens()
2365 def __init__(self
, n
, **kwargs
):
2367 jordan_product
= lambda x
,y
: x
2368 inner_product
= lambda x
,y
: x
2370 def jordan_product(x
,y
):
2372 return P( xi
*yi
for (xi
,yi
) in zip(x
,y
) )
2374 def inner_product(x
,y
):
2377 # New defaults for keyword arguments. Don't orthonormalize
2378 # because our basis is already orthonormal with respect to our
2379 # inner-product. Don't check the axioms, because we know this
2380 # is a valid EJA... but do double-check if the user passes
2381 # check_axioms=True. Note: we DON'T override the "check_field"
2382 # default here, because the user can pass in a field!
2383 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2384 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2386 column_basis
= tuple( b
.column() for b
in FreeModule(ZZ
, n
).basis() )
2387 super().__init
__(column_basis
, jordan_product
, inner_product
, **kwargs
)
2388 self
.rank
.set_cache(n
)
2391 self
.one
.set_cache( self
.zero() )
2393 self
.one
.set_cache( sum(self
.gens()) )
2396 def _max_random_instance_size():
2398 The maximum dimension of a random HadamardEJA.
2403 def random_instance(cls
, **kwargs
):
2405 Return a random instance of this type of algebra.
2407 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2408 return cls(n
, **kwargs
)
2411 class BilinearFormEJA(ConcreteEJA
):
2413 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2414 with the half-trace inner product and jordan product ``x*y =
2415 (<Bx,y>,y_bar>, x0*y_bar + y0*x_bar)`` where `B = 1 \times B22` is
2416 a symmetric positive-definite "bilinear form" matrix. Its
2417 dimension is the size of `B`, and it has rank two in dimensions
2418 larger than two. It reduces to the ``JordanSpinEJA`` when `B` is
2419 the identity matrix of order ``n``.
2421 We insist that the one-by-one upper-left identity block of `B` be
2422 passed in as well so that we can be passed a matrix of size zero
2423 to construct a trivial algebra.
2427 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
2428 ....: JordanSpinEJA)
2432 When no bilinear form is specified, the identity matrix is used,
2433 and the resulting algebra is the Jordan spin algebra::
2435 sage: B = matrix.identity(AA,3)
2436 sage: J0 = BilinearFormEJA(B)
2437 sage: J1 = JordanSpinEJA(3)
2438 sage: J0.multiplication_table() == J0.multiplication_table()
2441 An error is raised if the matrix `B` does not correspond to a
2442 positive-definite bilinear form::
2444 sage: B = matrix.random(QQ,2,3)
2445 sage: J = BilinearFormEJA(B)
2446 Traceback (most recent call last):
2448 ValueError: bilinear form is not positive-definite
2449 sage: B = matrix.zero(QQ,3)
2450 sage: J = BilinearFormEJA(B)
2451 Traceback (most recent call last):
2453 ValueError: bilinear form is not positive-definite
2457 We can create a zero-dimensional algebra::
2459 sage: B = matrix.identity(AA,0)
2460 sage: J = BilinearFormEJA(B)
2464 We can check the multiplication condition given in the Jordan, von
2465 Neumann, and Wigner paper (and also discussed on my "On the
2466 symmetry..." paper). Note that this relies heavily on the standard
2467 choice of basis, as does anything utilizing the bilinear form
2468 matrix. We opt not to orthonormalize the basis, because if we
2469 did, we would have to normalize the `s_{i}` in a similar manner::
2471 sage: set_random_seed()
2472 sage: n = ZZ.random_element(5)
2473 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2474 sage: B11 = matrix.identity(QQ,1)
2475 sage: B22 = M.transpose()*M
2476 sage: B = block_matrix(2,2,[ [B11,0 ],
2478 sage: J = BilinearFormEJA(B, orthonormalize=False)
2479 sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
2480 sage: V = J.vector_space()
2481 sage: sis = [ J( V([0] + (M.inverse()*ei).list()).column() )
2482 ....: for ei in eis ]
2483 sage: actual = [ sis[i]*sis[j]
2484 ....: for i in range(n-1)
2485 ....: for j in range(n-1) ]
2486 sage: expected = [ J.one() if i == j else J.zero()
2487 ....: for i in range(n-1)
2488 ....: for j in range(n-1) ]
2489 sage: actual == expected
2493 def __init__(self
, B
, **kwargs
):
2494 # The matrix "B" is supplied by the user in most cases,
2495 # so it makes sense to check whether or not its positive-
2496 # definite unless we are specifically asked not to...
2497 if ("check_axioms" not in kwargs
) or kwargs
["check_axioms"]:
2498 if not B
.is_positive_definite():
2499 raise ValueError("bilinear form is not positive-definite")
2501 # However, all of the other data for this EJA is computed
2502 # by us in manner that guarantees the axioms are
2503 # satisfied. So, again, unless we are specifically asked to
2504 # verify things, we'll skip the rest of the checks.
2505 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2507 def inner_product(x
,y
):
2508 return (y
.T
*B
*x
)[0,0]
2510 def jordan_product(x
,y
):
2516 z0
= inner_product(y
,x
)
2517 zbar
= y0
*xbar
+ x0
*ybar
2518 return P([z0
] + zbar
.list())
2521 column_basis
= tuple( b
.column() for b
in FreeModule(ZZ
, n
).basis() )
2522 super(BilinearFormEJA
, self
).__init
__(column_basis
,
2527 # The rank of this algebra is two, unless we're in a
2528 # one-dimensional ambient space (because the rank is bounded
2529 # by the ambient dimension).
2530 self
.rank
.set_cache(min(n
,2))
2533 self
.one
.set_cache( self
.zero() )
2535 self
.one
.set_cache( self
.monomial(0) )
2538 def _max_random_instance_size():
2540 The maximum dimension of a random BilinearFormEJA.
2545 def random_instance(cls
, **kwargs
):
2547 Return a random instance of this algebra.
2549 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2551 B
= matrix
.identity(ZZ
, n
)
2552 return cls(B
, **kwargs
)
2554 B11
= matrix
.identity(ZZ
, 1)
2555 M
= matrix
.random(ZZ
, n
-1)
2556 I
= matrix
.identity(ZZ
, n
-1)
2558 while alpha
.is_zero():
2559 alpha
= ZZ
.random_element().abs()
2560 B22
= M
.transpose()*M
+ alpha
*I
2562 from sage
.matrix
.special
import block_matrix
2563 B
= block_matrix(2,2, [ [B11
, ZZ(0) ],
2566 return cls(B
, **kwargs
)
2569 class JordanSpinEJA(BilinearFormEJA
):
2571 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2572 with the usual inner product and jordan product ``x*y =
2573 (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2578 sage: from mjo.eja.eja_algebra import JordanSpinEJA
2582 This multiplication table can be verified by hand::
2584 sage: J = JordanSpinEJA(4)
2585 sage: e0,e1,e2,e3 = J.gens()
2601 We can change the generator prefix::
2603 sage: JordanSpinEJA(2, prefix='B').gens()
2608 Ensure that we have the usual inner product on `R^n`::
2610 sage: set_random_seed()
2611 sage: J = JordanSpinEJA.random_instance()
2612 sage: x,y = J.random_elements(2)
2613 sage: actual = x.inner_product(y)
2614 sage: expected = x.to_vector().inner_product(y.to_vector())
2615 sage: actual == expected
2619 def __init__(self
, n
, **kwargs
):
2620 # This is a special case of the BilinearFormEJA with the
2621 # identity matrix as its bilinear form.
2622 B
= matrix
.identity(ZZ
, n
)
2624 # Don't orthonormalize because our basis is already
2625 # orthonormal with respect to our inner-product.
2626 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2628 # But also don't pass check_field=False here, because the user
2629 # can pass in a field!
2630 super(JordanSpinEJA
, self
).__init
__(B
, **kwargs
)
2633 def _max_random_instance_size():
2635 The maximum dimension of a random JordanSpinEJA.
2640 def random_instance(cls
, **kwargs
):
2642 Return a random instance of this type of algebra.
2644 Needed here to override the implementation for ``BilinearFormEJA``.
2646 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2647 return cls(n
, **kwargs
)
2650 class TrivialEJA(ConcreteEJA
):
2652 The trivial Euclidean Jordan algebra consisting of only a zero element.
2656 sage: from mjo.eja.eja_algebra import TrivialEJA
2660 sage: J = TrivialEJA()
2667 sage: 7*J.one()*12*J.one()
2669 sage: J.one().inner_product(J.one())
2671 sage: J.one().norm()
2673 sage: J.one().subalgebra_generated_by()
2674 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2679 def __init__(self
, **kwargs
):
2680 jordan_product
= lambda x
,y
: x
2681 inner_product
= lambda x
,y
: 0
2684 # New defaults for keyword arguments
2685 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2686 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2688 super(TrivialEJA
, self
).__init
__(basis
,
2692 # The rank is zero using my definition, namely the dimension of the
2693 # largest subalgebra generated by any element.
2694 self
.rank
.set_cache(0)
2695 self
.one
.set_cache( self
.zero() )
2698 def random_instance(cls
, **kwargs
):
2699 # We don't take a "size" argument so the superclass method is
2700 # inappropriate for us.
2701 return cls(**kwargs
)
2704 class CartesianProductEJA(CombinatorialFreeModule_CartesianProduct
,
2705 FiniteDimensionalEJA
):
2707 The external (orthogonal) direct sum of two or more Euclidean
2708 Jordan algebras. Every Euclidean Jordan algebra decomposes into an
2709 orthogonal direct sum of simple Euclidean Jordan algebras which is
2710 then isometric to a Cartesian product, so no generality is lost by
2711 providing only this construction.
2715 sage: from mjo.eja.eja_algebra import (random_eja,
2716 ....: CartesianProductEJA,
2718 ....: JordanSpinEJA,
2719 ....: RealSymmetricEJA)
2723 The Jordan product is inherited from our factors and implemented by
2724 our CombinatorialFreeModule Cartesian product superclass::
2726 sage: set_random_seed()
2727 sage: J1 = HadamardEJA(2)
2728 sage: J2 = RealSymmetricEJA(2)
2729 sage: J = cartesian_product([J1,J2])
2730 sage: x,y = J.random_elements(2)
2734 The ability to retrieve the original factors is implemented by our
2735 CombinatorialFreeModule Cartesian product superclass::
2737 sage: J1 = HadamardEJA(2, field=QQ)
2738 sage: J2 = JordanSpinEJA(3, field=QQ)
2739 sage: J = cartesian_product([J1,J2])
2740 sage: J.cartesian_factors()
2741 (Euclidean Jordan algebra of dimension 2 over Rational Field,
2742 Euclidean Jordan algebra of dimension 3 over Rational Field)
2744 You can provide more than two factors::
2746 sage: J1 = HadamardEJA(2)
2747 sage: J2 = JordanSpinEJA(3)
2748 sage: J3 = RealSymmetricEJA(3)
2749 sage: cartesian_product([J1,J2,J3])
2750 Euclidean Jordan algebra of dimension 2 over Algebraic Real
2751 Field (+) Euclidean Jordan algebra of dimension 3 over Algebraic
2752 Real Field (+) Euclidean Jordan algebra of dimension 6 over
2753 Algebraic Real Field
2755 Rank is additive on a Cartesian product::
2757 sage: J1 = HadamardEJA(1)
2758 sage: J2 = RealSymmetricEJA(2)
2759 sage: J = cartesian_product([J1,J2])
2760 sage: J1.rank.clear_cache()
2761 sage: J2.rank.clear_cache()
2762 sage: J.rank.clear_cache()
2765 sage: J.rank() == J1.rank() + J2.rank()
2768 The same rank computation works over the rationals, with whatever
2771 sage: J1 = HadamardEJA(1, field=QQ, orthonormalize=False)
2772 sage: J2 = RealSymmetricEJA(2, field=QQ, orthonormalize=False)
2773 sage: J = cartesian_product([J1,J2])
2774 sage: J1.rank.clear_cache()
2775 sage: J2.rank.clear_cache()
2776 sage: J.rank.clear_cache()
2779 sage: J.rank() == J1.rank() + J2.rank()
2784 All factors must share the same base field::
2786 sage: J1 = HadamardEJA(2, field=QQ)
2787 sage: J2 = RealSymmetricEJA(2)
2788 sage: CartesianProductEJA((J1,J2))
2789 Traceback (most recent call last):
2791 ValueError: all factors must share the same base field
2793 The "cached" Jordan and inner products are the componentwise
2796 sage: set_random_seed()
2797 sage: J1 = random_eja()
2798 sage: J2 = random_eja()
2799 sage: J = cartesian_product([J1,J2])
2800 sage: x,y = J.random_elements(2)
2801 sage: x*y == J.cartesian_jordan_product(x,y)
2803 sage: x.inner_product(y) == J.cartesian_inner_product(x,y)
2806 The cached unit element is the same one that would be computed::
2808 sage: set_random_seed() # long time
2809 sage: J1 = random_eja() # long time
2810 sage: J2 = random_eja() # long time
2811 sage: J = cartesian_product([J1,J2]) # long time
2812 sage: actual = J.one() # long time
2813 sage: J.one.clear_cache() # long time
2814 sage: expected = J.one() # long time
2815 sage: actual == expected # long time
2819 def __init__(self
, modules
, **kwargs
):
2820 CombinatorialFreeModule_CartesianProduct
.__init
__(self
,
2823 field
= modules
[0].base_ring()
2824 if not all( J
.base_ring() == field
for J
in modules
):
2825 raise ValueError("all factors must share the same base field")
2827 # The definition of matrix_space() and self.basis() relies
2828 # only on the stuff in the CFM_CartesianProduct class, which
2829 # we've already initialized.
2830 Js
= self
.cartesian_factors()
2832 MS
= self
.matrix_space()
2834 MS(tuple( self
.cartesian_projection(i
)(b
).to_matrix()
2835 for i
in range(m
) ))
2836 for b
in self
.basis()
2839 # Define jordan/inner products that operate on that matrix_basis.
2840 def jordan_product(x
,y
):
2842 (Js
[i
](x
[i
])*Js
[i
](y
[i
])).to_matrix() for i
in range(m
)
2845 def inner_product(x
, y
):
2847 Js
[i
](x
[i
]).inner_product(Js
[i
](y
[i
])) for i
in range(m
)
2850 # There's no need to check the field since it already came
2851 # from an EJA. Likewise the axioms are guaranteed to be
2852 # satisfied, unless the guy writing this class sucks.
2854 # If you want the basis to be orthonormalized, orthonormalize
2856 FiniteDimensionalEJA
.__init
__(self
,
2861 orthonormalize
=False,
2862 cartesian_product
=True,
2866 ones
= tuple(J
.one() for J
in modules
)
2867 self
.one
.set_cache(self
._cartesian
_product
_of
_elements
(ones
))
2868 self
.rank
.set_cache(sum(J
.rank() for J
in modules
))
2870 def matrix_space(self
):
2872 Return the space that our matrix basis lives in as a Cartesian
2877 sage: from mjo.eja.eja_algebra import (HadamardEJA,
2878 ....: RealSymmetricEJA)
2882 sage: J1 = HadamardEJA(1)
2883 sage: J2 = RealSymmetricEJA(2)
2884 sage: J = cartesian_product([J1,J2])
2885 sage: J.matrix_space()
2886 The Cartesian product of (Full MatrixSpace of 1 by 1 dense
2887 matrices over Algebraic Real Field, Full MatrixSpace of 2
2888 by 2 dense matrices over Algebraic Real Field)
2891 from sage
.categories
.cartesian_product
import cartesian_product
2892 return cartesian_product( [J
.matrix_space()
2893 for J
in self
.cartesian_factors()] )
2896 def cartesian_projection(self
, i
):
2900 sage: from mjo.eja.eja_algebra import (random_eja,
2901 ....: JordanSpinEJA,
2903 ....: RealSymmetricEJA,
2904 ....: ComplexHermitianEJA)
2908 The projection morphisms are Euclidean Jordan algebra
2911 sage: J1 = HadamardEJA(2)
2912 sage: J2 = RealSymmetricEJA(2)
2913 sage: J = cartesian_product([J1,J2])
2914 sage: J.cartesian_projection(0)
2915 Linear operator between finite-dimensional Euclidean Jordan
2916 algebras represented by the matrix:
2919 Domain: Euclidean Jordan algebra of dimension 2 over Algebraic
2920 Real Field (+) Euclidean Jordan algebra of dimension 3 over
2921 Algebraic Real Field
2922 Codomain: Euclidean Jordan algebra of dimension 2 over Algebraic
2924 sage: J.cartesian_projection(1)
2925 Linear operator between finite-dimensional Euclidean Jordan
2926 algebras represented by the matrix:
2930 Domain: Euclidean Jordan algebra of dimension 2 over Algebraic
2931 Real Field (+) Euclidean Jordan algebra of dimension 3 over
2932 Algebraic Real Field
2933 Codomain: Euclidean Jordan algebra of dimension 3 over Algebraic
2936 The projections work the way you'd expect on the vector
2937 representation of an element::
2939 sage: J1 = JordanSpinEJA(2)
2940 sage: J2 = ComplexHermitianEJA(2)
2941 sage: J = cartesian_product([J1,J2])
2942 sage: pi_left = J.cartesian_projection(0)
2943 sage: pi_right = J.cartesian_projection(1)
2944 sage: pi_left(J.one()).to_vector()
2946 sage: pi_right(J.one()).to_vector()
2948 sage: J.one().to_vector()
2953 The answer never changes::
2955 sage: set_random_seed()
2956 sage: J1 = random_eja()
2957 sage: J2 = random_eja()
2958 sage: J = cartesian_product([J1,J2])
2959 sage: P0 = J.cartesian_projection(0)
2960 sage: P1 = J.cartesian_projection(0)
2965 Ji
= self
.cartesian_factors()[i
]
2966 # Requires the fix on Trac 31421/31422 to work!
2967 Pi
= super().cartesian_projection(i
)
2968 return FiniteDimensionalEJAOperator(self
,Ji
,Pi
.matrix())
2971 def cartesian_embedding(self
, i
):
2975 sage: from mjo.eja.eja_algebra import (random_eja,
2976 ....: JordanSpinEJA,
2978 ....: RealSymmetricEJA)
2982 The embedding morphisms are Euclidean Jordan algebra
2985 sage: J1 = HadamardEJA(2)
2986 sage: J2 = RealSymmetricEJA(2)
2987 sage: J = cartesian_product([J1,J2])
2988 sage: J.cartesian_embedding(0)
2989 Linear operator between finite-dimensional Euclidean Jordan
2990 algebras represented by the matrix:
2996 Domain: Euclidean Jordan algebra of dimension 2 over
2997 Algebraic Real Field
2998 Codomain: Euclidean Jordan algebra of dimension 2 over
2999 Algebraic Real Field (+) Euclidean Jordan algebra of
3000 dimension 3 over Algebraic Real Field
3001 sage: J.cartesian_embedding(1)
3002 Linear operator between finite-dimensional Euclidean Jordan
3003 algebras represented by the matrix:
3009 Domain: Euclidean Jordan algebra of dimension 3 over
3010 Algebraic Real Field
3011 Codomain: Euclidean Jordan algebra of dimension 2 over
3012 Algebraic Real Field (+) Euclidean Jordan algebra of
3013 dimension 3 over Algebraic Real Field
3015 The embeddings work the way you'd expect on the vector
3016 representation of an element::
3018 sage: J1 = JordanSpinEJA(3)
3019 sage: J2 = RealSymmetricEJA(2)
3020 sage: J = cartesian_product([J1,J2])
3021 sage: iota_left = J.cartesian_embedding(0)
3022 sage: iota_right = J.cartesian_embedding(1)
3023 sage: iota_left(J1.zero()) == J.zero()
3025 sage: iota_right(J2.zero()) == J.zero()
3027 sage: J1.one().to_vector()
3029 sage: iota_left(J1.one()).to_vector()
3031 sage: J2.one().to_vector()
3033 sage: iota_right(J2.one()).to_vector()
3035 sage: J.one().to_vector()
3040 The answer never changes::
3042 sage: set_random_seed()
3043 sage: J1 = random_eja()
3044 sage: J2 = random_eja()
3045 sage: J = cartesian_product([J1,J2])
3046 sage: E0 = J.cartesian_embedding(0)
3047 sage: E1 = J.cartesian_embedding(0)
3051 Composing a projection with the corresponding inclusion should
3052 produce the identity map, and mismatching them should produce
3055 sage: set_random_seed()
3056 sage: J1 = random_eja()
3057 sage: J2 = random_eja()
3058 sage: J = cartesian_product([J1,J2])
3059 sage: iota_left = J.cartesian_embedding(0)
3060 sage: iota_right = J.cartesian_embedding(1)
3061 sage: pi_left = J.cartesian_projection(0)
3062 sage: pi_right = J.cartesian_projection(1)
3063 sage: pi_left*iota_left == J1.one().operator()
3065 sage: pi_right*iota_right == J2.one().operator()
3067 sage: (pi_left*iota_right).is_zero()
3069 sage: (pi_right*iota_left).is_zero()
3073 Ji
= self
.cartesian_factors()[i
]
3074 # Requires the fix on Trac 31421/31422 to work!
3075 Ei
= super().cartesian_embedding(i
)
3076 return FiniteDimensionalEJAOperator(Ji
,self
,Ei
.matrix())
3079 def cartesian_jordan_product(self
, x
, y
):
3081 The componentwise Jordan product.
3083 We project ``x`` and ``y`` onto our factors, and add up the
3084 Jordan products from the subalgebras. This may still be useful
3085 after (if) the default Jordan product in the Cartesian product
3086 algebra is overridden.
3090 sage: from mjo.eja.eja_algebra import (HadamardEJA,
3091 ....: JordanSpinEJA)
3095 sage: J1 = HadamardEJA(3)
3096 sage: J2 = JordanSpinEJA(3)
3097 sage: J = cartesian_product([J1,J2])
3098 sage: x1 = J1.from_vector(vector(QQ,(1,2,1)))
3099 sage: y1 = J1.from_vector(vector(QQ,(1,0,2)))
3100 sage: x2 = J2.from_vector(vector(QQ,(1,2,3)))
3101 sage: y2 = J2.from_vector(vector(QQ,(1,1,1)))
3102 sage: z1 = J.from_vector(vector(QQ,(1,2,1,1,2,3)))
3103 sage: z2 = J.from_vector(vector(QQ,(1,0,2,1,1,1)))
3104 sage: (x1*y1).to_vector()
3106 sage: (x2*y2).to_vector()
3108 sage: J.cartesian_jordan_product(z1,z2).to_vector()
3112 m
= len(self
.cartesian_factors())
3113 projections
= ( self
.cartesian_projection(i
) for i
in range(m
) )
3114 products
= ( P(x
)*P(y
) for P
in projections
)
3115 return self
._cartesian
_product
_of
_elements
(tuple(products
))
3117 def cartesian_inner_product(self
, x
, y
):
3119 The standard componentwise Cartesian inner-product.
3121 We project ``x`` and ``y`` onto our factors, and add up the
3122 inner-products from the subalgebras. This may still be useful
3123 after (if) the default inner product in the Cartesian product
3124 algebra is overridden.
3128 sage: from mjo.eja.eja_algebra import (HadamardEJA,
3129 ....: QuaternionHermitianEJA)
3133 sage: J1 = HadamardEJA(3,field=QQ)
3134 sage: J2 = QuaternionHermitianEJA(2,field=QQ,orthonormalize=False)
3135 sage: J = cartesian_product([J1,J2])
3140 sage: x1.inner_product(x2)
3142 sage: y1.inner_product(y2)
3144 sage: z1 = J._cartesian_product_of_elements((x1,y1))
3145 sage: z2 = J._cartesian_product_of_elements((x2,y2))
3146 sage: J.cartesian_inner_product(z1,z2)
3150 m
= len(self
.cartesian_factors())
3151 projections
= ( self
.cartesian_projection(i
) for i
in range(m
) )
3152 return sum( P(x
).inner_product(P(y
)) for P
in projections
)
3155 def _element_constructor_(self
, elt
):
3157 Construct an element of this algebra from an ordered tuple.
3159 We just apply the element constructor from each of our factors
3160 to the corresponding component of the tuple, and package up
3165 sage: from mjo.eja.eja_algebra import (HadamardEJA,
3166 ....: RealSymmetricEJA)
3170 sage: J1 = HadamardEJA(3)
3171 sage: J2 = RealSymmetricEJA(2)
3172 sage: J = cartesian_product([J1,J2])
3173 sage: J( (J1.matrix_basis()[1], J2.matrix_basis()[2]) )
3176 m
= len(self
.cartesian_factors())
3178 z
= tuple( self
.cartesian_factors()[i
](elt
[i
]) for i
in range(m
) )
3179 return self
._cartesian
_product
_of
_elements
(z
)
3181 raise ValueError("not an element of this algebra")
3183 Element
= CartesianProductEJAElement
3186 FiniteDimensionalEJA
.CartesianProduct
= CartesianProductEJA
3188 random_eja
= ConcreteEJA
.random_instance
3189 #def random_eja(*args, **kwargs):
3190 # from sage.categories.cartesian_product import cartesian_product
3191 # J1 = HadamardEJA(1, **kwargs)
3192 # J2 = RealSymmetricEJA(2, **kwargs)
3193 # J = cartesian_product([J1,J2])