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].matrix_space()
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
1105 n
= self
.dimension()
1106 R
= self
.coordinate_polynomial_ring()
1108 F
= R
.fraction_field()
1111 # From a result in my book, these are the entries of the
1112 # basis representation of L_x.
1113 return sum( vars[k
]*self
.monomial(k
).operator().matrix()[i
,j
]
1116 L_x
= matrix(F
, n
, n
, L_x_i_j
)
1119 if self
.rank
.is_in_cache():
1121 # There's no need to pad the system with redundant
1122 # columns if we *know* they'll be redundant.
1125 # Compute an extra power in case the rank is equal to
1126 # the dimension (otherwise, we would stop at x^(r-1)).
1127 x_powers
= [ (L_x
**k
)*self
.one().to_vector()
1128 for k
in range(n
+1) ]
1129 A
= matrix
.column(F
, x_powers
[:n
])
1130 AE
= A
.extended_echelon_form()
1137 # The theory says that only the first "r" coefficients are
1138 # nonzero, and they actually live in the original polynomial
1139 # ring and not the fraction field. We negate them because
1140 # in the actual characteristic polynomial, they get moved
1141 # to the other side where x^r lives.
1142 return -A_rref
.solve_right(E
*b
).change_ring(R
)[:r
]
1147 Return the rank of this EJA.
1149 This is a cached method because we know the rank a priori for
1150 all of the algebras we can construct. Thus we can avoid the
1151 expensive ``_charpoly_coefficients()`` call unless we truly
1152 need to compute the whole characteristic polynomial.
1156 sage: from mjo.eja.eja_algebra import (HadamardEJA,
1157 ....: JordanSpinEJA,
1158 ....: RealSymmetricEJA,
1159 ....: ComplexHermitianEJA,
1160 ....: QuaternionHermitianEJA,
1165 The rank of the Jordan spin algebra is always two::
1167 sage: JordanSpinEJA(2).rank()
1169 sage: JordanSpinEJA(3).rank()
1171 sage: JordanSpinEJA(4).rank()
1174 The rank of the `n`-by-`n` Hermitian real, complex, or
1175 quaternion matrices is `n`::
1177 sage: RealSymmetricEJA(4).rank()
1179 sage: ComplexHermitianEJA(3).rank()
1181 sage: QuaternionHermitianEJA(2).rank()
1186 Ensure that every EJA that we know how to construct has a
1187 positive integer rank, unless the algebra is trivial in
1188 which case its rank will be zero::
1190 sage: set_random_seed()
1191 sage: J = random_eja()
1195 sage: r > 0 or (r == 0 and J.is_trivial())
1198 Ensure that computing the rank actually works, since the ranks
1199 of all simple algebras are known and will be cached by default::
1201 sage: set_random_seed() # long time
1202 sage: J = random_eja() # long time
1203 sage: caches = J.rank() # long time
1204 sage: J.rank.clear_cache() # long time
1205 sage: J.rank() == cached # long time
1209 return len(self
._charpoly
_coefficients
())
1212 def vector_space(self
):
1214 Return the vector space that underlies this algebra.
1218 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1222 sage: J = RealSymmetricEJA(2)
1223 sage: J.vector_space()
1224 Vector space of dimension 3 over...
1227 return self
.zero().to_vector().parent().ambient_vector_space()
1230 Element
= FiniteDimensionalEJAElement
1232 class RationalBasisEJA(FiniteDimensionalEJA
):
1234 New class for algebras whose supplied basis elements have all rational entries.
1238 sage: from mjo.eja.eja_algebra import BilinearFormEJA
1242 The supplied basis is orthonormalized by default::
1244 sage: B = matrix(QQ, [[1, 0, 0], [0, 25, -32], [0, -32, 41]])
1245 sage: J = BilinearFormEJA(B)
1246 sage: J.matrix_basis()
1263 # Abuse the check_field parameter to check that the entries of
1264 # out basis (in ambient coordinates) are in the field QQ.
1265 if not all( all(b_i
in QQ
for b_i
in b
.list()) for b
in basis
):
1266 raise TypeError("basis not rational")
1268 self
._rational
_algebra
= None
1270 # There's no point in constructing the extra algebra if this
1271 # one is already rational.
1273 # Note: the same Jordan and inner-products work here,
1274 # because they are necessarily defined with respect to
1275 # ambient coordinates and not any particular basis.
1276 self
._rational
_algebra
= FiniteDimensionalEJA(
1281 orthonormalize
=False,
1285 super().__init
__(basis
,
1289 check_field
=check_field
,
1293 def _charpoly_coefficients(self
):
1297 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
1298 ....: JordanSpinEJA)
1302 The base ring of the resulting polynomial coefficients is what
1303 it should be, and not the rationals (unless the algebra was
1304 already over the rationals)::
1306 sage: J = JordanSpinEJA(3)
1307 sage: J._charpoly_coefficients()
1308 (X1^2 - X2^2 - X3^2, -2*X1)
1309 sage: a0 = J._charpoly_coefficients()[0]
1311 Algebraic Real Field
1312 sage: a0.base_ring()
1313 Algebraic Real Field
1316 if self
._rational
_algebra
is None:
1317 # There's no need to construct *another* algebra over the
1318 # rationals if this one is already over the
1319 # rationals. Likewise, if we never orthonormalized our
1320 # basis, we might as well just use the given one.
1321 return super()._charpoly
_coefficients
()
1323 # Do the computation over the rationals. The answer will be
1324 # the same, because all we've done is a change of basis.
1325 # Then, change back from QQ to our real base ring
1326 a
= ( a_i
.change_ring(self
.base_ring())
1327 for a_i
in self
._rational
_algebra
._charpoly
_coefficients
() )
1329 if self
._deortho
_matrix
is None:
1330 # This can happen if our base ring was, say, AA and we
1331 # chose not to (or didn't need to) orthonormalize. It's
1332 # still faster to do the computations over QQ even if
1333 # the numbers in the boxes stay the same.
1336 # Otherwise, convert the coordinate variables back to the
1337 # deorthonormalized ones.
1338 R
= self
.coordinate_polynomial_ring()
1339 from sage
.modules
.free_module_element
import vector
1340 X
= vector(R
, R
.gens())
1341 BX
= self
._deortho
_matrix
*X
1343 subs_dict
= { X[i]: BX[i] for i in range(len(X)) }
1344 return tuple( a_i
.subs(subs_dict
) for a_i
in a
)
1346 class ConcreteEJA(RationalBasisEJA
):
1348 A class for the Euclidean Jordan algebras that we know by name.
1350 These are the Jordan algebras whose basis, multiplication table,
1351 rank, and so on are known a priori. More to the point, they are
1352 the Euclidean Jordan algebras for which we are able to conjure up
1353 a "random instance."
1357 sage: from mjo.eja.eja_algebra import ConcreteEJA
1361 Our basis is normalized with respect to the algebra's inner
1362 product, unless we specify otherwise::
1364 sage: set_random_seed()
1365 sage: J = ConcreteEJA.random_instance()
1366 sage: all( b.norm() == 1 for b in J.gens() )
1369 Since our basis is orthonormal with respect to the algebra's inner
1370 product, and since we know that this algebra is an EJA, any
1371 left-multiplication operator's matrix will be symmetric because
1372 natural->EJA basis representation is an isometry and within the
1373 EJA the operator is self-adjoint by the Jordan axiom::
1375 sage: set_random_seed()
1376 sage: J = ConcreteEJA.random_instance()
1377 sage: x = J.random_element()
1378 sage: x.operator().is_self_adjoint()
1383 def _max_random_instance_size():
1385 Return an integer "size" that is an upper bound on the size of
1386 this algebra when it is used in a random test
1387 case. Unfortunately, the term "size" is ambiguous -- when
1388 dealing with `R^n` under either the Hadamard or Jordan spin
1389 product, the "size" refers to the dimension `n`. When dealing
1390 with a matrix algebra (real symmetric or complex/quaternion
1391 Hermitian), it refers to the size of the matrix, which is far
1392 less than the dimension of the underlying vector space.
1394 This method must be implemented in each subclass.
1396 raise NotImplementedError
1399 def random_instance(cls
, *args
, **kwargs
):
1401 Return a random instance of this type of algebra.
1403 This method should be implemented in each subclass.
1405 from sage
.misc
.prandom
import choice
1406 eja_class
= choice(cls
.__subclasses
__())
1408 # These all bubble up to the RationalBasisEJA superclass
1409 # constructor, so any (kw)args valid there are also valid
1411 return eja_class
.random_instance(*args
, **kwargs
)
1416 def dimension_over_reals():
1418 The dimension of this matrix's base ring over the reals.
1420 The reals are dimension one over themselves, obviously; that's
1421 just `\mathbb{R}^{1}`. Likewise, the complex numbers `a + bi`
1422 have dimension two. Finally, the quaternions have dimension
1423 four over the reals.
1425 This is used to determine the size of the matrix returned from
1426 :meth:`real_embed`, among other things.
1428 raise NotImplementedError
1431 def real_embed(cls
,M
):
1433 Embed the matrix ``M`` into a space of real matrices.
1435 The matrix ``M`` can have entries in any field at the moment:
1436 the real numbers, complex numbers, or quaternions. And although
1437 they are not a field, we can probably support octonions at some
1438 point, too. This function returns a real matrix that "acts like"
1439 the original with respect to matrix multiplication; i.e.
1441 real_embed(M*N) = real_embed(M)*real_embed(N)
1444 if M
.ncols() != M
.nrows():
1445 raise ValueError("the matrix 'M' must be square")
1450 def real_unembed(cls
,M
):
1452 The inverse of :meth:`real_embed`.
1454 if M
.ncols() != M
.nrows():
1455 raise ValueError("the matrix 'M' must be square")
1456 if not ZZ(M
.nrows()).mod(cls
.dimension_over_reals()).is_zero():
1457 raise ValueError("the matrix 'M' must be a real embedding")
1461 def jordan_product(X
,Y
):
1462 return (X
*Y
+ Y
*X
)/2
1465 def trace_inner_product(cls
,X
,Y
):
1467 Compute the trace inner-product of two real-embeddings.
1471 sage: from mjo.eja.eja_algebra import (RealSymmetricEJA,
1472 ....: ComplexHermitianEJA,
1473 ....: QuaternionHermitianEJA)
1477 This gives the same answer as it would if we computed the trace
1478 from the unembedded (original) matrices::
1480 sage: set_random_seed()
1481 sage: J = RealSymmetricEJA.random_instance()
1482 sage: x,y = J.random_elements(2)
1483 sage: Xe = x.to_matrix()
1484 sage: Ye = y.to_matrix()
1485 sage: X = J.real_unembed(Xe)
1486 sage: Y = J.real_unembed(Ye)
1487 sage: expected = (X*Y).trace()
1488 sage: actual = J.trace_inner_product(Xe,Ye)
1489 sage: actual == expected
1494 sage: set_random_seed()
1495 sage: J = ComplexHermitianEJA.random_instance()
1496 sage: x,y = J.random_elements(2)
1497 sage: Xe = x.to_matrix()
1498 sage: Ye = y.to_matrix()
1499 sage: X = J.real_unembed(Xe)
1500 sage: Y = J.real_unembed(Ye)
1501 sage: expected = (X*Y).trace().real()
1502 sage: actual = J.trace_inner_product(Xe,Ye)
1503 sage: actual == expected
1508 sage: set_random_seed()
1509 sage: J = QuaternionHermitianEJA.random_instance()
1510 sage: x,y = J.random_elements(2)
1511 sage: Xe = x.to_matrix()
1512 sage: Ye = y.to_matrix()
1513 sage: X = J.real_unembed(Xe)
1514 sage: Y = J.real_unembed(Ye)
1515 sage: expected = (X*Y).trace().coefficient_tuple()[0]
1516 sage: actual = J.trace_inner_product(Xe,Ye)
1517 sage: actual == expected
1521 Xu
= cls
.real_unembed(X
)
1522 Yu
= cls
.real_unembed(Y
)
1523 tr
= (Xu
*Yu
).trace()
1526 # Works in QQ, AA, RDF, et cetera.
1528 except AttributeError:
1529 # A quaternion doesn't have a real() method, but does
1530 # have coefficient_tuple() method that returns the
1531 # coefficients of 1, i, j, and k -- in that order.
1532 return tr
.coefficient_tuple()[0]
1535 class RealMatrixEJA(MatrixEJA
):
1537 def dimension_over_reals():
1541 class RealSymmetricEJA(ConcreteEJA
, RealMatrixEJA
):
1543 The rank-n simple EJA consisting of real symmetric n-by-n
1544 matrices, the usual symmetric Jordan product, and the trace inner
1545 product. It has dimension `(n^2 + n)/2` over the reals.
1549 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1553 sage: J = RealSymmetricEJA(2)
1554 sage: e0, e1, e2 = J.gens()
1562 In theory, our "field" can be any subfield of the reals::
1564 sage: RealSymmetricEJA(2, field=RDF)
1565 Euclidean Jordan algebra of dimension 3 over Real Double Field
1566 sage: RealSymmetricEJA(2, field=RR)
1567 Euclidean Jordan algebra of dimension 3 over Real Field with
1568 53 bits of precision
1572 The dimension of this algebra is `(n^2 + n) / 2`::
1574 sage: set_random_seed()
1575 sage: n_max = RealSymmetricEJA._max_random_instance_size()
1576 sage: n = ZZ.random_element(1, n_max)
1577 sage: J = RealSymmetricEJA(n)
1578 sage: J.dimension() == (n^2 + n)/2
1581 The Jordan multiplication is what we think it is::
1583 sage: set_random_seed()
1584 sage: J = RealSymmetricEJA.random_instance()
1585 sage: x,y = J.random_elements(2)
1586 sage: actual = (x*y).to_matrix()
1587 sage: X = x.to_matrix()
1588 sage: Y = y.to_matrix()
1589 sage: expected = (X*Y + Y*X)/2
1590 sage: actual == expected
1592 sage: J(expected) == x*y
1595 We can change the generator prefix::
1597 sage: RealSymmetricEJA(3, prefix='q').gens()
1598 (q0, q1, q2, q3, q4, q5)
1600 We can construct the (trivial) algebra of rank zero::
1602 sage: RealSymmetricEJA(0)
1603 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1607 def _denormalized_basis(cls
, n
):
1609 Return a basis for the space of real symmetric n-by-n matrices.
1613 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1617 sage: set_random_seed()
1618 sage: n = ZZ.random_element(1,5)
1619 sage: B = RealSymmetricEJA._denormalized_basis(n)
1620 sage: all( M.is_symmetric() for M in B)
1624 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1628 for j
in range(i
+1):
1629 Eij
= matrix(ZZ
, n
, lambda k
,l
: k
==i
and l
==j
)
1633 Sij
= Eij
+ Eij
.transpose()
1639 def _max_random_instance_size():
1640 return 4 # Dimension 10
1643 def random_instance(cls
, **kwargs
):
1645 Return a random instance of this type of algebra.
1647 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1648 return cls(n
, **kwargs
)
1650 def __init__(self
, n
, **kwargs
):
1651 # We know this is a valid EJA, but will double-check
1652 # if the user passes check_axioms=True.
1653 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
1655 super(RealSymmetricEJA
, self
).__init
__(self
._denormalized
_basis
(n
),
1656 self
.jordan_product
,
1657 self
.trace_inner_product
,
1660 # TODO: this could be factored out somehow, but is left here
1661 # because the MatrixEJA is not presently a subclass of the
1662 # FDEJA class that defines rank() and one().
1663 self
.rank
.set_cache(n
)
1664 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
1665 self
.one
.set_cache(self(idV
))
1669 class ComplexMatrixEJA(MatrixEJA
):
1671 def dimension_over_reals():
1675 def real_embed(cls
,M
):
1677 Embed the n-by-n complex matrix ``M`` into the space of real
1678 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1679 bi` to the block matrix ``[[a,b],[-b,a]]``.
1683 sage: from mjo.eja.eja_algebra import ComplexMatrixEJA
1687 sage: F = QuadraticField(-1, 'I')
1688 sage: x1 = F(4 - 2*i)
1689 sage: x2 = F(1 + 2*i)
1692 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1693 sage: ComplexMatrixEJA.real_embed(M)
1702 Embedding is a homomorphism (isomorphism, in fact)::
1704 sage: set_random_seed()
1705 sage: n = ZZ.random_element(3)
1706 sage: F = QuadraticField(-1, 'I')
1707 sage: X = random_matrix(F, n)
1708 sage: Y = random_matrix(F, n)
1709 sage: Xe = ComplexMatrixEJA.real_embed(X)
1710 sage: Ye = ComplexMatrixEJA.real_embed(Y)
1711 sage: XYe = ComplexMatrixEJA.real_embed(X*Y)
1716 super(ComplexMatrixEJA
,cls
).real_embed(M
)
1719 # We don't need any adjoined elements...
1720 field
= M
.base_ring().base_ring()
1726 blocks
.append(matrix(field
, 2, [ [ a
, b
],
1729 return matrix
.block(field
, n
, blocks
)
1733 def real_unembed(cls
,M
):
1735 The inverse of _embed_complex_matrix().
1739 sage: from mjo.eja.eja_algebra import ComplexMatrixEJA
1743 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1744 ....: [-2, 1, -4, 3],
1745 ....: [ 9, 10, 11, 12],
1746 ....: [-10, 9, -12, 11] ])
1747 sage: ComplexMatrixEJA.real_unembed(A)
1749 [ 10*I + 9 12*I + 11]
1753 Unembedding is the inverse of embedding::
1755 sage: set_random_seed()
1756 sage: F = QuadraticField(-1, 'I')
1757 sage: M = random_matrix(F, 3)
1758 sage: Me = ComplexMatrixEJA.real_embed(M)
1759 sage: ComplexMatrixEJA.real_unembed(Me) == M
1763 super(ComplexMatrixEJA
,cls
).real_unembed(M
)
1765 d
= cls
.dimension_over_reals()
1767 # If "M" was normalized, its base ring might have roots
1768 # adjoined and they can stick around after unembedding.
1769 field
= M
.base_ring()
1770 R
= PolynomialRing(field
, 'z')
1773 # Sage doesn't know how to adjoin the complex "i" (the root of
1774 # x^2 + 1) to a field in a general way. Here, we just enumerate
1775 # all of the cases that I have cared to support so far.
1777 # Sage doesn't know how to embed AA into QQbar, i.e. how
1778 # to adjoin sqrt(-1) to AA.
1780 elif not field
.is_exact():
1782 F
= field
.complex_field()
1784 # Works for QQ and... maybe some other fields.
1785 F
= field
.extension(z
**2 + 1, 'I', embedding
=CLF(-1).sqrt())
1788 # Go top-left to bottom-right (reading order), converting every
1789 # 2-by-2 block we see to a single complex element.
1791 for k
in range(n
/d
):
1792 for j
in range(n
/d
):
1793 submat
= M
[d
*k
:d
*k
+d
,d
*j
:d
*j
+d
]
1794 if submat
[0,0] != submat
[1,1]:
1795 raise ValueError('bad on-diagonal submatrix')
1796 if submat
[0,1] != -submat
[1,0]:
1797 raise ValueError('bad off-diagonal submatrix')
1798 z
= submat
[0,0] + submat
[0,1]*i
1801 return matrix(F
, n
/d
, elements
)
1804 class ComplexHermitianEJA(ConcreteEJA
, ComplexMatrixEJA
):
1806 The rank-n simple EJA consisting of complex Hermitian n-by-n
1807 matrices over the real numbers, the usual symmetric Jordan product,
1808 and the real-part-of-trace inner product. It has dimension `n^2` over
1813 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1817 In theory, our "field" can be any subfield of the reals::
1819 sage: ComplexHermitianEJA(2, field=RDF)
1820 Euclidean Jordan algebra of dimension 4 over Real Double Field
1821 sage: ComplexHermitianEJA(2, field=RR)
1822 Euclidean Jordan algebra of dimension 4 over Real Field with
1823 53 bits of precision
1827 The dimension of this algebra is `n^2`::
1829 sage: set_random_seed()
1830 sage: n_max = ComplexHermitianEJA._max_random_instance_size()
1831 sage: n = ZZ.random_element(1, n_max)
1832 sage: J = ComplexHermitianEJA(n)
1833 sage: J.dimension() == n^2
1836 The Jordan multiplication is what we think it is::
1838 sage: set_random_seed()
1839 sage: J = ComplexHermitianEJA.random_instance()
1840 sage: x,y = J.random_elements(2)
1841 sage: actual = (x*y).to_matrix()
1842 sage: X = x.to_matrix()
1843 sage: Y = y.to_matrix()
1844 sage: expected = (X*Y + Y*X)/2
1845 sage: actual == expected
1847 sage: J(expected) == x*y
1850 We can change the generator prefix::
1852 sage: ComplexHermitianEJA(2, prefix='z').gens()
1855 We can construct the (trivial) algebra of rank zero::
1857 sage: ComplexHermitianEJA(0)
1858 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1863 def _denormalized_basis(cls
, n
):
1865 Returns a basis for the space of complex Hermitian n-by-n matrices.
1867 Why do we embed these? Basically, because all of numerical linear
1868 algebra assumes that you're working with vectors consisting of `n`
1869 entries from a field and scalars from the same field. There's no way
1870 to tell SageMath that (for example) the vectors contain complex
1871 numbers, while the scalar field is real.
1875 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1879 sage: set_random_seed()
1880 sage: n = ZZ.random_element(1,5)
1881 sage: B = ComplexHermitianEJA._denormalized_basis(n)
1882 sage: all( M.is_symmetric() for M in B)
1887 R
= PolynomialRing(field
, 'z')
1889 F
= field
.extension(z
**2 + 1, 'I')
1892 # This is like the symmetric case, but we need to be careful:
1894 # * We want conjugate-symmetry, not just symmetry.
1895 # * The diagonal will (as a result) be real.
1898 Eij
= matrix
.zero(F
,n
)
1900 for j
in range(i
+1):
1904 Sij
= cls
.real_embed(Eij
)
1907 # The second one has a minus because it's conjugated.
1908 Eij
[j
,i
] = 1 # Eij = Eij + Eij.transpose()
1909 Sij_real
= cls
.real_embed(Eij
)
1911 # Eij = I*Eij - I*Eij.transpose()
1914 Sij_imag
= cls
.real_embed(Eij
)
1920 # Since we embedded these, we can drop back to the "field" that we
1921 # started with instead of the complex extension "F".
1922 return tuple( s
.change_ring(field
) for s
in S
)
1925 def __init__(self
, n
, **kwargs
):
1926 # We know this is a valid EJA, but will double-check
1927 # if the user passes check_axioms=True.
1928 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
1930 super(ComplexHermitianEJA
, self
).__init
__(self
._denormalized
_basis
(n
),
1931 self
.jordan_product
,
1932 self
.trace_inner_product
,
1934 # TODO: this could be factored out somehow, but is left here
1935 # because the MatrixEJA is not presently a subclass of the
1936 # FDEJA class that defines rank() and one().
1937 self
.rank
.set_cache(n
)
1938 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
1939 self
.one
.set_cache(self(idV
))
1942 def _max_random_instance_size():
1943 return 3 # Dimension 9
1946 def random_instance(cls
, **kwargs
):
1948 Return a random instance of this type of algebra.
1950 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1951 return cls(n
, **kwargs
)
1953 class QuaternionMatrixEJA(MatrixEJA
):
1955 def dimension_over_reals():
1959 def real_embed(cls
,M
):
1961 Embed the n-by-n quaternion matrix ``M`` into the space of real
1962 matrices of size 4n-by-4n by first sending each quaternion entry `z
1963 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
1964 c+di],[-c + di, a-bi]]`, and then embedding those into a real
1969 sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA
1973 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1974 sage: i,j,k = Q.gens()
1975 sage: x = 1 + 2*i + 3*j + 4*k
1976 sage: M = matrix(Q, 1, [[x]])
1977 sage: QuaternionMatrixEJA.real_embed(M)
1983 Embedding is a homomorphism (isomorphism, in fact)::
1985 sage: set_random_seed()
1986 sage: n = ZZ.random_element(2)
1987 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1988 sage: X = random_matrix(Q, n)
1989 sage: Y = random_matrix(Q, n)
1990 sage: Xe = QuaternionMatrixEJA.real_embed(X)
1991 sage: Ye = QuaternionMatrixEJA.real_embed(Y)
1992 sage: XYe = QuaternionMatrixEJA.real_embed(X*Y)
1997 super(QuaternionMatrixEJA
,cls
).real_embed(M
)
1998 quaternions
= M
.base_ring()
2001 F
= QuadraticField(-1, 'I')
2006 t
= z
.coefficient_tuple()
2011 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
2012 [-c
+ d
*i
, a
- b
*i
]])
2013 realM
= ComplexMatrixEJA
.real_embed(cplxM
)
2014 blocks
.append(realM
)
2016 # We should have real entries by now, so use the realest field
2017 # we've got for the return value.
2018 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
2023 def real_unembed(cls
,M
):
2025 The inverse of _embed_quaternion_matrix().
2029 sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA
2033 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
2034 ....: [-2, 1, -4, 3],
2035 ....: [-3, 4, 1, -2],
2036 ....: [-4, -3, 2, 1]])
2037 sage: QuaternionMatrixEJA.real_unembed(M)
2038 [1 + 2*i + 3*j + 4*k]
2042 Unembedding is the inverse of embedding::
2044 sage: set_random_seed()
2045 sage: Q = QuaternionAlgebra(QQ, -1, -1)
2046 sage: M = random_matrix(Q, 3)
2047 sage: Me = QuaternionMatrixEJA.real_embed(M)
2048 sage: QuaternionMatrixEJA.real_unembed(Me) == M
2052 super(QuaternionMatrixEJA
,cls
).real_unembed(M
)
2054 d
= cls
.dimension_over_reals()
2056 # Use the base ring of the matrix to ensure that its entries can be
2057 # multiplied by elements of the quaternion algebra.
2058 field
= M
.base_ring()
2059 Q
= QuaternionAlgebra(field
,-1,-1)
2062 # Go top-left to bottom-right (reading order), converting every
2063 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
2066 for l
in range(n
/d
):
2067 for m
in range(n
/d
):
2068 submat
= ComplexMatrixEJA
.real_unembed(
2069 M
[d
*l
:d
*l
+d
,d
*m
:d
*m
+d
] )
2070 if submat
[0,0] != submat
[1,1].conjugate():
2071 raise ValueError('bad on-diagonal submatrix')
2072 if submat
[0,1] != -submat
[1,0].conjugate():
2073 raise ValueError('bad off-diagonal submatrix')
2074 z
= submat
[0,0].real()
2075 z
+= submat
[0,0].imag()*i
2076 z
+= submat
[0,1].real()*j
2077 z
+= submat
[0,1].imag()*k
2080 return matrix(Q
, n
/d
, elements
)
2083 class QuaternionHermitianEJA(ConcreteEJA
, QuaternionMatrixEJA
):
2085 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
2086 matrices, the usual symmetric Jordan product, and the
2087 real-part-of-trace inner product. It has dimension `2n^2 - n` over
2092 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2096 In theory, our "field" can be any subfield of the reals::
2098 sage: QuaternionHermitianEJA(2, field=RDF)
2099 Euclidean Jordan algebra of dimension 6 over Real Double Field
2100 sage: QuaternionHermitianEJA(2, field=RR)
2101 Euclidean Jordan algebra of dimension 6 over Real Field with
2102 53 bits of precision
2106 The dimension of this algebra is `2*n^2 - n`::
2108 sage: set_random_seed()
2109 sage: n_max = QuaternionHermitianEJA._max_random_instance_size()
2110 sage: n = ZZ.random_element(1, n_max)
2111 sage: J = QuaternionHermitianEJA(n)
2112 sage: J.dimension() == 2*(n^2) - n
2115 The Jordan multiplication is what we think it is::
2117 sage: set_random_seed()
2118 sage: J = QuaternionHermitianEJA.random_instance()
2119 sage: x,y = J.random_elements(2)
2120 sage: actual = (x*y).to_matrix()
2121 sage: X = x.to_matrix()
2122 sage: Y = y.to_matrix()
2123 sage: expected = (X*Y + Y*X)/2
2124 sage: actual == expected
2126 sage: J(expected) == x*y
2129 We can change the generator prefix::
2131 sage: QuaternionHermitianEJA(2, prefix='a').gens()
2132 (a0, a1, a2, a3, a4, a5)
2134 We can construct the (trivial) algebra of rank zero::
2136 sage: QuaternionHermitianEJA(0)
2137 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2141 def _denormalized_basis(cls
, n
):
2143 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
2145 Why do we embed these? Basically, because all of numerical
2146 linear algebra assumes that you're working with vectors consisting
2147 of `n` entries from a field and scalars from the same field. There's
2148 no way to tell SageMath that (for example) the vectors contain
2149 complex numbers, while the scalar field is real.
2153 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2157 sage: set_random_seed()
2158 sage: n = ZZ.random_element(1,5)
2159 sage: B = QuaternionHermitianEJA._denormalized_basis(n)
2160 sage: all( M.is_symmetric() for M in B )
2165 Q
= QuaternionAlgebra(QQ
,-1,-1)
2168 # This is like the symmetric case, but we need to be careful:
2170 # * We want conjugate-symmetry, not just symmetry.
2171 # * The diagonal will (as a result) be real.
2174 Eij
= matrix
.zero(Q
,n
)
2176 for j
in range(i
+1):
2180 Sij
= cls
.real_embed(Eij
)
2183 # The second, third, and fourth ones have a minus
2184 # because they're conjugated.
2185 # Eij = Eij + Eij.transpose()
2187 Sij_real
= cls
.real_embed(Eij
)
2189 # Eij = I*(Eij - Eij.transpose())
2192 Sij_I
= cls
.real_embed(Eij
)
2194 # Eij = J*(Eij - Eij.transpose())
2197 Sij_J
= cls
.real_embed(Eij
)
2199 # Eij = K*(Eij - Eij.transpose())
2202 Sij_K
= cls
.real_embed(Eij
)
2208 # Since we embedded these, we can drop back to the "field" that we
2209 # started with instead of the quaternion algebra "Q".
2210 return tuple( s
.change_ring(field
) for s
in S
)
2213 def __init__(self
, n
, **kwargs
):
2214 # We know this is a valid EJA, but will double-check
2215 # if the user passes check_axioms=True.
2216 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2218 super(QuaternionHermitianEJA
, self
).__init
__(self
._denormalized
_basis
(n
),
2219 self
.jordan_product
,
2220 self
.trace_inner_product
,
2222 # TODO: this could be factored out somehow, but is left here
2223 # because the MatrixEJA is not presently a subclass of the
2224 # FDEJA class that defines rank() and one().
2225 self
.rank
.set_cache(n
)
2226 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
2227 self
.one
.set_cache(self(idV
))
2231 def _max_random_instance_size():
2233 The maximum rank of a random QuaternionHermitianEJA.
2235 return 2 # Dimension 6
2238 def random_instance(cls
, **kwargs
):
2240 Return a random instance of this type of algebra.
2242 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2243 return cls(n
, **kwargs
)
2246 class HadamardEJA(ConcreteEJA
):
2248 Return the Euclidean Jordan Algebra corresponding to the set
2249 `R^n` under the Hadamard product.
2251 Note: this is nothing more than the Cartesian product of ``n``
2252 copies of the spin algebra. Once Cartesian product algebras
2253 are implemented, this can go.
2257 sage: from mjo.eja.eja_algebra import HadamardEJA
2261 This multiplication table can be verified by hand::
2263 sage: J = HadamardEJA(3)
2264 sage: e0,e1,e2 = J.gens()
2280 We can change the generator prefix::
2282 sage: HadamardEJA(3, prefix='r').gens()
2286 def __init__(self
, n
, **kwargs
):
2288 jordan_product
= lambda x
,y
: x
2289 inner_product
= lambda x
,y
: x
2291 def jordan_product(x
,y
):
2293 return P( xi
*yi
for (xi
,yi
) in zip(x
,y
) )
2295 def inner_product(x
,y
):
2298 # New defaults for keyword arguments. Don't orthonormalize
2299 # because our basis is already orthonormal with respect to our
2300 # inner-product. Don't check the axioms, because we know this
2301 # is a valid EJA... but do double-check if the user passes
2302 # check_axioms=True. Note: we DON'T override the "check_field"
2303 # default here, because the user can pass in a field!
2304 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2305 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2307 column_basis
= tuple( b
.column() for b
in FreeModule(ZZ
, n
).basis() )
2308 super().__init
__(column_basis
, jordan_product
, inner_product
, **kwargs
)
2309 self
.rank
.set_cache(n
)
2312 self
.one
.set_cache( self
.zero() )
2314 self
.one
.set_cache( sum(self
.gens()) )
2317 def _max_random_instance_size():
2319 The maximum dimension of a random HadamardEJA.
2324 def random_instance(cls
, **kwargs
):
2326 Return a random instance of this type of algebra.
2328 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2329 return cls(n
, **kwargs
)
2332 class BilinearFormEJA(ConcreteEJA
):
2334 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2335 with the half-trace inner product and jordan product ``x*y =
2336 (<Bx,y>,y_bar>, x0*y_bar + y0*x_bar)`` where `B = 1 \times B22` is
2337 a symmetric positive-definite "bilinear form" matrix. Its
2338 dimension is the size of `B`, and it has rank two in dimensions
2339 larger than two. It reduces to the ``JordanSpinEJA`` when `B` is
2340 the identity matrix of order ``n``.
2342 We insist that the one-by-one upper-left identity block of `B` be
2343 passed in as well so that we can be passed a matrix of size zero
2344 to construct a trivial algebra.
2348 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
2349 ....: JordanSpinEJA)
2353 When no bilinear form is specified, the identity matrix is used,
2354 and the resulting algebra is the Jordan spin algebra::
2356 sage: B = matrix.identity(AA,3)
2357 sage: J0 = BilinearFormEJA(B)
2358 sage: J1 = JordanSpinEJA(3)
2359 sage: J0.multiplication_table() == J0.multiplication_table()
2362 An error is raised if the matrix `B` does not correspond to a
2363 positive-definite bilinear form::
2365 sage: B = matrix.random(QQ,2,3)
2366 sage: J = BilinearFormEJA(B)
2367 Traceback (most recent call last):
2369 ValueError: bilinear form is not positive-definite
2370 sage: B = matrix.zero(QQ,3)
2371 sage: J = BilinearFormEJA(B)
2372 Traceback (most recent call last):
2374 ValueError: bilinear form is not positive-definite
2378 We can create a zero-dimensional algebra::
2380 sage: B = matrix.identity(AA,0)
2381 sage: J = BilinearFormEJA(B)
2385 We can check the multiplication condition given in the Jordan, von
2386 Neumann, and Wigner paper (and also discussed on my "On the
2387 symmetry..." paper). Note that this relies heavily on the standard
2388 choice of basis, as does anything utilizing the bilinear form
2389 matrix. We opt not to orthonormalize the basis, because if we
2390 did, we would have to normalize the `s_{i}` in a similar manner::
2392 sage: set_random_seed()
2393 sage: n = ZZ.random_element(5)
2394 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2395 sage: B11 = matrix.identity(QQ,1)
2396 sage: B22 = M.transpose()*M
2397 sage: B = block_matrix(2,2,[ [B11,0 ],
2399 sage: J = BilinearFormEJA(B, orthonormalize=False)
2400 sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
2401 sage: V = J.vector_space()
2402 sage: sis = [ J( V([0] + (M.inverse()*ei).list()).column() )
2403 ....: for ei in eis ]
2404 sage: actual = [ sis[i]*sis[j]
2405 ....: for i in range(n-1)
2406 ....: for j in range(n-1) ]
2407 sage: expected = [ J.one() if i == j else J.zero()
2408 ....: for i in range(n-1)
2409 ....: for j in range(n-1) ]
2410 sage: actual == expected
2414 def __init__(self
, B
, **kwargs
):
2415 # The matrix "B" is supplied by the user in most cases,
2416 # so it makes sense to check whether or not its positive-
2417 # definite unless we are specifically asked not to...
2418 if ("check_axioms" not in kwargs
) or kwargs
["check_axioms"]:
2419 if not B
.is_positive_definite():
2420 raise ValueError("bilinear form is not positive-definite")
2422 # However, all of the other data for this EJA is computed
2423 # by us in manner that guarantees the axioms are
2424 # satisfied. So, again, unless we are specifically asked to
2425 # verify things, we'll skip the rest of the checks.
2426 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2428 def inner_product(x
,y
):
2429 return (y
.T
*B
*x
)[0,0]
2431 def jordan_product(x
,y
):
2437 z0
= inner_product(y
,x
)
2438 zbar
= y0
*xbar
+ x0
*ybar
2439 return P([z0
] + zbar
.list())
2442 column_basis
= tuple( b
.column() for b
in FreeModule(ZZ
, n
).basis() )
2443 super(BilinearFormEJA
, self
).__init
__(column_basis
,
2448 # The rank of this algebra is two, unless we're in a
2449 # one-dimensional ambient space (because the rank is bounded
2450 # by the ambient dimension).
2451 self
.rank
.set_cache(min(n
,2))
2454 self
.one
.set_cache( self
.zero() )
2456 self
.one
.set_cache( self
.monomial(0) )
2459 def _max_random_instance_size():
2461 The maximum dimension of a random BilinearFormEJA.
2466 def random_instance(cls
, **kwargs
):
2468 Return a random instance of this algebra.
2470 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2472 B
= matrix
.identity(ZZ
, n
)
2473 return cls(B
, **kwargs
)
2475 B11
= matrix
.identity(ZZ
, 1)
2476 M
= matrix
.random(ZZ
, n
-1)
2477 I
= matrix
.identity(ZZ
, n
-1)
2479 while alpha
.is_zero():
2480 alpha
= ZZ
.random_element().abs()
2481 B22
= M
.transpose()*M
+ alpha
*I
2483 from sage
.matrix
.special
import block_matrix
2484 B
= block_matrix(2,2, [ [B11
, ZZ(0) ],
2487 return cls(B
, **kwargs
)
2490 class JordanSpinEJA(BilinearFormEJA
):
2492 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2493 with the usual inner product and jordan product ``x*y =
2494 (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2499 sage: from mjo.eja.eja_algebra import JordanSpinEJA
2503 This multiplication table can be verified by hand::
2505 sage: J = JordanSpinEJA(4)
2506 sage: e0,e1,e2,e3 = J.gens()
2522 We can change the generator prefix::
2524 sage: JordanSpinEJA(2, prefix='B').gens()
2529 Ensure that we have the usual inner product on `R^n`::
2531 sage: set_random_seed()
2532 sage: J = JordanSpinEJA.random_instance()
2533 sage: x,y = J.random_elements(2)
2534 sage: actual = x.inner_product(y)
2535 sage: expected = x.to_vector().inner_product(y.to_vector())
2536 sage: actual == expected
2540 def __init__(self
, n
, **kwargs
):
2541 # This is a special case of the BilinearFormEJA with the
2542 # identity matrix as its bilinear form.
2543 B
= matrix
.identity(ZZ
, n
)
2545 # Don't orthonormalize because our basis is already
2546 # orthonormal with respect to our inner-product.
2547 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2549 # But also don't pass check_field=False here, because the user
2550 # can pass in a field!
2551 super(JordanSpinEJA
, self
).__init
__(B
, **kwargs
)
2554 def _max_random_instance_size():
2556 The maximum dimension of a random JordanSpinEJA.
2561 def random_instance(cls
, **kwargs
):
2563 Return a random instance of this type of algebra.
2565 Needed here to override the implementation for ``BilinearFormEJA``.
2567 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2568 return cls(n
, **kwargs
)
2571 class TrivialEJA(ConcreteEJA
):
2573 The trivial Euclidean Jordan algebra consisting of only a zero element.
2577 sage: from mjo.eja.eja_algebra import TrivialEJA
2581 sage: J = TrivialEJA()
2588 sage: 7*J.one()*12*J.one()
2590 sage: J.one().inner_product(J.one())
2592 sage: J.one().norm()
2594 sage: J.one().subalgebra_generated_by()
2595 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2600 def __init__(self
, **kwargs
):
2601 jordan_product
= lambda x
,y
: x
2602 inner_product
= lambda x
,y
: 0
2605 # New defaults for keyword arguments
2606 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2607 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2609 super(TrivialEJA
, self
).__init
__(basis
,
2613 # The rank is zero using my definition, namely the dimension of the
2614 # largest subalgebra generated by any element.
2615 self
.rank
.set_cache(0)
2616 self
.one
.set_cache( self
.zero() )
2619 def random_instance(cls
, **kwargs
):
2620 # We don't take a "size" argument so the superclass method is
2621 # inappropriate for us.
2622 return cls(**kwargs
)
2624 # class DirectSumEJA(ConcreteEJA):
2626 # The external (orthogonal) direct sum of two other Euclidean Jordan
2627 # algebras. Essentially the Cartesian product of its two factors.
2628 # Every Euclidean Jordan algebra decomposes into an orthogonal
2629 # direct sum of simple Euclidean Jordan algebras, so no generality
2630 # is lost by providing only this construction.
2634 # sage: from mjo.eja.eja_algebra import (random_eja,
2635 # ....: HadamardEJA,
2636 # ....: RealSymmetricEJA,
2637 # ....: DirectSumEJA)
2641 # sage: J1 = HadamardEJA(2)
2642 # sage: J2 = RealSymmetricEJA(3)
2643 # sage: J = DirectSumEJA(J1,J2)
2644 # sage: J.dimension()
2651 # The external direct sum construction is only valid when the two factors
2652 # have the same base ring; an error is raised otherwise::
2654 # sage: set_random_seed()
2655 # sage: J1 = random_eja(field=AA)
2656 # sage: J2 = random_eja(field=QQ,orthonormalize=False)
2657 # sage: J = DirectSumEJA(J1,J2)
2658 # Traceback (most recent call last):
2660 # ValueError: algebras must share the same base field
2663 # def __init__(self, J1, J2, **kwargs):
2664 # if J1.base_ring() != J2.base_ring():
2665 # raise ValueError("algebras must share the same base field")
2666 # field = J1.base_ring()
2668 # self._factors = (J1, J2)
2669 # n1 = J1.dimension()
2670 # n2 = J2.dimension()
2672 # V = VectorSpace(field, n)
2673 # mult_table = [ [ V.zero() for j in range(i+1) ]
2674 # for i in range(n) ]
2675 # for i in range(n1):
2676 # for j in range(i+1):
2677 # p = (J1.monomial(i)*J1.monomial(j)).to_vector()
2678 # mult_table[i][j] = V(p.list() + [field.zero()]*n2)
2680 # for i in range(n2):
2681 # for j in range(i+1):
2682 # p = (J2.monomial(i)*J2.monomial(j)).to_vector()
2683 # mult_table[n1+i][n1+j] = V([field.zero()]*n1 + p.list())
2685 # # TODO: build the IP table here from the two constituent IP
2686 # # matrices (it'll be block diagonal, I think).
2687 # ip_table = [ [ field.zero() for j in range(i+1) ]
2688 # for i in range(n) ]
2689 # super(DirectSumEJA, self).__init__(field,
2692 # check_axioms=False,
2694 # self.rank.set_cache(J1.rank() + J2.rank())
2697 # def factors(self):
2699 # Return the pair of this algebra's factors.
2703 # sage: from mjo.eja.eja_algebra import (HadamardEJA,
2704 # ....: JordanSpinEJA,
2705 # ....: DirectSumEJA)
2709 # sage: J1 = HadamardEJA(2, field=QQ)
2710 # sage: J2 = JordanSpinEJA(3, field=QQ)
2711 # sage: J = DirectSumEJA(J1,J2)
2713 # (Euclidean Jordan algebra of dimension 2 over Rational Field,
2714 # Euclidean Jordan algebra of dimension 3 over Rational Field)
2717 # return self._factors
2719 # def projections(self):
2721 # Return a pair of projections onto this algebra's factors.
2725 # sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
2726 # ....: ComplexHermitianEJA,
2727 # ....: DirectSumEJA)
2731 # sage: J1 = JordanSpinEJA(2)
2732 # sage: J2 = ComplexHermitianEJA(2)
2733 # sage: J = DirectSumEJA(J1,J2)
2734 # sage: (pi_left, pi_right) = J.projections()
2735 # sage: J.one().to_vector()
2736 # (1, 0, 1, 0, 0, 1)
2737 # sage: pi_left(J.one()).to_vector()
2739 # sage: pi_right(J.one()).to_vector()
2743 # (J1,J2) = self.factors()
2744 # m = J1.dimension()
2745 # n = J2.dimension()
2746 # V_basis = self.vector_space().basis()
2747 # # Need to specify the dimensions explicitly so that we don't
2748 # # wind up with a zero-by-zero matrix when we want e.g. a
2749 # # zero-by-two matrix (important for composing things).
2750 # P1 = matrix(self.base_ring(), m, m+n, V_basis[:m])
2751 # P2 = matrix(self.base_ring(), n, m+n, V_basis[m:])
2752 # pi_left = FiniteDimensionalEJAOperator(self,J1,P1)
2753 # pi_right = FiniteDimensionalEJAOperator(self,J2,P2)
2754 # return (pi_left, pi_right)
2756 # def inclusions(self):
2758 # Return the pair of inclusion maps from our factors into us.
2762 # sage: from mjo.eja.eja_algebra import (random_eja,
2763 # ....: JordanSpinEJA,
2764 # ....: RealSymmetricEJA,
2765 # ....: DirectSumEJA)
2769 # sage: J1 = JordanSpinEJA(3)
2770 # sage: J2 = RealSymmetricEJA(2)
2771 # sage: J = DirectSumEJA(J1,J2)
2772 # sage: (iota_left, iota_right) = J.inclusions()
2773 # sage: iota_left(J1.zero()) == J.zero()
2775 # sage: iota_right(J2.zero()) == J.zero()
2777 # sage: J1.one().to_vector()
2779 # sage: iota_left(J1.one()).to_vector()
2780 # (1, 0, 0, 0, 0, 0)
2781 # sage: J2.one().to_vector()
2783 # sage: iota_right(J2.one()).to_vector()
2784 # (0, 0, 0, 1, 0, 1)
2785 # sage: J.one().to_vector()
2786 # (1, 0, 0, 1, 0, 1)
2790 # Composing a projection with the corresponding inclusion should
2791 # produce the identity map, and mismatching them should produce
2794 # sage: set_random_seed()
2795 # sage: J1 = random_eja()
2796 # sage: J2 = random_eja()
2797 # sage: J = DirectSumEJA(J1,J2)
2798 # sage: (iota_left, iota_right) = J.inclusions()
2799 # sage: (pi_left, pi_right) = J.projections()
2800 # sage: pi_left*iota_left == J1.one().operator()
2802 # sage: pi_right*iota_right == J2.one().operator()
2804 # sage: (pi_left*iota_right).is_zero()
2806 # sage: (pi_right*iota_left).is_zero()
2810 # (J1,J2) = self.factors()
2811 # m = J1.dimension()
2812 # n = J2.dimension()
2813 # V_basis = self.vector_space().basis()
2814 # # Need to specify the dimensions explicitly so that we don't
2815 # # wind up with a zero-by-zero matrix when we want e.g. a
2816 # # two-by-zero matrix (important for composing things).
2817 # I1 = matrix.column(self.base_ring(), m, m+n, V_basis[:m])
2818 # I2 = matrix.column(self.base_ring(), n, m+n, V_basis[m:])
2819 # iota_left = FiniteDimensionalEJAOperator(J1,self,I1)
2820 # iota_right = FiniteDimensionalEJAOperator(J2,self,I2)
2821 # return (iota_left, iota_right)
2823 # def inner_product(self, x, y):
2825 # The standard Cartesian inner-product.
2827 # We project ``x`` and ``y`` onto our factors, and add up the
2828 # inner-products from the subalgebras.
2833 # sage: from mjo.eja.eja_algebra import (HadamardEJA,
2834 # ....: QuaternionHermitianEJA,
2835 # ....: DirectSumEJA)
2839 # sage: J1 = HadamardEJA(3,field=QQ)
2840 # sage: J2 = QuaternionHermitianEJA(2,field=QQ,orthonormalize=False)
2841 # sage: J = DirectSumEJA(J1,J2)
2842 # sage: x1 = J1.one()
2844 # sage: y1 = J2.one()
2846 # sage: x1.inner_product(x2)
2848 # sage: y1.inner_product(y2)
2850 # sage: J.one().inner_product(J.one())
2854 # (pi_left, pi_right) = self.projections()
2860 # return (x1.inner_product(y1) + x2.inner_product(y2))
2864 random_eja
= ConcreteEJA
.random_instance