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()
1724 a
= z
.list()[0] # real part, I guess
1725 b
= z
.list()[1] # imag part, I guess
1726 blocks
.append(matrix(field
, 2, [[a
,b
],[-b
,a
]]))
1728 return matrix
.block(field
, n
, blocks
)
1732 def real_unembed(cls
,M
):
1734 The inverse of _embed_complex_matrix().
1738 sage: from mjo.eja.eja_algebra import ComplexMatrixEJA
1742 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1743 ....: [-2, 1, -4, 3],
1744 ....: [ 9, 10, 11, 12],
1745 ....: [-10, 9, -12, 11] ])
1746 sage: ComplexMatrixEJA.real_unembed(A)
1748 [ 10*I + 9 12*I + 11]
1752 Unembedding is the inverse of embedding::
1754 sage: set_random_seed()
1755 sage: F = QuadraticField(-1, 'I')
1756 sage: M = random_matrix(F, 3)
1757 sage: Me = ComplexMatrixEJA.real_embed(M)
1758 sage: ComplexMatrixEJA.real_unembed(Me) == M
1762 super(ComplexMatrixEJA
,cls
).real_unembed(M
)
1764 d
= cls
.dimension_over_reals()
1766 # If "M" was normalized, its base ring might have roots
1767 # adjoined and they can stick around after unembedding.
1768 field
= M
.base_ring()
1769 R
= PolynomialRing(field
, 'z')
1772 # Sage doesn't know how to adjoin the complex "i" (the root of
1773 # x^2 + 1) to a field in a general way. Here, we just enumerate
1774 # all of the cases that I have cared to support so far.
1776 # Sage doesn't know how to embed AA into QQbar, i.e. how
1777 # to adjoin sqrt(-1) to AA.
1779 elif not field
.is_exact():
1781 F
= field
.complex_field()
1783 # Works for QQ and... maybe some other fields.
1784 F
= field
.extension(z
**2 + 1, 'I', embedding
=CLF(-1).sqrt())
1787 # Go top-left to bottom-right (reading order), converting every
1788 # 2-by-2 block we see to a single complex element.
1790 for k
in range(n
/d
):
1791 for j
in range(n
/d
):
1792 submat
= M
[d
*k
:d
*k
+d
,d
*j
:d
*j
+d
]
1793 if submat
[0,0] != submat
[1,1]:
1794 raise ValueError('bad on-diagonal submatrix')
1795 if submat
[0,1] != -submat
[1,0]:
1796 raise ValueError('bad off-diagonal submatrix')
1797 z
= submat
[0,0] + submat
[0,1]*i
1800 return matrix(F
, n
/d
, elements
)
1803 class ComplexHermitianEJA(ConcreteEJA
, ComplexMatrixEJA
):
1805 The rank-n simple EJA consisting of complex Hermitian n-by-n
1806 matrices over the real numbers, the usual symmetric Jordan product,
1807 and the real-part-of-trace inner product. It has dimension `n^2` over
1812 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1816 In theory, our "field" can be any subfield of the reals::
1818 sage: ComplexHermitianEJA(2, field=RDF)
1819 Euclidean Jordan algebra of dimension 4 over Real Double Field
1820 sage: ComplexHermitianEJA(2, field=RR)
1821 Euclidean Jordan algebra of dimension 4 over Real Field with
1822 53 bits of precision
1826 The dimension of this algebra is `n^2`::
1828 sage: set_random_seed()
1829 sage: n_max = ComplexHermitianEJA._max_random_instance_size()
1830 sage: n = ZZ.random_element(1, n_max)
1831 sage: J = ComplexHermitianEJA(n)
1832 sage: J.dimension() == n^2
1835 The Jordan multiplication is what we think it is::
1837 sage: set_random_seed()
1838 sage: J = ComplexHermitianEJA.random_instance()
1839 sage: x,y = J.random_elements(2)
1840 sage: actual = (x*y).to_matrix()
1841 sage: X = x.to_matrix()
1842 sage: Y = y.to_matrix()
1843 sage: expected = (X*Y + Y*X)/2
1844 sage: actual == expected
1846 sage: J(expected) == x*y
1849 We can change the generator prefix::
1851 sage: ComplexHermitianEJA(2, prefix='z').gens()
1854 We can construct the (trivial) algebra of rank zero::
1856 sage: ComplexHermitianEJA(0)
1857 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1862 def _denormalized_basis(cls
, n
):
1864 Returns a basis for the space of complex Hermitian n-by-n matrices.
1866 Why do we embed these? Basically, because all of numerical linear
1867 algebra assumes that you're working with vectors consisting of `n`
1868 entries from a field and scalars from the same field. There's no way
1869 to tell SageMath that (for example) the vectors contain complex
1870 numbers, while the scalar field is real.
1874 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1878 sage: set_random_seed()
1879 sage: n = ZZ.random_element(1,5)
1880 sage: B = ComplexHermitianEJA._denormalized_basis(n)
1881 sage: all( M.is_symmetric() for M in B)
1886 R
= PolynomialRing(field
, 'z')
1888 F
= field
.extension(z
**2 + 1, 'I')
1891 # This is like the symmetric case, but we need to be careful:
1893 # * We want conjugate-symmetry, not just symmetry.
1894 # * The diagonal will (as a result) be real.
1897 Eij
= matrix
.zero(F
,n
)
1899 for j
in range(i
+1):
1903 Sij
= cls
.real_embed(Eij
)
1906 # The second one has a minus because it's conjugated.
1907 Eij
[j
,i
] = 1 # Eij = Eij + Eij.transpose()
1908 Sij_real
= cls
.real_embed(Eij
)
1910 # Eij = I*Eij - I*Eij.transpose()
1913 Sij_imag
= cls
.real_embed(Eij
)
1919 # Since we embedded these, we can drop back to the "field" that we
1920 # started with instead of the complex extension "F".
1921 return tuple( s
.change_ring(field
) for s
in S
)
1924 def __init__(self
, n
, **kwargs
):
1925 # We know this is a valid EJA, but will double-check
1926 # if the user passes check_axioms=True.
1927 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
1929 super(ComplexHermitianEJA
, self
).__init
__(self
._denormalized
_basis
(n
),
1930 self
.jordan_product
,
1931 self
.trace_inner_product
,
1933 # TODO: this could be factored out somehow, but is left here
1934 # because the MatrixEJA is not presently a subclass of the
1935 # FDEJA class that defines rank() and one().
1936 self
.rank
.set_cache(n
)
1937 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
1938 self
.one
.set_cache(self(idV
))
1941 def _max_random_instance_size():
1942 return 3 # Dimension 9
1945 def random_instance(cls
, **kwargs
):
1947 Return a random instance of this type of algebra.
1949 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1950 return cls(n
, **kwargs
)
1952 class QuaternionMatrixEJA(MatrixEJA
):
1954 def dimension_over_reals():
1958 def real_embed(cls
,M
):
1960 Embed the n-by-n quaternion matrix ``M`` into the space of real
1961 matrices of size 4n-by-4n by first sending each quaternion entry `z
1962 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
1963 c+di],[-c + di, a-bi]]`, and then embedding those into a real
1968 sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA
1972 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1973 sage: i,j,k = Q.gens()
1974 sage: x = 1 + 2*i + 3*j + 4*k
1975 sage: M = matrix(Q, 1, [[x]])
1976 sage: QuaternionMatrixEJA.real_embed(M)
1982 Embedding is a homomorphism (isomorphism, in fact)::
1984 sage: set_random_seed()
1985 sage: n = ZZ.random_element(2)
1986 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1987 sage: X = random_matrix(Q, n)
1988 sage: Y = random_matrix(Q, n)
1989 sage: Xe = QuaternionMatrixEJA.real_embed(X)
1990 sage: Ye = QuaternionMatrixEJA.real_embed(Y)
1991 sage: XYe = QuaternionMatrixEJA.real_embed(X*Y)
1996 super(QuaternionMatrixEJA
,cls
).real_embed(M
)
1997 quaternions
= M
.base_ring()
2000 F
= QuadraticField(-1, 'I')
2005 t
= z
.coefficient_tuple()
2010 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
2011 [-c
+ d
*i
, a
- b
*i
]])
2012 realM
= ComplexMatrixEJA
.real_embed(cplxM
)
2013 blocks
.append(realM
)
2015 # We should have real entries by now, so use the realest field
2016 # we've got for the return value.
2017 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
2022 def real_unembed(cls
,M
):
2024 The inverse of _embed_quaternion_matrix().
2028 sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA
2032 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
2033 ....: [-2, 1, -4, 3],
2034 ....: [-3, 4, 1, -2],
2035 ....: [-4, -3, 2, 1]])
2036 sage: QuaternionMatrixEJA.real_unembed(M)
2037 [1 + 2*i + 3*j + 4*k]
2041 Unembedding is the inverse of embedding::
2043 sage: set_random_seed()
2044 sage: Q = QuaternionAlgebra(QQ, -1, -1)
2045 sage: M = random_matrix(Q, 3)
2046 sage: Me = QuaternionMatrixEJA.real_embed(M)
2047 sage: QuaternionMatrixEJA.real_unembed(Me) == M
2051 super(QuaternionMatrixEJA
,cls
).real_unembed(M
)
2053 d
= cls
.dimension_over_reals()
2055 # Use the base ring of the matrix to ensure that its entries can be
2056 # multiplied by elements of the quaternion algebra.
2057 field
= M
.base_ring()
2058 Q
= QuaternionAlgebra(field
,-1,-1)
2061 # Go top-left to bottom-right (reading order), converting every
2062 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
2065 for l
in range(n
/d
):
2066 for m
in range(n
/d
):
2067 submat
= ComplexMatrixEJA
.real_unembed(
2068 M
[d
*l
:d
*l
+d
,d
*m
:d
*m
+d
] )
2069 if submat
[0,0] != submat
[1,1].conjugate():
2070 raise ValueError('bad on-diagonal submatrix')
2071 if submat
[0,1] != -submat
[1,0].conjugate():
2072 raise ValueError('bad off-diagonal submatrix')
2073 z
= submat
[0,0].real()
2074 z
+= submat
[0,0].imag()*i
2075 z
+= submat
[0,1].real()*j
2076 z
+= submat
[0,1].imag()*k
2079 return matrix(Q
, n
/d
, elements
)
2082 class QuaternionHermitianEJA(ConcreteEJA
, QuaternionMatrixEJA
):
2084 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
2085 matrices, the usual symmetric Jordan product, and the
2086 real-part-of-trace inner product. It has dimension `2n^2 - n` over
2091 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2095 In theory, our "field" can be any subfield of the reals::
2097 sage: QuaternionHermitianEJA(2, field=RDF)
2098 Euclidean Jordan algebra of dimension 6 over Real Double Field
2099 sage: QuaternionHermitianEJA(2, field=RR)
2100 Euclidean Jordan algebra of dimension 6 over Real Field with
2101 53 bits of precision
2105 The dimension of this algebra is `2*n^2 - n`::
2107 sage: set_random_seed()
2108 sage: n_max = QuaternionHermitianEJA._max_random_instance_size()
2109 sage: n = ZZ.random_element(1, n_max)
2110 sage: J = QuaternionHermitianEJA(n)
2111 sage: J.dimension() == 2*(n^2) - n
2114 The Jordan multiplication is what we think it is::
2116 sage: set_random_seed()
2117 sage: J = QuaternionHermitianEJA.random_instance()
2118 sage: x,y = J.random_elements(2)
2119 sage: actual = (x*y).to_matrix()
2120 sage: X = x.to_matrix()
2121 sage: Y = y.to_matrix()
2122 sage: expected = (X*Y + Y*X)/2
2123 sage: actual == expected
2125 sage: J(expected) == x*y
2128 We can change the generator prefix::
2130 sage: QuaternionHermitianEJA(2, prefix='a').gens()
2131 (a0, a1, a2, a3, a4, a5)
2133 We can construct the (trivial) algebra of rank zero::
2135 sage: QuaternionHermitianEJA(0)
2136 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2140 def _denormalized_basis(cls
, n
):
2142 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
2144 Why do we embed these? Basically, because all of numerical
2145 linear algebra assumes that you're working with vectors consisting
2146 of `n` entries from a field and scalars from the same field. There's
2147 no way to tell SageMath that (for example) the vectors contain
2148 complex numbers, while the scalar field is real.
2152 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2156 sage: set_random_seed()
2157 sage: n = ZZ.random_element(1,5)
2158 sage: B = QuaternionHermitianEJA._denormalized_basis(n)
2159 sage: all( M.is_symmetric() for M in B )
2164 Q
= QuaternionAlgebra(QQ
,-1,-1)
2167 # This is like the symmetric case, but we need to be careful:
2169 # * We want conjugate-symmetry, not just symmetry.
2170 # * The diagonal will (as a result) be real.
2173 Eij
= matrix
.zero(Q
,n
)
2175 for j
in range(i
+1):
2179 Sij
= cls
.real_embed(Eij
)
2182 # The second, third, and fourth ones have a minus
2183 # because they're conjugated.
2184 # Eij = Eij + Eij.transpose()
2186 Sij_real
= cls
.real_embed(Eij
)
2188 # Eij = I*(Eij - Eij.transpose())
2191 Sij_I
= cls
.real_embed(Eij
)
2193 # Eij = J*(Eij - Eij.transpose())
2196 Sij_J
= cls
.real_embed(Eij
)
2198 # Eij = K*(Eij - Eij.transpose())
2201 Sij_K
= cls
.real_embed(Eij
)
2207 # Since we embedded these, we can drop back to the "field" that we
2208 # started with instead of the quaternion algebra "Q".
2209 return tuple( s
.change_ring(field
) for s
in S
)
2212 def __init__(self
, n
, **kwargs
):
2213 # We know this is a valid EJA, but will double-check
2214 # if the user passes check_axioms=True.
2215 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2217 super(QuaternionHermitianEJA
, self
).__init
__(self
._denormalized
_basis
(n
),
2218 self
.jordan_product
,
2219 self
.trace_inner_product
,
2221 # TODO: this could be factored out somehow, but is left here
2222 # because the MatrixEJA is not presently a subclass of the
2223 # FDEJA class that defines rank() and one().
2224 self
.rank
.set_cache(n
)
2225 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
2226 self
.one
.set_cache(self(idV
))
2230 def _max_random_instance_size():
2232 The maximum rank of a random QuaternionHermitianEJA.
2234 return 2 # Dimension 6
2237 def random_instance(cls
, **kwargs
):
2239 Return a random instance of this type of algebra.
2241 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2242 return cls(n
, **kwargs
)
2245 class HadamardEJA(ConcreteEJA
):
2247 Return the Euclidean Jordan Algebra corresponding to the set
2248 `R^n` under the Hadamard product.
2250 Note: this is nothing more than the Cartesian product of ``n``
2251 copies of the spin algebra. Once Cartesian product algebras
2252 are implemented, this can go.
2256 sage: from mjo.eja.eja_algebra import HadamardEJA
2260 This multiplication table can be verified by hand::
2262 sage: J = HadamardEJA(3)
2263 sage: e0,e1,e2 = J.gens()
2279 We can change the generator prefix::
2281 sage: HadamardEJA(3, prefix='r').gens()
2285 def __init__(self
, n
, **kwargs
):
2287 jordan_product
= lambda x
,y
: x
2288 inner_product
= lambda x
,y
: x
2290 def jordan_product(x
,y
):
2292 return P( xi
*yi
for (xi
,yi
) in zip(x
,y
) )
2294 def inner_product(x
,y
):
2297 # New defaults for keyword arguments. Don't orthonormalize
2298 # because our basis is already orthonormal with respect to our
2299 # inner-product. Don't check the axioms, because we know this
2300 # is a valid EJA... but do double-check if the user passes
2301 # check_axioms=True. Note: we DON'T override the "check_field"
2302 # default here, because the user can pass in a field!
2303 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2304 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2306 column_basis
= tuple( b
.column() for b
in FreeModule(ZZ
, n
).basis() )
2307 super().__init
__(column_basis
, jordan_product
, inner_product
, **kwargs
)
2308 self
.rank
.set_cache(n
)
2311 self
.one
.set_cache( self
.zero() )
2313 self
.one
.set_cache( sum(self
.gens()) )
2316 def _max_random_instance_size():
2318 The maximum dimension of a random HadamardEJA.
2323 def random_instance(cls
, **kwargs
):
2325 Return a random instance of this type of algebra.
2327 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2328 return cls(n
, **kwargs
)
2331 class BilinearFormEJA(ConcreteEJA
):
2333 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2334 with the half-trace inner product and jordan product ``x*y =
2335 (<Bx,y>,y_bar>, x0*y_bar + y0*x_bar)`` where `B = 1 \times B22` is
2336 a symmetric positive-definite "bilinear form" matrix. Its
2337 dimension is the size of `B`, and it has rank two in dimensions
2338 larger than two. It reduces to the ``JordanSpinEJA`` when `B` is
2339 the identity matrix of order ``n``.
2341 We insist that the one-by-one upper-left identity block of `B` be
2342 passed in as well so that we can be passed a matrix of size zero
2343 to construct a trivial algebra.
2347 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
2348 ....: JordanSpinEJA)
2352 When no bilinear form is specified, the identity matrix is used,
2353 and the resulting algebra is the Jordan spin algebra::
2355 sage: B = matrix.identity(AA,3)
2356 sage: J0 = BilinearFormEJA(B)
2357 sage: J1 = JordanSpinEJA(3)
2358 sage: J0.multiplication_table() == J0.multiplication_table()
2361 An error is raised if the matrix `B` does not correspond to a
2362 positive-definite bilinear form::
2364 sage: B = matrix.random(QQ,2,3)
2365 sage: J = BilinearFormEJA(B)
2366 Traceback (most recent call last):
2368 ValueError: bilinear form is not positive-definite
2369 sage: B = matrix.zero(QQ,3)
2370 sage: J = BilinearFormEJA(B)
2371 Traceback (most recent call last):
2373 ValueError: bilinear form is not positive-definite
2377 We can create a zero-dimensional algebra::
2379 sage: B = matrix.identity(AA,0)
2380 sage: J = BilinearFormEJA(B)
2384 We can check the multiplication condition given in the Jordan, von
2385 Neumann, and Wigner paper (and also discussed on my "On the
2386 symmetry..." paper). Note that this relies heavily on the standard
2387 choice of basis, as does anything utilizing the bilinear form
2388 matrix. We opt not to orthonormalize the basis, because if we
2389 did, we would have to normalize the `s_{i}` in a similar manner::
2391 sage: set_random_seed()
2392 sage: n = ZZ.random_element(5)
2393 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2394 sage: B11 = matrix.identity(QQ,1)
2395 sage: B22 = M.transpose()*M
2396 sage: B = block_matrix(2,2,[ [B11,0 ],
2398 sage: J = BilinearFormEJA(B, orthonormalize=False)
2399 sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
2400 sage: V = J.vector_space()
2401 sage: sis = [ J( V([0] + (M.inverse()*ei).list()).column() )
2402 ....: for ei in eis ]
2403 sage: actual = [ sis[i]*sis[j]
2404 ....: for i in range(n-1)
2405 ....: for j in range(n-1) ]
2406 sage: expected = [ J.one() if i == j else J.zero()
2407 ....: for i in range(n-1)
2408 ....: for j in range(n-1) ]
2409 sage: actual == expected
2413 def __init__(self
, B
, **kwargs
):
2414 # The matrix "B" is supplied by the user in most cases,
2415 # so it makes sense to check whether or not its positive-
2416 # definite unless we are specifically asked not to...
2417 if ("check_axioms" not in kwargs
) or kwargs
["check_axioms"]:
2418 if not B
.is_positive_definite():
2419 raise ValueError("bilinear form is not positive-definite")
2421 # However, all of the other data for this EJA is computed
2422 # by us in manner that guarantees the axioms are
2423 # satisfied. So, again, unless we are specifically asked to
2424 # verify things, we'll skip the rest of the checks.
2425 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2427 def inner_product(x
,y
):
2428 return (y
.T
*B
*x
)[0,0]
2430 def jordan_product(x
,y
):
2436 z0
= inner_product(y
,x
)
2437 zbar
= y0
*xbar
+ x0
*ybar
2438 return P([z0
] + zbar
.list())
2441 column_basis
= tuple( b
.column() for b
in FreeModule(ZZ
, n
).basis() )
2442 super(BilinearFormEJA
, self
).__init
__(column_basis
,
2447 # The rank of this algebra is two, unless we're in a
2448 # one-dimensional ambient space (because the rank is bounded
2449 # by the ambient dimension).
2450 self
.rank
.set_cache(min(n
,2))
2453 self
.one
.set_cache( self
.zero() )
2455 self
.one
.set_cache( self
.monomial(0) )
2458 def _max_random_instance_size():
2460 The maximum dimension of a random BilinearFormEJA.
2465 def random_instance(cls
, **kwargs
):
2467 Return a random instance of this algebra.
2469 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2471 B
= matrix
.identity(ZZ
, n
)
2472 return cls(B
, **kwargs
)
2474 B11
= matrix
.identity(ZZ
, 1)
2475 M
= matrix
.random(ZZ
, n
-1)
2476 I
= matrix
.identity(ZZ
, n
-1)
2478 while alpha
.is_zero():
2479 alpha
= ZZ
.random_element().abs()
2480 B22
= M
.transpose()*M
+ alpha
*I
2482 from sage
.matrix
.special
import block_matrix
2483 B
= block_matrix(2,2, [ [B11
, ZZ(0) ],
2486 return cls(B
, **kwargs
)
2489 class JordanSpinEJA(BilinearFormEJA
):
2491 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2492 with the usual inner product and jordan product ``x*y =
2493 (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2498 sage: from mjo.eja.eja_algebra import JordanSpinEJA
2502 This multiplication table can be verified by hand::
2504 sage: J = JordanSpinEJA(4)
2505 sage: e0,e1,e2,e3 = J.gens()
2521 We can change the generator prefix::
2523 sage: JordanSpinEJA(2, prefix='B').gens()
2528 Ensure that we have the usual inner product on `R^n`::
2530 sage: set_random_seed()
2531 sage: J = JordanSpinEJA.random_instance()
2532 sage: x,y = J.random_elements(2)
2533 sage: actual = x.inner_product(y)
2534 sage: expected = x.to_vector().inner_product(y.to_vector())
2535 sage: actual == expected
2539 def __init__(self
, n
, **kwargs
):
2540 # This is a special case of the BilinearFormEJA with the
2541 # identity matrix as its bilinear form.
2542 B
= matrix
.identity(ZZ
, n
)
2544 # Don't orthonormalize because our basis is already
2545 # orthonormal with respect to our inner-product.
2546 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2548 # But also don't pass check_field=False here, because the user
2549 # can pass in a field!
2550 super(JordanSpinEJA
, self
).__init
__(B
, **kwargs
)
2553 def _max_random_instance_size():
2555 The maximum dimension of a random JordanSpinEJA.
2560 def random_instance(cls
, **kwargs
):
2562 Return a random instance of this type of algebra.
2564 Needed here to override the implementation for ``BilinearFormEJA``.
2566 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2567 return cls(n
, **kwargs
)
2570 class TrivialEJA(ConcreteEJA
):
2572 The trivial Euclidean Jordan algebra consisting of only a zero element.
2576 sage: from mjo.eja.eja_algebra import TrivialEJA
2580 sage: J = TrivialEJA()
2587 sage: 7*J.one()*12*J.one()
2589 sage: J.one().inner_product(J.one())
2591 sage: J.one().norm()
2593 sage: J.one().subalgebra_generated_by()
2594 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2599 def __init__(self
, **kwargs
):
2600 jordan_product
= lambda x
,y
: x
2601 inner_product
= lambda x
,y
: 0
2604 # New defaults for keyword arguments
2605 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2606 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2608 super(TrivialEJA
, self
).__init
__(basis
,
2612 # The rank is zero using my definition, namely the dimension of the
2613 # largest subalgebra generated by any element.
2614 self
.rank
.set_cache(0)
2615 self
.one
.set_cache( self
.zero() )
2618 def random_instance(cls
, **kwargs
):
2619 # We don't take a "size" argument so the superclass method is
2620 # inappropriate for us.
2621 return cls(**kwargs
)
2623 # class DirectSumEJA(ConcreteEJA):
2625 # The external (orthogonal) direct sum of two other Euclidean Jordan
2626 # algebras. Essentially the Cartesian product of its two factors.
2627 # Every Euclidean Jordan algebra decomposes into an orthogonal
2628 # direct sum of simple Euclidean Jordan algebras, so no generality
2629 # is lost by providing only this construction.
2633 # sage: from mjo.eja.eja_algebra import (random_eja,
2634 # ....: HadamardEJA,
2635 # ....: RealSymmetricEJA,
2636 # ....: DirectSumEJA)
2640 # sage: J1 = HadamardEJA(2)
2641 # sage: J2 = RealSymmetricEJA(3)
2642 # sage: J = DirectSumEJA(J1,J2)
2643 # sage: J.dimension()
2650 # The external direct sum construction is only valid when the two factors
2651 # have the same base ring; an error is raised otherwise::
2653 # sage: set_random_seed()
2654 # sage: J1 = random_eja(field=AA)
2655 # sage: J2 = random_eja(field=QQ,orthonormalize=False)
2656 # sage: J = DirectSumEJA(J1,J2)
2657 # Traceback (most recent call last):
2659 # ValueError: algebras must share the same base field
2662 # def __init__(self, J1, J2, **kwargs):
2663 # if J1.base_ring() != J2.base_ring():
2664 # raise ValueError("algebras must share the same base field")
2665 # field = J1.base_ring()
2667 # self._factors = (J1, J2)
2668 # n1 = J1.dimension()
2669 # n2 = J2.dimension()
2671 # V = VectorSpace(field, n)
2672 # mult_table = [ [ V.zero() for j in range(i+1) ]
2673 # for i in range(n) ]
2674 # for i in range(n1):
2675 # for j in range(i+1):
2676 # p = (J1.monomial(i)*J1.monomial(j)).to_vector()
2677 # mult_table[i][j] = V(p.list() + [field.zero()]*n2)
2679 # for i in range(n2):
2680 # for j in range(i+1):
2681 # p = (J2.monomial(i)*J2.monomial(j)).to_vector()
2682 # mult_table[n1+i][n1+j] = V([field.zero()]*n1 + p.list())
2684 # # TODO: build the IP table here from the two constituent IP
2685 # # matrices (it'll be block diagonal, I think).
2686 # ip_table = [ [ field.zero() for j in range(i+1) ]
2687 # for i in range(n) ]
2688 # super(DirectSumEJA, self).__init__(field,
2691 # check_axioms=False,
2693 # self.rank.set_cache(J1.rank() + J2.rank())
2696 # def factors(self):
2698 # Return the pair of this algebra's factors.
2702 # sage: from mjo.eja.eja_algebra import (HadamardEJA,
2703 # ....: JordanSpinEJA,
2704 # ....: DirectSumEJA)
2708 # sage: J1 = HadamardEJA(2, field=QQ)
2709 # sage: J2 = JordanSpinEJA(3, field=QQ)
2710 # sage: J = DirectSumEJA(J1,J2)
2712 # (Euclidean Jordan algebra of dimension 2 over Rational Field,
2713 # Euclidean Jordan algebra of dimension 3 over Rational Field)
2716 # return self._factors
2718 # def projections(self):
2720 # Return a pair of projections onto this algebra's factors.
2724 # sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
2725 # ....: ComplexHermitianEJA,
2726 # ....: DirectSumEJA)
2730 # sage: J1 = JordanSpinEJA(2)
2731 # sage: J2 = ComplexHermitianEJA(2)
2732 # sage: J = DirectSumEJA(J1,J2)
2733 # sage: (pi_left, pi_right) = J.projections()
2734 # sage: J.one().to_vector()
2735 # (1, 0, 1, 0, 0, 1)
2736 # sage: pi_left(J.one()).to_vector()
2738 # sage: pi_right(J.one()).to_vector()
2742 # (J1,J2) = self.factors()
2743 # m = J1.dimension()
2744 # n = J2.dimension()
2745 # V_basis = self.vector_space().basis()
2746 # # Need to specify the dimensions explicitly so that we don't
2747 # # wind up with a zero-by-zero matrix when we want e.g. a
2748 # # zero-by-two matrix (important for composing things).
2749 # P1 = matrix(self.base_ring(), m, m+n, V_basis[:m])
2750 # P2 = matrix(self.base_ring(), n, m+n, V_basis[m:])
2751 # pi_left = FiniteDimensionalEJAOperator(self,J1,P1)
2752 # pi_right = FiniteDimensionalEJAOperator(self,J2,P2)
2753 # return (pi_left, pi_right)
2755 # def inclusions(self):
2757 # Return the pair of inclusion maps from our factors into us.
2761 # sage: from mjo.eja.eja_algebra import (random_eja,
2762 # ....: JordanSpinEJA,
2763 # ....: RealSymmetricEJA,
2764 # ....: DirectSumEJA)
2768 # sage: J1 = JordanSpinEJA(3)
2769 # sage: J2 = RealSymmetricEJA(2)
2770 # sage: J = DirectSumEJA(J1,J2)
2771 # sage: (iota_left, iota_right) = J.inclusions()
2772 # sage: iota_left(J1.zero()) == J.zero()
2774 # sage: iota_right(J2.zero()) == J.zero()
2776 # sage: J1.one().to_vector()
2778 # sage: iota_left(J1.one()).to_vector()
2779 # (1, 0, 0, 0, 0, 0)
2780 # sage: J2.one().to_vector()
2782 # sage: iota_right(J2.one()).to_vector()
2783 # (0, 0, 0, 1, 0, 1)
2784 # sage: J.one().to_vector()
2785 # (1, 0, 0, 1, 0, 1)
2789 # Composing a projection with the corresponding inclusion should
2790 # produce the identity map, and mismatching them should produce
2793 # sage: set_random_seed()
2794 # sage: J1 = random_eja()
2795 # sage: J2 = random_eja()
2796 # sage: J = DirectSumEJA(J1,J2)
2797 # sage: (iota_left, iota_right) = J.inclusions()
2798 # sage: (pi_left, pi_right) = J.projections()
2799 # sage: pi_left*iota_left == J1.one().operator()
2801 # sage: pi_right*iota_right == J2.one().operator()
2803 # sage: (pi_left*iota_right).is_zero()
2805 # sage: (pi_right*iota_left).is_zero()
2809 # (J1,J2) = self.factors()
2810 # m = J1.dimension()
2811 # n = J2.dimension()
2812 # V_basis = self.vector_space().basis()
2813 # # Need to specify the dimensions explicitly so that we don't
2814 # # wind up with a zero-by-zero matrix when we want e.g. a
2815 # # two-by-zero matrix (important for composing things).
2816 # I1 = matrix.column(self.base_ring(), m, m+n, V_basis[:m])
2817 # I2 = matrix.column(self.base_ring(), n, m+n, V_basis[m:])
2818 # iota_left = FiniteDimensionalEJAOperator(J1,self,I1)
2819 # iota_right = FiniteDimensionalEJAOperator(J2,self,I2)
2820 # return (iota_left, iota_right)
2822 # def inner_product(self, x, y):
2824 # The standard Cartesian inner-product.
2826 # We project ``x`` and ``y`` onto our factors, and add up the
2827 # inner-products from the subalgebras.
2832 # sage: from mjo.eja.eja_algebra import (HadamardEJA,
2833 # ....: QuaternionHermitianEJA,
2834 # ....: DirectSumEJA)
2838 # sage: J1 = HadamardEJA(3,field=QQ)
2839 # sage: J2 = QuaternionHermitianEJA(2,field=QQ,orthonormalize=False)
2840 # sage: J = DirectSumEJA(J1,J2)
2841 # sage: x1 = J1.one()
2843 # sage: y1 = J2.one()
2845 # sage: x1.inner_product(x2)
2847 # sage: y1.inner_product(y2)
2849 # sage: J.one().inner_product(J.one())
2853 # (pi_left, pi_right) = self.projections()
2859 # return (x1.inner_product(y1) + x2.inner_product(y2))
2863 random_eja
= ConcreteEJA
.random_instance