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
.combinat
.free_module
import CombinatorialFreeModule
24 from sage
.matrix
.constructor
import matrix
25 from sage
.matrix
.matrix_space
import MatrixSpace
26 from sage
.misc
.cachefunc
import cached_method
27 from sage
.misc
.table
import table
28 from sage
.modules
.free_module
import FreeModule
, VectorSpace
29 from sage
.rings
.all
import (ZZ
, QQ
, AA
, QQbar
, RR
, RLF
, CLF
,
32 from mjo
.eja
.eja_element
import FiniteDimensionalEJAElement
33 from mjo
.eja
.eja_operator
import FiniteDimensionalEJAOperator
34 from mjo
.eja
.eja_utils
import _mat2vec
36 class FiniteDimensionalEJA(CombinatorialFreeModule
):
38 A finite-dimensional Euclidean Jordan algebra.
42 - basis -- a tuple of basis elements in their matrix form.
44 - jordan_product -- function of two elements (in matrix form)
45 that returns their jordan product in this algebra; this will
46 be applied to ``basis`` to compute a multiplication table for
49 - inner_product -- function of two elements (in matrix form) that
50 returns their inner product. This will be applied to ``basis`` to
51 compute an inner-product table (basically a matrix) for this algebra.
54 Element
= FiniteDimensionalEJAElement
68 if not field
.is_subring(RR
):
69 # Note: this does return true for the real algebraic
70 # field, the rationals, and any quadratic field where
71 # we've specified a real embedding.
72 raise ValueError("scalar field is not real")
74 # If the basis given to us wasn't over the field that it's
75 # supposed to be over, fix that. Or, you know, crash.
76 basis
= tuple( b
.change_ring(field
) for b
in basis
)
79 # Check commutativity of the Jordan and inner-products.
80 # This has to be done before we build the multiplication
81 # and inner-product tables/matrices, because we take
82 # advantage of symmetry in the process.
83 if not all( jordan_product(bi
,bj
) == jordan_product(bj
,bi
)
86 raise ValueError("Jordan product is not commutative")
88 if not all( inner_product(bi
,bj
) == inner_product(bj
,bi
)
91 raise ValueError("inner-product is not commutative")
94 category
= MagmaticAlgebras(field
).FiniteDimensional()
95 category
= category
.WithBasis().Unital()
97 # Element subalgebras can take advantage of this.
98 category
= category
.Associative()
100 # Call the superclass constructor so that we can use its from_vector()
101 # method to build our multiplication table.
103 super().__init
__(field
,
109 # Now comes all of the hard work. We'll be constructing an
110 # ambient vector space V that our (vectorized) basis lives in,
111 # as well as a subspace W of V spanned by those (vectorized)
112 # basis elements. The W-coordinates are the coefficients that
113 # we see in things like x = 1*e1 + 2*e2.
118 # Works on both column and square matrices...
119 degree
= len(basis
[0].list())
121 # Build an ambient space that fits our matrix basis when
122 # written out as "long vectors."
123 V
= VectorSpace(field
, degree
)
125 # The matrix that will hole the orthonormal -> unorthonormal
126 # coordinate transformation.
127 self
._deortho
_matrix
= None
130 # Save a copy of the un-orthonormalized basis for later.
131 # Convert it to ambient V (vector) coordinates while we're
132 # at it, because we'd have to do it later anyway.
133 deortho_vector_basis
= tuple( V(b
.list()) for b
in basis
)
135 from mjo
.eja
.eja_utils
import gram_schmidt
136 basis
= tuple(gram_schmidt(basis
, inner_product
))
138 # Save the (possibly orthonormalized) matrix basis for
140 self
._matrix
_basis
= basis
142 # Now create the vector space for the algebra, which will have
143 # its own set of non-ambient coordinates (in terms of the
145 vector_basis
= tuple( V(b
.list()) for b
in basis
)
146 W
= V
.span_of_basis( vector_basis
, check
=check_axioms
)
149 # Now "W" is the vector space of our algebra coordinates. The
150 # variables "X1", "X2",... refer to the entries of vectors in
151 # W. Thus to convert back and forth between the orthonormal
152 # coordinates and the given ones, we need to stick the original
154 U
= V
.span_of_basis( deortho_vector_basis
, check
=check_axioms
)
155 self
._deortho
_matrix
= matrix( U
.coordinate_vector(q
)
156 for q
in vector_basis
)
159 # Now we actually compute the multiplication and inner-product
160 # tables/matrices using the possibly-orthonormalized basis.
161 self
._inner
_product
_matrix
= matrix
.identity(field
, n
)
162 self
._multiplication
_table
= [ [0 for j
in range(i
+1)]
165 # Note: the Jordan and inner-products are defined in terms
166 # of the ambient basis. It's important that their arguments
167 # are in ambient coordinates as well.
170 # ortho basis w.r.t. ambient coords
174 # The jordan product returns a matrixy answer, so we
175 # have to convert it to the algebra coordinates.
176 elt
= jordan_product(q_i
, q_j
)
177 elt
= W
.coordinate_vector(V(elt
.list()))
178 self
._multiplication
_table
[i
][j
] = self
.from_vector(elt
)
180 if not orthonormalize
:
181 # If we're orthonormalizing the basis with respect
182 # to an inner-product, then the inner-product
183 # matrix with respect to the resulting basis is
184 # just going to be the identity.
185 ip
= inner_product(q_i
, q_j
)
186 self
._inner
_product
_matrix
[i
,j
] = ip
187 self
._inner
_product
_matrix
[j
,i
] = ip
189 self
._inner
_product
_matrix
._cache
= {'hermitian': True}
190 self
._inner
_product
_matrix
.set_immutable()
193 if not self
._is
_jordanian
():
194 raise ValueError("Jordan identity does not hold")
195 if not self
._inner
_product
_is
_associative
():
196 raise ValueError("inner product is not associative")
199 def _coerce_map_from_base_ring(self
):
201 Disable the map from the base ring into the algebra.
203 Performing a nonsense conversion like this automatically
204 is counterpedagogical. The fallback is to try the usual
205 element constructor, which should also fail.
209 sage: from mjo.eja.eja_algebra import random_eja
213 sage: set_random_seed()
214 sage: J = random_eja()
216 Traceback (most recent call last):
218 ValueError: not an element of this algebra
224 def product_on_basis(self
, i
, j
):
225 # We only stored the lower-triangular portion of the
226 # multiplication table.
228 return self
._multiplication
_table
[i
][j
]
230 return self
._multiplication
_table
[j
][i
]
232 def inner_product(self
, x
, y
):
234 The inner product associated with this Euclidean Jordan algebra.
236 Defaults to the trace inner product, but can be overridden by
237 subclasses if they are sure that the necessary properties are
242 sage: from mjo.eja.eja_algebra import (random_eja,
244 ....: BilinearFormEJA)
248 Our inner product is "associative," which means the following for
249 a symmetric bilinear form::
251 sage: set_random_seed()
252 sage: J = random_eja()
253 sage: x,y,z = J.random_elements(3)
254 sage: (x*y).inner_product(z) == y.inner_product(x*z)
259 Ensure that this is the usual inner product for the algebras
262 sage: set_random_seed()
263 sage: J = HadamardEJA.random_instance()
264 sage: x,y = J.random_elements(2)
265 sage: actual = x.inner_product(y)
266 sage: expected = x.to_vector().inner_product(y.to_vector())
267 sage: actual == expected
270 Ensure that this is one-half of the trace inner-product in a
271 BilinearFormEJA that isn't just the reals (when ``n`` isn't
272 one). This is in Faraut and Koranyi, and also my "On the
275 sage: set_random_seed()
276 sage: J = BilinearFormEJA.random_instance()
277 sage: n = J.dimension()
278 sage: x = J.random_element()
279 sage: y = J.random_element()
280 sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
283 B
= self
._inner
_product
_matrix
284 return (B
*x
.to_vector()).inner_product(y
.to_vector())
287 def _is_commutative(self
):
289 Whether or not this algebra's multiplication table is commutative.
291 This method should of course always return ``True``, unless
292 this algebra was constructed with ``check_axioms=False`` and
293 passed an invalid multiplication table.
295 return all( self
.product_on_basis(i
,j
) == self
.product_on_basis(i
,j
)
296 for i
in range(self
.dimension())
297 for j
in range(self
.dimension()) )
299 def _is_jordanian(self
):
301 Whether or not this algebra's multiplication table respects the
302 Jordan identity `(x^{2})(xy) = x(x^{2}y)`.
304 We only check one arrangement of `x` and `y`, so for a
305 ``True`` result to be truly true, you should also check
306 :meth:`_is_commutative`. This method should of course always
307 return ``True``, unless this algebra was constructed with
308 ``check_axioms=False`` and passed an invalid multiplication table.
310 return all( (self
.monomial(i
)**2)*(self
.monomial(i
)*self
.monomial(j
))
312 (self
.monomial(i
))*((self
.monomial(i
)**2)*self
.monomial(j
))
313 for i
in range(self
.dimension())
314 for j
in range(self
.dimension()) )
316 def _inner_product_is_associative(self
):
318 Return whether or not this algebra's inner product `B` is
319 associative; that is, whether or not `B(xy,z) = B(x,yz)`.
321 This method should of course always return ``True``, unless
322 this algebra was constructed with ``check_axioms=False`` and
323 passed an invalid Jordan or inner-product.
326 # Used to check whether or not something is zero in an inexact
327 # ring. This number is sufficient to allow the construction of
328 # QuaternionHermitianEJA(2, field=RDF) with check_axioms=True.
331 for i
in range(self
.dimension()):
332 for j
in range(self
.dimension()):
333 for k
in range(self
.dimension()):
337 diff
= (x
*y
).inner_product(z
) - x
.inner_product(y
*z
)
339 if self
.base_ring().is_exact():
343 if diff
.abs() > epsilon
:
348 def _element_constructor_(self
, elt
):
350 Construct an element of this algebra from its vector or matrix
353 This gets called only after the parent element _call_ method
354 fails to find a coercion for the argument.
358 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
360 ....: RealSymmetricEJA)
364 The identity in `S^n` is converted to the identity in the EJA::
366 sage: J = RealSymmetricEJA(3)
367 sage: I = matrix.identity(QQ,3)
368 sage: J(I) == J.one()
371 This skew-symmetric matrix can't be represented in the EJA::
373 sage: J = RealSymmetricEJA(3)
374 sage: A = matrix(QQ,3, lambda i,j: i-j)
376 Traceback (most recent call last):
378 ValueError: not an element of this algebra
382 Ensure that we can convert any element of the two non-matrix
383 simple algebras (whose matrix representations are columns)
384 back and forth faithfully::
386 sage: set_random_seed()
387 sage: J = HadamardEJA.random_instance()
388 sage: x = J.random_element()
389 sage: J(x.to_vector().column()) == x
391 sage: J = JordanSpinEJA.random_instance()
392 sage: x = J.random_element()
393 sage: J(x.to_vector().column()) == x
397 msg
= "not an element of this algebra"
399 # The superclass implementation of random_element()
400 # needs to be able to coerce "0" into the algebra.
402 elif elt
in self
.base_ring():
403 # Ensure that no base ring -> algebra coercion is performed
404 # by this method. There's some stupidity in sage that would
405 # otherwise propagate to this method; for example, sage thinks
406 # that the integer 3 belongs to the space of 2-by-2 matrices.
407 raise ValueError(msg
)
411 except (AttributeError, TypeError):
412 # Try to convert a vector into a column-matrix
415 if elt
not in self
.matrix_space():
416 raise ValueError(msg
)
418 # Thanks for nothing! Matrix spaces aren't vector spaces in
419 # Sage, so we have to figure out its matrix-basis coordinates
420 # ourselves. We use the basis space's ring instead of the
421 # element's ring because the basis space might be an algebraic
422 # closure whereas the base ring of the 3-by-3 identity matrix
423 # could be QQ instead of QQbar.
425 # We pass check=False because the matrix basis is "guaranteed"
426 # to be linearly independent... right? Ha ha.
427 V
= VectorSpace(self
.base_ring(), elt
.nrows()*elt
.ncols())
428 W
= V
.span_of_basis( (_mat2vec(s
) for s
in self
.matrix_basis()),
432 coords
= W
.coordinate_vector(_mat2vec(elt
))
433 except ArithmeticError: # vector is not in free module
434 raise ValueError(msg
)
436 return self
.from_vector(coords
)
440 Return a string representation of ``self``.
444 sage: from mjo.eja.eja_algebra import JordanSpinEJA
448 Ensure that it says what we think it says::
450 sage: JordanSpinEJA(2, field=AA)
451 Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
452 sage: JordanSpinEJA(3, field=RDF)
453 Euclidean Jordan algebra of dimension 3 over Real Double Field
456 fmt
= "Euclidean Jordan algebra of dimension {} over {}"
457 return fmt
.format(self
.dimension(), self
.base_ring())
461 def characteristic_polynomial_of(self
):
463 Return the algebra's "characteristic polynomial of" function,
464 which is itself a multivariate polynomial that, when evaluated
465 at the coordinates of some algebra element, returns that
466 element's characteristic polynomial.
468 The resulting polynomial has `n+1` variables, where `n` is the
469 dimension of this algebra. The first `n` variables correspond to
470 the coordinates of an algebra element: when evaluated at the
471 coordinates of an algebra element with respect to a certain
472 basis, the result is a univariate polynomial (in the one
473 remaining variable ``t``), namely the characteristic polynomial
478 sage: from mjo.eja.eja_algebra import JordanSpinEJA, TrivialEJA
482 The characteristic polynomial in the spin algebra is given in
483 Alizadeh, Example 11.11::
485 sage: J = JordanSpinEJA(3)
486 sage: p = J.characteristic_polynomial_of(); p
487 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
488 sage: xvec = J.one().to_vector()
492 By definition, the characteristic polynomial is a monic
493 degree-zero polynomial in a rank-zero algebra. Note that
494 Cayley-Hamilton is indeed satisfied since the polynomial
495 ``1`` evaluates to the identity element of the algebra on
498 sage: J = TrivialEJA()
499 sage: J.characteristic_polynomial_of()
506 # The list of coefficient polynomials a_0, a_1, a_2, ..., a_(r-1).
507 a
= self
._charpoly
_coefficients
()
509 # We go to a bit of trouble here to reorder the
510 # indeterminates, so that it's easier to evaluate the
511 # characteristic polynomial at x's coordinates and get back
512 # something in terms of t, which is what we want.
513 S
= PolynomialRing(self
.base_ring(),'t')
517 S
= PolynomialRing(S
, R
.variable_names())
520 return (t
**r
+ sum( a
[k
]*(t
**k
) for k
in range(r
) ))
522 def coordinate_polynomial_ring(self
):
524 The multivariate polynomial ring in which this algebra's
525 :meth:`characteristic_polynomial_of` lives.
529 sage: from mjo.eja.eja_algebra import (HadamardEJA,
530 ....: RealSymmetricEJA)
534 sage: J = HadamardEJA(2)
535 sage: J.coordinate_polynomial_ring()
536 Multivariate Polynomial Ring in X1, X2...
537 sage: J = RealSymmetricEJA(3,field=QQ,orthonormalize=False)
538 sage: J.coordinate_polynomial_ring()
539 Multivariate Polynomial Ring in X1, X2, X3, X4, X5, X6...
542 var_names
= tuple( "X%d" % z
for z
in range(1, self
.dimension()+1) )
543 return PolynomialRing(self
.base_ring(), var_names
)
545 def inner_product(self
, x
, y
):
547 The inner product associated with this Euclidean Jordan algebra.
549 Defaults to the trace inner product, but can be overridden by
550 subclasses if they are sure that the necessary properties are
555 sage: from mjo.eja.eja_algebra import (random_eja,
557 ....: BilinearFormEJA)
561 Our inner product is "associative," which means the following for
562 a symmetric bilinear form::
564 sage: set_random_seed()
565 sage: J = random_eja()
566 sage: x,y,z = J.random_elements(3)
567 sage: (x*y).inner_product(z) == y.inner_product(x*z)
572 Ensure that this is the usual inner product for the algebras
575 sage: set_random_seed()
576 sage: J = HadamardEJA.random_instance()
577 sage: x,y = J.random_elements(2)
578 sage: actual = x.inner_product(y)
579 sage: expected = x.to_vector().inner_product(y.to_vector())
580 sage: actual == expected
583 Ensure that this is one-half of the trace inner-product in a
584 BilinearFormEJA that isn't just the reals (when ``n`` isn't
585 one). This is in Faraut and Koranyi, and also my "On the
588 sage: set_random_seed()
589 sage: J = BilinearFormEJA.random_instance()
590 sage: n = J.dimension()
591 sage: x = J.random_element()
592 sage: y = J.random_element()
593 sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
596 B
= self
._inner
_product
_matrix
597 return (B
*x
.to_vector()).inner_product(y
.to_vector())
600 def is_trivial(self
):
602 Return whether or not this algebra is trivial.
604 A trivial algebra contains only the zero element.
608 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
613 sage: J = ComplexHermitianEJA(3)
619 sage: J = TrivialEJA()
624 return self
.dimension() == 0
627 def multiplication_table(self
):
629 Return a visual representation of this algebra's multiplication
630 table (on basis elements).
634 sage: from mjo.eja.eja_algebra import JordanSpinEJA
638 sage: J = JordanSpinEJA(4)
639 sage: J.multiplication_table()
640 +----++----+----+----+----+
641 | * || e0 | e1 | e2 | e3 |
642 +====++====+====+====+====+
643 | e0 || e0 | e1 | e2 | e3 |
644 +----++----+----+----+----+
645 | e1 || e1 | e0 | 0 | 0 |
646 +----++----+----+----+----+
647 | e2 || e2 | 0 | e0 | 0 |
648 +----++----+----+----+----+
649 | e3 || e3 | 0 | 0 | e0 |
650 +----++----+----+----+----+
654 # Prepend the header row.
655 M
= [["*"] + list(self
.gens())]
657 # And to each subsequent row, prepend an entry that belongs to
658 # the left-side "header column."
659 M
+= [ [self
.monomial(i
)] + [ self
.product_on_basis(i
,j
)
663 return table(M
, header_row
=True, header_column
=True, frame
=True)
666 def matrix_basis(self
):
668 Return an (often more natural) representation of this algebras
669 basis as an ordered tuple of matrices.
671 Every finite-dimensional Euclidean Jordan Algebra is a, up to
672 Jordan isomorphism, a direct sum of five simple
673 algebras---four of which comprise Hermitian matrices. And the
674 last type of algebra can of course be thought of as `n`-by-`1`
675 column matrices (ambiguusly called column vectors) to avoid
676 special cases. As a result, matrices (and column vectors) are
677 a natural representation format for Euclidean Jordan algebra
680 But, when we construct an algebra from a basis of matrices,
681 those matrix representations are lost in favor of coordinate
682 vectors *with respect to* that basis. We could eventually
683 convert back if we tried hard enough, but having the original
684 representations handy is valuable enough that we simply store
685 them and return them from this method.
687 Why implement this for non-matrix algebras? Avoiding special
688 cases for the :class:`BilinearFormEJA` pays with simplicity in
689 its own right. But mainly, we would like to be able to assume
690 that elements of a :class:`DirectSumEJA` can be displayed
691 nicely, without having to have special classes for direct sums
692 one of whose components was a matrix algebra.
696 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
697 ....: RealSymmetricEJA)
701 sage: J = RealSymmetricEJA(2)
703 Finite family {0: e0, 1: e1, 2: e2}
704 sage: J.matrix_basis()
706 [1 0] [ 0 0.7071067811865475?] [0 0]
707 [0 0], [0.7071067811865475? 0], [0 1]
712 sage: J = JordanSpinEJA(2)
714 Finite family {0: e0, 1: e1}
715 sage: J.matrix_basis()
721 return self
._matrix
_basis
724 def matrix_space(self
):
726 Return the matrix space in which this algebra's elements live, if
727 we think of them as matrices (including column vectors of the
730 Generally this will be an `n`-by-`1` column-vector space,
731 except when the algebra is trivial. There it's `n`-by-`n`
732 (where `n` is zero), to ensure that two elements of the matrix
733 space (empty matrices) can be multiplied.
735 Matrix algebras override this with something more useful.
737 if self
.is_trivial():
738 return MatrixSpace(self
.base_ring(), 0)
740 return self
.matrix_basis()[0].parent()
746 Return the unit element of this algebra.
750 sage: from mjo.eja.eja_algebra import (HadamardEJA,
755 We can compute unit element in the Hadamard EJA::
757 sage: J = HadamardEJA(5)
759 e0 + e1 + e2 + e3 + e4
761 The unit element in the Hadamard EJA is inherited in the
762 subalgebras generated by its elements::
764 sage: J = HadamardEJA(5)
766 e0 + e1 + e2 + e3 + e4
767 sage: x = sum(J.gens())
768 sage: A = x.subalgebra_generated_by(orthonormalize=False)
771 sage: A.one().superalgebra_element()
772 e0 + e1 + e2 + e3 + e4
776 The identity element acts like the identity, regardless of
777 whether or not we orthonormalize::
779 sage: set_random_seed()
780 sage: J = random_eja()
781 sage: x = J.random_element()
782 sage: J.one()*x == x and x*J.one() == x
784 sage: A = x.subalgebra_generated_by()
785 sage: y = A.random_element()
786 sage: A.one()*y == y and y*A.one() == y
791 sage: set_random_seed()
792 sage: J = random_eja(field=QQ, orthonormalize=False)
793 sage: x = J.random_element()
794 sage: J.one()*x == x and x*J.one() == x
796 sage: A = x.subalgebra_generated_by(orthonormalize=False)
797 sage: y = A.random_element()
798 sage: A.one()*y == y and y*A.one() == y
801 The matrix of the unit element's operator is the identity,
802 regardless of the base field and whether or not we
805 sage: set_random_seed()
806 sage: J = random_eja()
807 sage: actual = J.one().operator().matrix()
808 sage: expected = matrix.identity(J.base_ring(), J.dimension())
809 sage: actual == expected
811 sage: x = J.random_element()
812 sage: A = x.subalgebra_generated_by()
813 sage: actual = A.one().operator().matrix()
814 sage: expected = matrix.identity(A.base_ring(), A.dimension())
815 sage: actual == expected
820 sage: set_random_seed()
821 sage: J = random_eja(field=QQ, orthonormalize=False)
822 sage: actual = J.one().operator().matrix()
823 sage: expected = matrix.identity(J.base_ring(), J.dimension())
824 sage: actual == expected
826 sage: x = J.random_element()
827 sage: A = x.subalgebra_generated_by(orthonormalize=False)
828 sage: actual = A.one().operator().matrix()
829 sage: expected = matrix.identity(A.base_ring(), A.dimension())
830 sage: actual == expected
833 Ensure that the cached unit element (often precomputed by
834 hand) agrees with the computed one::
836 sage: set_random_seed()
837 sage: J = random_eja()
838 sage: cached = J.one()
839 sage: J.one.clear_cache()
840 sage: J.one() == cached
845 sage: set_random_seed()
846 sage: J = random_eja(field=QQ, orthonormalize=False)
847 sage: cached = J.one()
848 sage: J.one.clear_cache()
849 sage: J.one() == cached
853 # We can brute-force compute the matrices of the operators
854 # that correspond to the basis elements of this algebra.
855 # If some linear combination of those basis elements is the
856 # algebra identity, then the same linear combination of
857 # their matrices has to be the identity matrix.
859 # Of course, matrices aren't vectors in sage, so we have to
860 # appeal to the "long vectors" isometry.
861 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
863 # Now we use basic linear algebra to find the coefficients,
864 # of the matrices-as-vectors-linear-combination, which should
865 # work for the original algebra basis too.
866 A
= matrix(self
.base_ring(), oper_vecs
)
868 # We used the isometry on the left-hand side already, but we
869 # still need to do it for the right-hand side. Recall that we
870 # wanted something that summed to the identity matrix.
871 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
873 # Now if there's an identity element in the algebra, this
874 # should work. We solve on the left to avoid having to
875 # transpose the matrix "A".
876 return self
.from_vector(A
.solve_left(b
))
879 def peirce_decomposition(self
, c
):
881 The Peirce decomposition of this algebra relative to the
884 In the future, this can be extended to a complete system of
885 orthogonal idempotents.
889 - ``c`` -- an idempotent of this algebra.
893 A triple (J0, J5, J1) containing two subalgebras and one subspace
896 - ``J0`` -- the algebra on the eigenspace of ``c.operator()``
897 corresponding to the eigenvalue zero.
899 - ``J5`` -- the eigenspace (NOT a subalgebra) of ``c.operator()``
900 corresponding to the eigenvalue one-half.
902 - ``J1`` -- the algebra on the eigenspace of ``c.operator()``
903 corresponding to the eigenvalue one.
905 These are the only possible eigenspaces for that operator, and this
906 algebra is a direct sum of them. The spaces ``J0`` and ``J1`` are
907 orthogonal, and are subalgebras of this algebra with the appropriate
912 sage: from mjo.eja.eja_algebra import random_eja, RealSymmetricEJA
916 The canonical example comes from the symmetric matrices, which
917 decompose into diagonal and off-diagonal parts::
919 sage: J = RealSymmetricEJA(3)
920 sage: C = matrix(QQ, [ [1,0,0],
924 sage: J0,J5,J1 = J.peirce_decomposition(c)
926 Euclidean Jordan algebra of dimension 1...
928 Vector space of degree 6 and dimension 2...
930 Euclidean Jordan algebra of dimension 3...
931 sage: J0.one().to_matrix()
935 sage: orig_df = AA.options.display_format
936 sage: AA.options.display_format = 'radical'
937 sage: J.from_vector(J5.basis()[0]).to_matrix()
941 sage: J.from_vector(J5.basis()[1]).to_matrix()
945 sage: AA.options.display_format = orig_df
946 sage: J1.one().to_matrix()
953 Every algebra decomposes trivially with respect to its identity
956 sage: set_random_seed()
957 sage: J = random_eja()
958 sage: J0,J5,J1 = J.peirce_decomposition(J.one())
959 sage: J0.dimension() == 0 and J5.dimension() == 0
961 sage: J1.superalgebra() == J and J1.dimension() == J.dimension()
964 The decomposition is into eigenspaces, and its components are
965 therefore necessarily orthogonal. Moreover, the identity
966 elements in the two subalgebras are the projections onto their
967 respective subspaces of the superalgebra's identity element::
969 sage: set_random_seed()
970 sage: J = random_eja()
971 sage: x = J.random_element()
972 sage: if not J.is_trivial():
973 ....: while x.is_nilpotent():
974 ....: x = J.random_element()
975 sage: c = x.subalgebra_idempotent()
976 sage: J0,J5,J1 = J.peirce_decomposition(c)
978 sage: for (w,y,z) in zip(J0.basis(), J5.basis(), J1.basis()):
979 ....: w = w.superalgebra_element()
980 ....: y = J.from_vector(y)
981 ....: z = z.superalgebra_element()
982 ....: ipsum += w.inner_product(y).abs()
983 ....: ipsum += w.inner_product(z).abs()
984 ....: ipsum += y.inner_product(z).abs()
987 sage: J1(c) == J1.one()
989 sage: J0(J.one() - c) == J0.one()
993 if not c
.is_idempotent():
994 raise ValueError("element is not idempotent: %s" % c
)
996 from mjo
.eja
.eja_subalgebra
import FiniteDimensionalEJASubalgebra
998 # Default these to what they should be if they turn out to be
999 # trivial, because eigenspaces_left() won't return eigenvalues
1000 # corresponding to trivial spaces (e.g. it returns only the
1001 # eigenspace corresponding to lambda=1 if you take the
1002 # decomposition relative to the identity element).
1003 trivial
= FiniteDimensionalEJASubalgebra(self
, ())
1004 J0
= trivial
# eigenvalue zero
1005 J5
= VectorSpace(self
.base_ring(), 0) # eigenvalue one-half
1006 J1
= trivial
# eigenvalue one
1008 for (eigval
, eigspace
) in c
.operator().matrix().right_eigenspaces():
1009 if eigval
== ~
(self
.base_ring()(2)):
1012 gens
= tuple( self
.from_vector(b
) for b
in eigspace
.basis() )
1013 subalg
= FiniteDimensionalEJASubalgebra(self
,
1021 raise ValueError("unexpected eigenvalue: %s" % eigval
)
1026 def random_element(self
, thorough
=False):
1028 Return a random element of this algebra.
1030 Our algebra superclass method only returns a linear
1031 combination of at most two basis elements. We instead
1032 want the vector space "random element" method that
1033 returns a more diverse selection.
1037 - ``thorough`` -- (boolean; default False) whether or not we
1038 should generate irrational coefficients for the random
1039 element when our base ring is irrational; this slows the
1040 algebra operations to a crawl, but any truly random method
1044 # For a general base ring... maybe we can trust this to do the
1045 # right thing? Unlikely, but.
1046 V
= self
.vector_space()
1047 v
= V
.random_element()
1049 if self
.base_ring() is AA
:
1050 # The "random element" method of the algebraic reals is
1051 # stupid at the moment, and only returns integers between
1052 # -2 and 2, inclusive:
1054 # https://trac.sagemath.org/ticket/30875
1056 # Instead, we implement our own "random vector" method,
1057 # and then coerce that into the algebra. We use the vector
1058 # space degree here instead of the dimension because a
1059 # subalgebra could (for example) be spanned by only two
1060 # vectors, each with five coordinates. We need to
1061 # generate all five coordinates.
1063 v
*= QQbar
.random_element().real()
1065 v
*= QQ
.random_element()
1067 return self
.from_vector(V
.coordinate_vector(v
))
1069 def random_elements(self
, count
, thorough
=False):
1071 Return ``count`` random elements as a tuple.
1075 - ``thorough`` -- (boolean; default False) whether or not we
1076 should generate irrational coefficients for the random
1077 elements when our base ring is irrational; this slows the
1078 algebra operations to a crawl, but any truly random method
1083 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1087 sage: J = JordanSpinEJA(3)
1088 sage: x,y,z = J.random_elements(3)
1089 sage: all( [ x in J, y in J, z in J ])
1091 sage: len( J.random_elements(10) ) == 10
1095 return tuple( self
.random_element(thorough
)
1096 for idx
in range(count
) )
1100 def _charpoly_coefficients(self
):
1102 The `r` polynomial coefficients of the "characteristic polynomial
1107 sage: from mjo.eja.eja_algebra import random_eja
1111 The theory shows that these are all homogeneous polynomials of
1114 sage: set_random_seed()
1115 sage: J = random_eja()
1116 sage: all(p.is_homogeneous() for p in J._charpoly_coefficients())
1120 n
= self
.dimension()
1121 R
= self
.coordinate_polynomial_ring()
1123 F
= R
.fraction_field()
1126 # From a result in my book, these are the entries of the
1127 # basis representation of L_x.
1128 return sum( vars[k
]*self
.monomial(k
).operator().matrix()[i
,j
]
1131 L_x
= matrix(F
, n
, n
, L_x_i_j
)
1134 if self
.rank
.is_in_cache():
1136 # There's no need to pad the system with redundant
1137 # columns if we *know* they'll be redundant.
1140 # Compute an extra power in case the rank is equal to
1141 # the dimension (otherwise, we would stop at x^(r-1)).
1142 x_powers
= [ (L_x
**k
)*self
.one().to_vector()
1143 for k
in range(n
+1) ]
1144 A
= matrix
.column(F
, x_powers
[:n
])
1145 AE
= A
.extended_echelon_form()
1152 # The theory says that only the first "r" coefficients are
1153 # nonzero, and they actually live in the original polynomial
1154 # ring and not the fraction field. We negate them because
1155 # in the actual characteristic polynomial, they get moved
1156 # to the other side where x^r lives.
1157 return -A_rref
.solve_right(E
*b
).change_ring(R
)[:r
]
1162 Return the rank of this EJA.
1164 This is a cached method because we know the rank a priori for
1165 all of the algebras we can construct. Thus we can avoid the
1166 expensive ``_charpoly_coefficients()`` call unless we truly
1167 need to compute the whole characteristic polynomial.
1171 sage: from mjo.eja.eja_algebra import (HadamardEJA,
1172 ....: JordanSpinEJA,
1173 ....: RealSymmetricEJA,
1174 ....: ComplexHermitianEJA,
1175 ....: QuaternionHermitianEJA,
1180 The rank of the Jordan spin algebra is always two::
1182 sage: JordanSpinEJA(2).rank()
1184 sage: JordanSpinEJA(3).rank()
1186 sage: JordanSpinEJA(4).rank()
1189 The rank of the `n`-by-`n` Hermitian real, complex, or
1190 quaternion matrices is `n`::
1192 sage: RealSymmetricEJA(4).rank()
1194 sage: ComplexHermitianEJA(3).rank()
1196 sage: QuaternionHermitianEJA(2).rank()
1201 Ensure that every EJA that we know how to construct has a
1202 positive integer rank, unless the algebra is trivial in
1203 which case its rank will be zero::
1205 sage: set_random_seed()
1206 sage: J = random_eja()
1210 sage: r > 0 or (r == 0 and J.is_trivial())
1213 Ensure that computing the rank actually works, since the ranks
1214 of all simple algebras are known and will be cached by default::
1216 sage: set_random_seed() # long time
1217 sage: J = random_eja() # long time
1218 sage: caches = J.rank() # long time
1219 sage: J.rank.clear_cache() # long time
1220 sage: J.rank() == cached # long time
1224 return len(self
._charpoly
_coefficients
())
1227 def vector_space(self
):
1229 Return the vector space that underlies this algebra.
1233 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1237 sage: J = RealSymmetricEJA(2)
1238 sage: J.vector_space()
1239 Vector space of dimension 3 over...
1242 return self
.zero().to_vector().parent().ambient_vector_space()
1245 Element
= FiniteDimensionalEJAElement
1247 class RationalBasisEJA(FiniteDimensionalEJA
):
1249 New class for algebras whose supplied basis elements have all rational entries.
1253 sage: from mjo.eja.eja_algebra import BilinearFormEJA
1257 The supplied basis is orthonormalized by default::
1259 sage: B = matrix(QQ, [[1, 0, 0], [0, 25, -32], [0, -32, 41]])
1260 sage: J = BilinearFormEJA(B)
1261 sage: J.matrix_basis()
1278 # Abuse the check_field parameter to check that the entries of
1279 # out basis (in ambient coordinates) are in the field QQ.
1280 if not all( all(b_i
in QQ
for b_i
in b
.list()) for b
in basis
):
1281 raise TypeError("basis not rational")
1283 self
._rational
_algebra
= None
1285 # There's no point in constructing the extra algebra if this
1286 # one is already rational.
1288 # Note: the same Jordan and inner-products work here,
1289 # because they are necessarily defined with respect to
1290 # ambient coordinates and not any particular basis.
1291 self
._rational
_algebra
= FiniteDimensionalEJA(
1296 orthonormalize
=False,
1300 super().__init
__(basis
,
1304 check_field
=check_field
,
1308 def _charpoly_coefficients(self
):
1312 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
1313 ....: JordanSpinEJA)
1317 The base ring of the resulting polynomial coefficients is what
1318 it should be, and not the rationals (unless the algebra was
1319 already over the rationals)::
1321 sage: J = JordanSpinEJA(3)
1322 sage: J._charpoly_coefficients()
1323 (X1^2 - X2^2 - X3^2, -2*X1)
1324 sage: a0 = J._charpoly_coefficients()[0]
1326 Algebraic Real Field
1327 sage: a0.base_ring()
1328 Algebraic Real Field
1331 if self
._rational
_algebra
is None:
1332 # There's no need to construct *another* algebra over the
1333 # rationals if this one is already over the
1334 # rationals. Likewise, if we never orthonormalized our
1335 # basis, we might as well just use the given one.
1336 return super()._charpoly
_coefficients
()
1338 # Do the computation over the rationals. The answer will be
1339 # the same, because all we've done is a change of basis.
1340 # Then, change back from QQ to our real base ring
1341 a
= ( a_i
.change_ring(self
.base_ring())
1342 for a_i
in self
._rational
_algebra
._charpoly
_coefficients
() )
1344 if self
._deortho
_matrix
is None:
1345 # This can happen if our base ring was, say, AA and we
1346 # chose not to (or didn't need to) orthonormalize. It's
1347 # still faster to do the computations over QQ even if
1348 # the numbers in the boxes stay the same.
1351 # Otherwise, convert the coordinate variables back to the
1352 # deorthonormalized ones.
1353 R
= self
.coordinate_polynomial_ring()
1354 from sage
.modules
.free_module_element
import vector
1355 X
= vector(R
, R
.gens())
1356 BX
= self
._deortho
_matrix
*X
1358 subs_dict
= { X[i]: BX[i] for i in range(len(X)) }
1359 return tuple( a_i
.subs(subs_dict
) for a_i
in a
)
1361 class ConcreteEJA(RationalBasisEJA
):
1363 A class for the Euclidean Jordan algebras that we know by name.
1365 These are the Jordan algebras whose basis, multiplication table,
1366 rank, and so on are known a priori. More to the point, they are
1367 the Euclidean Jordan algebras for which we are able to conjure up
1368 a "random instance."
1372 sage: from mjo.eja.eja_algebra import ConcreteEJA
1376 Our basis is normalized with respect to the algebra's inner
1377 product, unless we specify otherwise::
1379 sage: set_random_seed()
1380 sage: J = ConcreteEJA.random_instance()
1381 sage: all( b.norm() == 1 for b in J.gens() )
1384 Since our basis is orthonormal with respect to the algebra's inner
1385 product, and since we know that this algebra is an EJA, any
1386 left-multiplication operator's matrix will be symmetric because
1387 natural->EJA basis representation is an isometry and within the
1388 EJA the operator is self-adjoint by the Jordan axiom::
1390 sage: set_random_seed()
1391 sage: J = ConcreteEJA.random_instance()
1392 sage: x = J.random_element()
1393 sage: x.operator().is_self_adjoint()
1398 def _max_random_instance_size():
1400 Return an integer "size" that is an upper bound on the size of
1401 this algebra when it is used in a random test
1402 case. Unfortunately, the term "size" is ambiguous -- when
1403 dealing with `R^n` under either the Hadamard or Jordan spin
1404 product, the "size" refers to the dimension `n`. When dealing
1405 with a matrix algebra (real symmetric or complex/quaternion
1406 Hermitian), it refers to the size of the matrix, which is far
1407 less than the dimension of the underlying vector space.
1409 This method must be implemented in each subclass.
1411 raise NotImplementedError
1414 def random_instance(cls
, *args
, **kwargs
):
1416 Return a random instance of this type of algebra.
1418 This method should be implemented in each subclass.
1420 from sage
.misc
.prandom
import choice
1421 eja_class
= choice(cls
.__subclasses
__())
1423 # These all bubble up to the RationalBasisEJA superclass
1424 # constructor, so any (kw)args valid there are also valid
1426 return eja_class
.random_instance(*args
, **kwargs
)
1431 def dimension_over_reals():
1433 The dimension of this matrix's base ring over the reals.
1435 The reals are dimension one over themselves, obviously; that's
1436 just `\mathbb{R}^{1}`. Likewise, the complex numbers `a + bi`
1437 have dimension two. Finally, the quaternions have dimension
1438 four over the reals.
1440 This is used to determine the size of the matrix returned from
1441 :meth:`real_embed`, among other things.
1443 raise NotImplementedError
1446 def real_embed(cls
,M
):
1448 Embed the matrix ``M`` into a space of real matrices.
1450 The matrix ``M`` can have entries in any field at the moment:
1451 the real numbers, complex numbers, or quaternions. And although
1452 they are not a field, we can probably support octonions at some
1453 point, too. This function returns a real matrix that "acts like"
1454 the original with respect to matrix multiplication; i.e.
1456 real_embed(M*N) = real_embed(M)*real_embed(N)
1459 if M
.ncols() != M
.nrows():
1460 raise ValueError("the matrix 'M' must be square")
1465 def real_unembed(cls
,M
):
1467 The inverse of :meth:`real_embed`.
1469 if M
.ncols() != M
.nrows():
1470 raise ValueError("the matrix 'M' must be square")
1471 if not ZZ(M
.nrows()).mod(cls
.dimension_over_reals()).is_zero():
1472 raise ValueError("the matrix 'M' must be a real embedding")
1476 def jordan_product(X
,Y
):
1477 return (X
*Y
+ Y
*X
)/2
1480 def trace_inner_product(cls
,X
,Y
):
1482 Compute the trace inner-product of two real-embeddings.
1486 sage: from mjo.eja.eja_algebra import (RealSymmetricEJA,
1487 ....: ComplexHermitianEJA,
1488 ....: QuaternionHermitianEJA)
1492 This gives the same answer as it would if we computed the trace
1493 from the unembedded (original) matrices::
1495 sage: set_random_seed()
1496 sage: J = RealSymmetricEJA.random_instance()
1497 sage: x,y = J.random_elements(2)
1498 sage: Xe = x.to_matrix()
1499 sage: Ye = y.to_matrix()
1500 sage: X = J.real_unembed(Xe)
1501 sage: Y = J.real_unembed(Ye)
1502 sage: expected = (X*Y).trace()
1503 sage: actual = J.trace_inner_product(Xe,Ye)
1504 sage: actual == expected
1509 sage: set_random_seed()
1510 sage: J = ComplexHermitianEJA.random_instance()
1511 sage: x,y = J.random_elements(2)
1512 sage: Xe = x.to_matrix()
1513 sage: Ye = y.to_matrix()
1514 sage: X = J.real_unembed(Xe)
1515 sage: Y = J.real_unembed(Ye)
1516 sage: expected = (X*Y).trace().real()
1517 sage: actual = J.trace_inner_product(Xe,Ye)
1518 sage: actual == expected
1523 sage: set_random_seed()
1524 sage: J = QuaternionHermitianEJA.random_instance()
1525 sage: x,y = J.random_elements(2)
1526 sage: Xe = x.to_matrix()
1527 sage: Ye = y.to_matrix()
1528 sage: X = J.real_unembed(Xe)
1529 sage: Y = J.real_unembed(Ye)
1530 sage: expected = (X*Y).trace().coefficient_tuple()[0]
1531 sage: actual = J.trace_inner_product(Xe,Ye)
1532 sage: actual == expected
1536 Xu
= cls
.real_unembed(X
)
1537 Yu
= cls
.real_unembed(Y
)
1538 tr
= (Xu
*Yu
).trace()
1541 # Works in QQ, AA, RDF, et cetera.
1543 except AttributeError:
1544 # A quaternion doesn't have a real() method, but does
1545 # have coefficient_tuple() method that returns the
1546 # coefficients of 1, i, j, and k -- in that order.
1547 return tr
.coefficient_tuple()[0]
1550 class RealMatrixEJA(MatrixEJA
):
1552 def dimension_over_reals():
1556 class RealSymmetricEJA(ConcreteEJA
, RealMatrixEJA
):
1558 The rank-n simple EJA consisting of real symmetric n-by-n
1559 matrices, the usual symmetric Jordan product, and the trace inner
1560 product. It has dimension `(n^2 + n)/2` over the reals.
1564 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1568 sage: J = RealSymmetricEJA(2)
1569 sage: e0, e1, e2 = J.gens()
1577 In theory, our "field" can be any subfield of the reals::
1579 sage: RealSymmetricEJA(2, field=RDF)
1580 Euclidean Jordan algebra of dimension 3 over Real Double Field
1581 sage: RealSymmetricEJA(2, field=RR)
1582 Euclidean Jordan algebra of dimension 3 over Real Field with
1583 53 bits of precision
1587 The dimension of this algebra is `(n^2 + n) / 2`::
1589 sage: set_random_seed()
1590 sage: n_max = RealSymmetricEJA._max_random_instance_size()
1591 sage: n = ZZ.random_element(1, n_max)
1592 sage: J = RealSymmetricEJA(n)
1593 sage: J.dimension() == (n^2 + n)/2
1596 The Jordan multiplication is what we think it is::
1598 sage: set_random_seed()
1599 sage: J = RealSymmetricEJA.random_instance()
1600 sage: x,y = J.random_elements(2)
1601 sage: actual = (x*y).to_matrix()
1602 sage: X = x.to_matrix()
1603 sage: Y = y.to_matrix()
1604 sage: expected = (X*Y + Y*X)/2
1605 sage: actual == expected
1607 sage: J(expected) == x*y
1610 We can change the generator prefix::
1612 sage: RealSymmetricEJA(3, prefix='q').gens()
1613 (q0, q1, q2, q3, q4, q5)
1615 We can construct the (trivial) algebra of rank zero::
1617 sage: RealSymmetricEJA(0)
1618 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1622 def _denormalized_basis(cls
, n
):
1624 Return a basis for the space of real symmetric n-by-n matrices.
1628 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1632 sage: set_random_seed()
1633 sage: n = ZZ.random_element(1,5)
1634 sage: B = RealSymmetricEJA._denormalized_basis(n)
1635 sage: all( M.is_symmetric() for M in B)
1639 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1643 for j
in range(i
+1):
1644 Eij
= matrix(ZZ
, n
, lambda k
,l
: k
==i
and l
==j
)
1648 Sij
= Eij
+ Eij
.transpose()
1654 def _max_random_instance_size():
1655 return 4 # Dimension 10
1658 def random_instance(cls
, **kwargs
):
1660 Return a random instance of this type of algebra.
1662 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1663 return cls(n
, **kwargs
)
1665 def __init__(self
, n
, **kwargs
):
1666 # We know this is a valid EJA, but will double-check
1667 # if the user passes check_axioms=True.
1668 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
1670 super(RealSymmetricEJA
, self
).__init
__(self
._denormalized
_basis
(n
),
1671 self
.jordan_product
,
1672 self
.trace_inner_product
,
1675 # TODO: this could be factored out somehow, but is left here
1676 # because the MatrixEJA is not presently a subclass of the
1677 # FDEJA class that defines rank() and one().
1678 self
.rank
.set_cache(n
)
1679 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
1680 self
.one
.set_cache(self(idV
))
1684 class ComplexMatrixEJA(MatrixEJA
):
1685 # A manual dictionary-cache for the complex_extension() method,
1686 # since apparently @classmethods can't also be @cached_methods.
1687 _complex_extension
= {}
1690 def complex_extension(cls
,field
):
1692 The complex field that we embed/unembed, as an extension
1693 of the given ``field``.
1695 if field
in cls
._complex
_extension
:
1696 return cls
._complex
_extension
[field
]
1698 # Sage doesn't know how to adjoin the complex "i" (the root of
1699 # x^2 + 1) to a field in a general way. Here, we just enumerate
1700 # all of the cases that I have cared to support so far.
1702 # Sage doesn't know how to embed AA into QQbar, i.e. how
1703 # to adjoin sqrt(-1) to AA.
1705 elif not field
.is_exact():
1707 F
= field
.complex_field()
1709 # Works for QQ and... maybe some other fields.
1710 R
= PolynomialRing(field
, 'z')
1712 F
= field
.extension(z
**2 + 1, 'I', embedding
=CLF(-1).sqrt())
1714 cls
._complex
_extension
[field
] = F
1718 def dimension_over_reals():
1722 def real_embed(cls
,M
):
1724 Embed the n-by-n complex matrix ``M`` into the space of real
1725 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1726 bi` to the block matrix ``[[a,b],[-b,a]]``.
1730 sage: from mjo.eja.eja_algebra import ComplexMatrixEJA
1734 sage: F = QuadraticField(-1, 'I')
1735 sage: x1 = F(4 - 2*i)
1736 sage: x2 = F(1 + 2*i)
1739 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1740 sage: ComplexMatrixEJA.real_embed(M)
1749 Embedding is a homomorphism (isomorphism, in fact)::
1751 sage: set_random_seed()
1752 sage: n = ZZ.random_element(3)
1753 sage: F = QuadraticField(-1, 'I')
1754 sage: X = random_matrix(F, n)
1755 sage: Y = random_matrix(F, n)
1756 sage: Xe = ComplexMatrixEJA.real_embed(X)
1757 sage: Ye = ComplexMatrixEJA.real_embed(Y)
1758 sage: XYe = ComplexMatrixEJA.real_embed(X*Y)
1763 super(ComplexMatrixEJA
,cls
).real_embed(M
)
1766 # We don't need any adjoined elements...
1767 field
= M
.base_ring().base_ring()
1773 blocks
.append(matrix(field
, 2, [ [ a
, b
],
1776 return matrix
.block(field
, n
, blocks
)
1780 def real_unembed(cls
,M
):
1782 The inverse of _embed_complex_matrix().
1786 sage: from mjo.eja.eja_algebra import ComplexMatrixEJA
1790 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1791 ....: [-2, 1, -4, 3],
1792 ....: [ 9, 10, 11, 12],
1793 ....: [-10, 9, -12, 11] ])
1794 sage: ComplexMatrixEJA.real_unembed(A)
1796 [ 10*I + 9 12*I + 11]
1800 Unembedding is the inverse of embedding::
1802 sage: set_random_seed()
1803 sage: F = QuadraticField(-1, 'I')
1804 sage: M = random_matrix(F, 3)
1805 sage: Me = ComplexMatrixEJA.real_embed(M)
1806 sage: ComplexMatrixEJA.real_unembed(Me) == M
1810 super(ComplexMatrixEJA
,cls
).real_unembed(M
)
1812 d
= cls
.dimension_over_reals()
1813 F
= cls
.complex_extension(M
.base_ring())
1816 # Go top-left to bottom-right (reading order), converting every
1817 # 2-by-2 block we see to a single complex element.
1819 for k
in range(n
/d
):
1820 for j
in range(n
/d
):
1821 submat
= M
[d
*k
:d
*k
+d
,d
*j
:d
*j
+d
]
1822 if submat
[0,0] != submat
[1,1]:
1823 raise ValueError('bad on-diagonal submatrix')
1824 if submat
[0,1] != -submat
[1,0]:
1825 raise ValueError('bad off-diagonal submatrix')
1826 z
= submat
[0,0] + submat
[0,1]*i
1829 return matrix(F
, n
/d
, elements
)
1832 class ComplexHermitianEJA(ConcreteEJA
, ComplexMatrixEJA
):
1834 The rank-n simple EJA consisting of complex Hermitian n-by-n
1835 matrices over the real numbers, the usual symmetric Jordan product,
1836 and the real-part-of-trace inner product. It has dimension `n^2` over
1841 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1845 In theory, our "field" can be any subfield of the reals::
1847 sage: ComplexHermitianEJA(2, field=RDF)
1848 Euclidean Jordan algebra of dimension 4 over Real Double Field
1849 sage: ComplexHermitianEJA(2, field=RR)
1850 Euclidean Jordan algebra of dimension 4 over Real Field with
1851 53 bits of precision
1855 The dimension of this algebra is `n^2`::
1857 sage: set_random_seed()
1858 sage: n_max = ComplexHermitianEJA._max_random_instance_size()
1859 sage: n = ZZ.random_element(1, n_max)
1860 sage: J = ComplexHermitianEJA(n)
1861 sage: J.dimension() == n^2
1864 The Jordan multiplication is what we think it is::
1866 sage: set_random_seed()
1867 sage: J = ComplexHermitianEJA.random_instance()
1868 sage: x,y = J.random_elements(2)
1869 sage: actual = (x*y).to_matrix()
1870 sage: X = x.to_matrix()
1871 sage: Y = y.to_matrix()
1872 sage: expected = (X*Y + Y*X)/2
1873 sage: actual == expected
1875 sage: J(expected) == x*y
1878 We can change the generator prefix::
1880 sage: ComplexHermitianEJA(2, prefix='z').gens()
1883 We can construct the (trivial) algebra of rank zero::
1885 sage: ComplexHermitianEJA(0)
1886 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1891 def _denormalized_basis(cls
, n
):
1893 Returns a basis for the space of complex Hermitian n-by-n matrices.
1895 Why do we embed these? Basically, because all of numerical linear
1896 algebra assumes that you're working with vectors consisting of `n`
1897 entries from a field and scalars from the same field. There's no way
1898 to tell SageMath that (for example) the vectors contain complex
1899 numbers, while the scalar field is real.
1903 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1907 sage: set_random_seed()
1908 sage: n = ZZ.random_element(1,5)
1909 sage: B = ComplexHermitianEJA._denormalized_basis(n)
1910 sage: all( M.is_symmetric() for M in B)
1915 R
= PolynomialRing(field
, 'z')
1917 F
= field
.extension(z
**2 + 1, 'I')
1920 # This is like the symmetric case, but we need to be careful:
1922 # * We want conjugate-symmetry, not just symmetry.
1923 # * The diagonal will (as a result) be real.
1926 Eij
= matrix
.zero(F
,n
)
1928 for j
in range(i
+1):
1932 Sij
= cls
.real_embed(Eij
)
1935 # The second one has a minus because it's conjugated.
1936 Eij
[j
,i
] = 1 # Eij = Eij + Eij.transpose()
1937 Sij_real
= cls
.real_embed(Eij
)
1939 # Eij = I*Eij - I*Eij.transpose()
1942 Sij_imag
= cls
.real_embed(Eij
)
1948 # Since we embedded these, we can drop back to the "field" that we
1949 # started with instead of the complex extension "F".
1950 return tuple( s
.change_ring(field
) for s
in S
)
1953 def __init__(self
, n
, **kwargs
):
1954 # We know this is a valid EJA, but will double-check
1955 # if the user passes check_axioms=True.
1956 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
1958 super(ComplexHermitianEJA
, self
).__init
__(self
._denormalized
_basis
(n
),
1959 self
.jordan_product
,
1960 self
.trace_inner_product
,
1962 # TODO: this could be factored out somehow, but is left here
1963 # because the MatrixEJA is not presently a subclass of the
1964 # FDEJA class that defines rank() and one().
1965 self
.rank
.set_cache(n
)
1966 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
1967 self
.one
.set_cache(self(idV
))
1970 def _max_random_instance_size():
1971 return 3 # Dimension 9
1974 def random_instance(cls
, **kwargs
):
1976 Return a random instance of this type of algebra.
1978 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1979 return cls(n
, **kwargs
)
1981 class QuaternionMatrixEJA(MatrixEJA
):
1983 # A manual dictionary-cache for the quaternion_extension() method,
1984 # since apparently @classmethods can't also be @cached_methods.
1985 _quaternion_extension
= {}
1988 def quaternion_extension(cls
,field
):
1990 The quaternion field that we embed/unembed, as an extension
1991 of the given ``field``.
1993 if field
in cls
._quaternion
_extension
:
1994 return cls
._quaternion
_extension
[field
]
1996 Q
= QuaternionAlgebra(field
,-1,-1)
1998 cls
._quaternion
_extension
[field
] = Q
2002 def dimension_over_reals():
2006 def real_embed(cls
,M
):
2008 Embed the n-by-n quaternion matrix ``M`` into the space of real
2009 matrices of size 4n-by-4n by first sending each quaternion entry `z
2010 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
2011 c+di],[-c + di, a-bi]]`, and then embedding those into a real
2016 sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA
2020 sage: Q = QuaternionAlgebra(QQ,-1,-1)
2021 sage: i,j,k = Q.gens()
2022 sage: x = 1 + 2*i + 3*j + 4*k
2023 sage: M = matrix(Q, 1, [[x]])
2024 sage: QuaternionMatrixEJA.real_embed(M)
2030 Embedding is a homomorphism (isomorphism, in fact)::
2032 sage: set_random_seed()
2033 sage: n = ZZ.random_element(2)
2034 sage: Q = QuaternionAlgebra(QQ,-1,-1)
2035 sage: X = random_matrix(Q, n)
2036 sage: Y = random_matrix(Q, n)
2037 sage: Xe = QuaternionMatrixEJA.real_embed(X)
2038 sage: Ye = QuaternionMatrixEJA.real_embed(Y)
2039 sage: XYe = QuaternionMatrixEJA.real_embed(X*Y)
2044 super(QuaternionMatrixEJA
,cls
).real_embed(M
)
2045 quaternions
= M
.base_ring()
2048 F
= QuadraticField(-1, 'I')
2053 t
= z
.coefficient_tuple()
2058 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
2059 [-c
+ d
*i
, a
- b
*i
]])
2060 realM
= ComplexMatrixEJA
.real_embed(cplxM
)
2061 blocks
.append(realM
)
2063 # We should have real entries by now, so use the realest field
2064 # we've got for the return value.
2065 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
2070 def real_unembed(cls
,M
):
2072 The inverse of _embed_quaternion_matrix().
2076 sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA
2080 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
2081 ....: [-2, 1, -4, 3],
2082 ....: [-3, 4, 1, -2],
2083 ....: [-4, -3, 2, 1]])
2084 sage: QuaternionMatrixEJA.real_unembed(M)
2085 [1 + 2*i + 3*j + 4*k]
2089 Unembedding is the inverse of embedding::
2091 sage: set_random_seed()
2092 sage: Q = QuaternionAlgebra(QQ, -1, -1)
2093 sage: M = random_matrix(Q, 3)
2094 sage: Me = QuaternionMatrixEJA.real_embed(M)
2095 sage: QuaternionMatrixEJA.real_unembed(Me) == M
2099 super(QuaternionMatrixEJA
,cls
).real_unembed(M
)
2101 d
= cls
.dimension_over_reals()
2103 # Use the base ring of the matrix to ensure that its entries can be
2104 # multiplied by elements of the quaternion algebra.
2105 Q
= cls
.quaternion_extension(M
.base_ring())
2108 # Go top-left to bottom-right (reading order), converting every
2109 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
2112 for l
in range(n
/d
):
2113 for m
in range(n
/d
):
2114 submat
= ComplexMatrixEJA
.real_unembed(
2115 M
[d
*l
:d
*l
+d
,d
*m
:d
*m
+d
] )
2116 if submat
[0,0] != submat
[1,1].conjugate():
2117 raise ValueError('bad on-diagonal submatrix')
2118 if submat
[0,1] != -submat
[1,0].conjugate():
2119 raise ValueError('bad off-diagonal submatrix')
2120 z
= submat
[0,0].real()
2121 z
+= submat
[0,0].imag()*i
2122 z
+= submat
[0,1].real()*j
2123 z
+= submat
[0,1].imag()*k
2126 return matrix(Q
, n
/d
, elements
)
2129 class QuaternionHermitianEJA(ConcreteEJA
, QuaternionMatrixEJA
):
2131 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
2132 matrices, the usual symmetric Jordan product, and the
2133 real-part-of-trace inner product. It has dimension `2n^2 - n` over
2138 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2142 In theory, our "field" can be any subfield of the reals::
2144 sage: QuaternionHermitianEJA(2, field=RDF)
2145 Euclidean Jordan algebra of dimension 6 over Real Double Field
2146 sage: QuaternionHermitianEJA(2, field=RR)
2147 Euclidean Jordan algebra of dimension 6 over Real Field with
2148 53 bits of precision
2152 The dimension of this algebra is `2*n^2 - n`::
2154 sage: set_random_seed()
2155 sage: n_max = QuaternionHermitianEJA._max_random_instance_size()
2156 sage: n = ZZ.random_element(1, n_max)
2157 sage: J = QuaternionHermitianEJA(n)
2158 sage: J.dimension() == 2*(n^2) - n
2161 The Jordan multiplication is what we think it is::
2163 sage: set_random_seed()
2164 sage: J = QuaternionHermitianEJA.random_instance()
2165 sage: x,y = J.random_elements(2)
2166 sage: actual = (x*y).to_matrix()
2167 sage: X = x.to_matrix()
2168 sage: Y = y.to_matrix()
2169 sage: expected = (X*Y + Y*X)/2
2170 sage: actual == expected
2172 sage: J(expected) == x*y
2175 We can change the generator prefix::
2177 sage: QuaternionHermitianEJA(2, prefix='a').gens()
2178 (a0, a1, a2, a3, a4, a5)
2180 We can construct the (trivial) algebra of rank zero::
2182 sage: QuaternionHermitianEJA(0)
2183 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2187 def _denormalized_basis(cls
, n
):
2189 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
2191 Why do we embed these? Basically, because all of numerical
2192 linear algebra assumes that you're working with vectors consisting
2193 of `n` entries from a field and scalars from the same field. There's
2194 no way to tell SageMath that (for example) the vectors contain
2195 complex numbers, while the scalar field is real.
2199 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2203 sage: set_random_seed()
2204 sage: n = ZZ.random_element(1,5)
2205 sage: B = QuaternionHermitianEJA._denormalized_basis(n)
2206 sage: all( M.is_symmetric() for M in B )
2211 Q
= QuaternionAlgebra(QQ
,-1,-1)
2214 # This is like the symmetric case, but we need to be careful:
2216 # * We want conjugate-symmetry, not just symmetry.
2217 # * The diagonal will (as a result) be real.
2220 Eij
= matrix
.zero(Q
,n
)
2222 for j
in range(i
+1):
2226 Sij
= cls
.real_embed(Eij
)
2229 # The second, third, and fourth ones have a minus
2230 # because they're conjugated.
2231 # Eij = Eij + Eij.transpose()
2233 Sij_real
= cls
.real_embed(Eij
)
2235 # Eij = I*(Eij - Eij.transpose())
2238 Sij_I
= cls
.real_embed(Eij
)
2240 # Eij = J*(Eij - Eij.transpose())
2243 Sij_J
= cls
.real_embed(Eij
)
2245 # Eij = K*(Eij - Eij.transpose())
2248 Sij_K
= cls
.real_embed(Eij
)
2254 # Since we embedded these, we can drop back to the "field" that we
2255 # started with instead of the quaternion algebra "Q".
2256 return tuple( s
.change_ring(field
) for s
in S
)
2259 def __init__(self
, n
, **kwargs
):
2260 # We know this is a valid EJA, but will double-check
2261 # if the user passes check_axioms=True.
2262 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2264 super(QuaternionHermitianEJA
, self
).__init
__(self
._denormalized
_basis
(n
),
2265 self
.jordan_product
,
2266 self
.trace_inner_product
,
2268 # TODO: this could be factored out somehow, but is left here
2269 # because the MatrixEJA is not presently a subclass of the
2270 # FDEJA class that defines rank() and one().
2271 self
.rank
.set_cache(n
)
2272 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
2273 self
.one
.set_cache(self(idV
))
2277 def _max_random_instance_size():
2279 The maximum rank of a random QuaternionHermitianEJA.
2281 return 2 # Dimension 6
2284 def random_instance(cls
, **kwargs
):
2286 Return a random instance of this type of algebra.
2288 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2289 return cls(n
, **kwargs
)
2292 class HadamardEJA(ConcreteEJA
):
2294 Return the Euclidean Jordan Algebra corresponding to the set
2295 `R^n` under the Hadamard product.
2297 Note: this is nothing more than the Cartesian product of ``n``
2298 copies of the spin algebra. Once Cartesian product algebras
2299 are implemented, this can go.
2303 sage: from mjo.eja.eja_algebra import HadamardEJA
2307 This multiplication table can be verified by hand::
2309 sage: J = HadamardEJA(3)
2310 sage: e0,e1,e2 = J.gens()
2326 We can change the generator prefix::
2328 sage: HadamardEJA(3, prefix='r').gens()
2332 def __init__(self
, n
, **kwargs
):
2334 jordan_product
= lambda x
,y
: x
2335 inner_product
= lambda x
,y
: x
2337 def jordan_product(x
,y
):
2339 return P( xi
*yi
for (xi
,yi
) in zip(x
,y
) )
2341 def inner_product(x
,y
):
2344 # New defaults for keyword arguments. Don't orthonormalize
2345 # because our basis is already orthonormal with respect to our
2346 # inner-product. Don't check the axioms, because we know this
2347 # is a valid EJA... but do double-check if the user passes
2348 # check_axioms=True. Note: we DON'T override the "check_field"
2349 # default here, because the user can pass in a field!
2350 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2351 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2353 column_basis
= tuple( b
.column() for b
in FreeModule(ZZ
, n
).basis() )
2354 super().__init
__(column_basis
, jordan_product
, inner_product
, **kwargs
)
2355 self
.rank
.set_cache(n
)
2358 self
.one
.set_cache( self
.zero() )
2360 self
.one
.set_cache( sum(self
.gens()) )
2363 def _max_random_instance_size():
2365 The maximum dimension of a random HadamardEJA.
2370 def random_instance(cls
, **kwargs
):
2372 Return a random instance of this type of algebra.
2374 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2375 return cls(n
, **kwargs
)
2378 class BilinearFormEJA(ConcreteEJA
):
2380 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2381 with the half-trace inner product and jordan product ``x*y =
2382 (<Bx,y>,y_bar>, x0*y_bar + y0*x_bar)`` where `B = 1 \times B22` is
2383 a symmetric positive-definite "bilinear form" matrix. Its
2384 dimension is the size of `B`, and it has rank two in dimensions
2385 larger than two. It reduces to the ``JordanSpinEJA`` when `B` is
2386 the identity matrix of order ``n``.
2388 We insist that the one-by-one upper-left identity block of `B` be
2389 passed in as well so that we can be passed a matrix of size zero
2390 to construct a trivial algebra.
2394 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
2395 ....: JordanSpinEJA)
2399 When no bilinear form is specified, the identity matrix is used,
2400 and the resulting algebra is the Jordan spin algebra::
2402 sage: B = matrix.identity(AA,3)
2403 sage: J0 = BilinearFormEJA(B)
2404 sage: J1 = JordanSpinEJA(3)
2405 sage: J0.multiplication_table() == J0.multiplication_table()
2408 An error is raised if the matrix `B` does not correspond to a
2409 positive-definite bilinear form::
2411 sage: B = matrix.random(QQ,2,3)
2412 sage: J = BilinearFormEJA(B)
2413 Traceback (most recent call last):
2415 ValueError: bilinear form is not positive-definite
2416 sage: B = matrix.zero(QQ,3)
2417 sage: J = BilinearFormEJA(B)
2418 Traceback (most recent call last):
2420 ValueError: bilinear form is not positive-definite
2424 We can create a zero-dimensional algebra::
2426 sage: B = matrix.identity(AA,0)
2427 sage: J = BilinearFormEJA(B)
2431 We can check the multiplication condition given in the Jordan, von
2432 Neumann, and Wigner paper (and also discussed on my "On the
2433 symmetry..." paper). Note that this relies heavily on the standard
2434 choice of basis, as does anything utilizing the bilinear form
2435 matrix. We opt not to orthonormalize the basis, because if we
2436 did, we would have to normalize the `s_{i}` in a similar manner::
2438 sage: set_random_seed()
2439 sage: n = ZZ.random_element(5)
2440 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2441 sage: B11 = matrix.identity(QQ,1)
2442 sage: B22 = M.transpose()*M
2443 sage: B = block_matrix(2,2,[ [B11,0 ],
2445 sage: J = BilinearFormEJA(B, orthonormalize=False)
2446 sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
2447 sage: V = J.vector_space()
2448 sage: sis = [ J( V([0] + (M.inverse()*ei).list()).column() )
2449 ....: for ei in eis ]
2450 sage: actual = [ sis[i]*sis[j]
2451 ....: for i in range(n-1)
2452 ....: for j in range(n-1) ]
2453 sage: expected = [ J.one() if i == j else J.zero()
2454 ....: for i in range(n-1)
2455 ....: for j in range(n-1) ]
2456 sage: actual == expected
2460 def __init__(self
, B
, **kwargs
):
2461 # The matrix "B" is supplied by the user in most cases,
2462 # so it makes sense to check whether or not its positive-
2463 # definite unless we are specifically asked not to...
2464 if ("check_axioms" not in kwargs
) or kwargs
["check_axioms"]:
2465 if not B
.is_positive_definite():
2466 raise ValueError("bilinear form is not positive-definite")
2468 # However, all of the other data for this EJA is computed
2469 # by us in manner that guarantees the axioms are
2470 # satisfied. So, again, unless we are specifically asked to
2471 # verify things, we'll skip the rest of the checks.
2472 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2474 def inner_product(x
,y
):
2475 return (y
.T
*B
*x
)[0,0]
2477 def jordan_product(x
,y
):
2483 z0
= inner_product(y
,x
)
2484 zbar
= y0
*xbar
+ x0
*ybar
2485 return P([z0
] + zbar
.list())
2488 column_basis
= tuple( b
.column() for b
in FreeModule(ZZ
, n
).basis() )
2489 super(BilinearFormEJA
, self
).__init
__(column_basis
,
2494 # The rank of this algebra is two, unless we're in a
2495 # one-dimensional ambient space (because the rank is bounded
2496 # by the ambient dimension).
2497 self
.rank
.set_cache(min(n
,2))
2500 self
.one
.set_cache( self
.zero() )
2502 self
.one
.set_cache( self
.monomial(0) )
2505 def _max_random_instance_size():
2507 The maximum dimension of a random BilinearFormEJA.
2512 def random_instance(cls
, **kwargs
):
2514 Return a random instance of this algebra.
2516 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2518 B
= matrix
.identity(ZZ
, n
)
2519 return cls(B
, **kwargs
)
2521 B11
= matrix
.identity(ZZ
, 1)
2522 M
= matrix
.random(ZZ
, n
-1)
2523 I
= matrix
.identity(ZZ
, n
-1)
2525 while alpha
.is_zero():
2526 alpha
= ZZ
.random_element().abs()
2527 B22
= M
.transpose()*M
+ alpha
*I
2529 from sage
.matrix
.special
import block_matrix
2530 B
= block_matrix(2,2, [ [B11
, ZZ(0) ],
2533 return cls(B
, **kwargs
)
2536 class JordanSpinEJA(BilinearFormEJA
):
2538 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2539 with the usual inner product and jordan product ``x*y =
2540 (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2545 sage: from mjo.eja.eja_algebra import JordanSpinEJA
2549 This multiplication table can be verified by hand::
2551 sage: J = JordanSpinEJA(4)
2552 sage: e0,e1,e2,e3 = J.gens()
2568 We can change the generator prefix::
2570 sage: JordanSpinEJA(2, prefix='B').gens()
2575 Ensure that we have the usual inner product on `R^n`::
2577 sage: set_random_seed()
2578 sage: J = JordanSpinEJA.random_instance()
2579 sage: x,y = J.random_elements(2)
2580 sage: actual = x.inner_product(y)
2581 sage: expected = x.to_vector().inner_product(y.to_vector())
2582 sage: actual == expected
2586 def __init__(self
, n
, **kwargs
):
2587 # This is a special case of the BilinearFormEJA with the
2588 # identity matrix as its bilinear form.
2589 B
= matrix
.identity(ZZ
, n
)
2591 # Don't orthonormalize because our basis is already
2592 # orthonormal with respect to our inner-product.
2593 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2595 # But also don't pass check_field=False here, because the user
2596 # can pass in a field!
2597 super(JordanSpinEJA
, self
).__init
__(B
, **kwargs
)
2600 def _max_random_instance_size():
2602 The maximum dimension of a random JordanSpinEJA.
2607 def random_instance(cls
, **kwargs
):
2609 Return a random instance of this type of algebra.
2611 Needed here to override the implementation for ``BilinearFormEJA``.
2613 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2614 return cls(n
, **kwargs
)
2617 class TrivialEJA(ConcreteEJA
):
2619 The trivial Euclidean Jordan algebra consisting of only a zero element.
2623 sage: from mjo.eja.eja_algebra import TrivialEJA
2627 sage: J = TrivialEJA()
2634 sage: 7*J.one()*12*J.one()
2636 sage: J.one().inner_product(J.one())
2638 sage: J.one().norm()
2640 sage: J.one().subalgebra_generated_by()
2641 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2646 def __init__(self
, **kwargs
):
2647 jordan_product
= lambda x
,y
: x
2648 inner_product
= lambda x
,y
: 0
2651 # New defaults for keyword arguments
2652 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2653 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2655 super(TrivialEJA
, self
).__init
__(basis
,
2659 # The rank is zero using my definition, namely the dimension of the
2660 # largest subalgebra generated by any element.
2661 self
.rank
.set_cache(0)
2662 self
.one
.set_cache( self
.zero() )
2665 def random_instance(cls
, **kwargs
):
2666 # We don't take a "size" argument so the superclass method is
2667 # inappropriate for us.
2668 return cls(**kwargs
)
2671 class DirectSumEJA(FiniteDimensionalEJA
):
2673 The external (orthogonal) direct sum of two other Euclidean Jordan
2674 algebras. Essentially the Cartesian product of its two factors.
2675 Every Euclidean Jordan algebra decomposes into an orthogonal
2676 direct sum of simple Euclidean Jordan algebras, so no generality
2677 is lost by providing only this construction.
2681 sage: from mjo.eja.eja_algebra import (random_eja,
2683 ....: RealSymmetricEJA,
2688 sage: J1 = HadamardEJA(2)
2689 sage: J2 = RealSymmetricEJA(3)
2690 sage: J = DirectSumEJA(J1,J2)
2695 sage: J.matrix_space()
2696 The Cartesian product of (Full MatrixSpace of 2 by 1 dense matrices
2697 over Algebraic Real Field, Full MatrixSpace of 3 by 3 dense matrices
2698 over Algebraic Real Field)
2702 The external direct sum construction is only valid when the two factors
2703 have the same base ring; an error is raised otherwise::
2705 sage: set_random_seed()
2706 sage: J1 = random_eja(field=AA)
2707 sage: J2 = random_eja(field=QQ,orthonormalize=False)
2708 sage: J = DirectSumEJA(J1,J2)
2709 Traceback (most recent call last):
2711 ValueError: algebras must share the same base field
2714 def __init__(self
, J1
, J2
, **kwargs
):
2715 if J1
.base_ring() != J2
.base_ring():
2716 raise ValueError("algebras must share the same base field")
2717 field
= J1
.base_ring()
2719 M
= J1
.matrix_space().cartesian_product(J2
.matrix_space())
2720 self
._cartprod
_algebra
= J1
.cartesian_product(J2
)
2722 self
._matrix
_basis
= tuple( [M((a
,0)) for a
in J1
.matrix_basis()] +
2723 [M((0,b
)) for b
in J2
.matrix_basis()] )
2725 n
= len(self
._matrix
_basis
)
2727 CombinatorialFreeModule
.__init
__(
2731 category
=self
._cartprod
_algebra
.category(),
2734 self
.rank
.set_cache(J1
.rank() + J2
.rank())
2738 def product(self
,x
,y
):
2742 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
2743 ....: ComplexHermitianEJA,
2748 sage: set_random_seed()
2749 sage: J1 = JordanSpinEJA(3, field=QQ)
2750 sage: J2 = ComplexHermitianEJA(2, field=QQ, orthonormalize=False)
2751 sage: J = DirectSumEJA(J1,J2)
2752 sage: J.random_element()*J.random_element() in J
2756 xv
= self
._cartprod
_algebra
.from_vector(x
.to_vector())
2757 yv
= self
._cartprod
_algebra
.from_vector(y
.to_vector())
2758 return self
.from_vector((xv
*yv
).to_vector())
2761 def cartesian_factors(self
):
2763 Return the pair of this algebra's factors.
2767 sage: from mjo.eja.eja_algebra import (HadamardEJA,
2768 ....: JordanSpinEJA,
2773 sage: J1 = HadamardEJA(2, field=QQ)
2774 sage: J2 = JordanSpinEJA(3, field=QQ)
2775 sage: J = DirectSumEJA(J1,J2)
2776 sage: J.cartesian_factors()
2777 (Euclidean Jordan algebra of dimension 2 over Rational Field,
2778 Euclidean Jordan algebra of dimension 3 over Rational Field)
2781 return self
._cartprod
_algebra
.cartesian_factors()
2784 # def projections(self):
2786 # Return a pair of projections onto this algebra's factors.
2790 # sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
2791 # ....: ComplexHermitianEJA,
2792 # ....: DirectSumEJA)
2796 # sage: J1 = JordanSpinEJA(2)
2797 # sage: J2 = ComplexHermitianEJA(2)
2798 # sage: J = DirectSumEJA(J1,J2)
2799 # sage: (pi_left, pi_right) = J.projections()
2800 # sage: J.one().to_vector()
2801 # (1, 0, 1, 0, 0, 1)
2802 # sage: pi_left(J.one()).to_vector()
2804 # sage: pi_right(J.one()).to_vector()
2808 # (J1,J2) = self.factors()
2809 # m = J1.dimension()
2810 # n = J2.dimension()
2811 # V_basis = self.vector_space().basis()
2812 # # Need to specify the dimensions explicitly so that we don't
2813 # # wind up with a zero-by-zero matrix when we want e.g. a
2814 # # zero-by-two matrix (important for composing things).
2815 # P1 = matrix(self.base_ring(), m, m+n, V_basis[:m])
2816 # P2 = matrix(self.base_ring(), n, m+n, V_basis[m:])
2817 # pi_left = FiniteDimensionalEJAOperator(self,J1,P1)
2818 # pi_right = FiniteDimensionalEJAOperator(self,J2,P2)
2819 # return (pi_left, pi_right)
2821 # def inclusions(self):
2823 # Return the pair of inclusion maps from our factors into us.
2827 # sage: from mjo.eja.eja_algebra import (random_eja,
2828 # ....: JordanSpinEJA,
2829 # ....: RealSymmetricEJA,
2830 # ....: DirectSumEJA)
2834 # sage: J1 = JordanSpinEJA(3)
2835 # sage: J2 = RealSymmetricEJA(2)
2836 # sage: J = DirectSumEJA(J1,J2)
2837 # sage: (iota_left, iota_right) = J.inclusions()
2838 # sage: iota_left(J1.zero()) == J.zero()
2840 # sage: iota_right(J2.zero()) == J.zero()
2842 # sage: J1.one().to_vector()
2844 # sage: iota_left(J1.one()).to_vector()
2845 # (1, 0, 0, 0, 0, 0)
2846 # sage: J2.one().to_vector()
2848 # sage: iota_right(J2.one()).to_vector()
2849 # (0, 0, 0, 1, 0, 1)
2850 # sage: J.one().to_vector()
2851 # (1, 0, 0, 1, 0, 1)
2855 # Composing a projection with the corresponding inclusion should
2856 # produce the identity map, and mismatching them should produce
2859 # sage: set_random_seed()
2860 # sage: J1 = random_eja()
2861 # sage: J2 = random_eja()
2862 # sage: J = DirectSumEJA(J1,J2)
2863 # sage: (iota_left, iota_right) = J.inclusions()
2864 # sage: (pi_left, pi_right) = J.projections()
2865 # sage: pi_left*iota_left == J1.one().operator()
2867 # sage: pi_right*iota_right == J2.one().operator()
2869 # sage: (pi_left*iota_right).is_zero()
2871 # sage: (pi_right*iota_left).is_zero()
2875 # (J1,J2) = self.factors()
2876 # m = J1.dimension()
2877 # n = J2.dimension()
2878 # V_basis = self.vector_space().basis()
2879 # # Need to specify the dimensions explicitly so that we don't
2880 # # wind up with a zero-by-zero matrix when we want e.g. a
2881 # # two-by-zero matrix (important for composing things).
2882 # I1 = matrix.column(self.base_ring(), m, m+n, V_basis[:m])
2883 # I2 = matrix.column(self.base_ring(), n, m+n, V_basis[m:])
2884 # iota_left = FiniteDimensionalEJAOperator(J1,self,I1)
2885 # iota_right = FiniteDimensionalEJAOperator(J2,self,I2)
2886 # return (iota_left, iota_right)
2888 # def inner_product(self, x, y):
2890 # The standard Cartesian inner-product.
2892 # We project ``x`` and ``y`` onto our factors, and add up the
2893 # inner-products from the subalgebras.
2898 # sage: from mjo.eja.eja_algebra import (HadamardEJA,
2899 # ....: QuaternionHermitianEJA,
2900 # ....: DirectSumEJA)
2904 # sage: J1 = HadamardEJA(3,field=QQ)
2905 # sage: J2 = QuaternionHermitianEJA(2,field=QQ,orthonormalize=False)
2906 # sage: J = DirectSumEJA(J1,J2)
2907 # sage: x1 = J1.one()
2909 # sage: y1 = J2.one()
2911 # sage: x1.inner_product(x2)
2913 # sage: y1.inner_product(y2)
2915 # sage: J.one().inner_product(J.one())
2919 # (pi_left, pi_right) = self.projections()
2925 # return (x1.inner_product(y1) + x2.inner_product(y2))
2929 random_eja
= ConcreteEJA
.random_instance