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: field = QuadraticField(2, 'sqrt2')
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.
1899 for j
in range(i
+1):
1900 Eij
= matrix(F
, n
, lambda k
,l
: k
==i
and l
==j
)
1902 Sij
= cls
.real_embed(Eij
)
1905 # The second one has a minus because it's conjugated.
1906 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
1908 Sij_imag
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
1911 # Since we embedded these, we can drop back to the "field" that we
1912 # started with instead of the complex extension "F".
1913 return tuple( s
.change_ring(field
) for s
in S
)
1916 def __init__(self
, n
, **kwargs
):
1917 # We know this is a valid EJA, but will double-check
1918 # if the user passes check_axioms=True.
1919 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
1921 super(ComplexHermitianEJA
, self
).__init
__(self
._denormalized
_basis
(n
),
1922 self
.jordan_product
,
1923 self
.trace_inner_product
,
1925 # TODO: this could be factored out somehow, but is left here
1926 # because the MatrixEJA is not presently a subclass of the
1927 # FDEJA class that defines rank() and one().
1928 self
.rank
.set_cache(n
)
1929 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
1930 self
.one
.set_cache(self(idV
))
1933 def _max_random_instance_size():
1934 return 3 # Dimension 9
1937 def random_instance(cls
, **kwargs
):
1939 Return a random instance of this type of algebra.
1941 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1942 return cls(n
, **kwargs
)
1944 class QuaternionMatrixEJA(MatrixEJA
):
1946 def dimension_over_reals():
1950 def real_embed(cls
,M
):
1952 Embed the n-by-n quaternion matrix ``M`` into the space of real
1953 matrices of size 4n-by-4n by first sending each quaternion entry `z
1954 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
1955 c+di],[-c + di, a-bi]]`, and then embedding those into a real
1960 sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA
1964 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1965 sage: i,j,k = Q.gens()
1966 sage: x = 1 + 2*i + 3*j + 4*k
1967 sage: M = matrix(Q, 1, [[x]])
1968 sage: QuaternionMatrixEJA.real_embed(M)
1974 Embedding is a homomorphism (isomorphism, in fact)::
1976 sage: set_random_seed()
1977 sage: n = ZZ.random_element(2)
1978 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1979 sage: X = random_matrix(Q, n)
1980 sage: Y = random_matrix(Q, n)
1981 sage: Xe = QuaternionMatrixEJA.real_embed(X)
1982 sage: Ye = QuaternionMatrixEJA.real_embed(Y)
1983 sage: XYe = QuaternionMatrixEJA.real_embed(X*Y)
1988 super(QuaternionMatrixEJA
,cls
).real_embed(M
)
1989 quaternions
= M
.base_ring()
1992 F
= QuadraticField(-1, 'I')
1997 t
= z
.coefficient_tuple()
2002 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
2003 [-c
+ d
*i
, a
- b
*i
]])
2004 realM
= ComplexMatrixEJA
.real_embed(cplxM
)
2005 blocks
.append(realM
)
2007 # We should have real entries by now, so use the realest field
2008 # we've got for the return value.
2009 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
2014 def real_unembed(cls
,M
):
2016 The inverse of _embed_quaternion_matrix().
2020 sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA
2024 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
2025 ....: [-2, 1, -4, 3],
2026 ....: [-3, 4, 1, -2],
2027 ....: [-4, -3, 2, 1]])
2028 sage: QuaternionMatrixEJA.real_unembed(M)
2029 [1 + 2*i + 3*j + 4*k]
2033 Unembedding is the inverse of embedding::
2035 sage: set_random_seed()
2036 sage: Q = QuaternionAlgebra(QQ, -1, -1)
2037 sage: M = random_matrix(Q, 3)
2038 sage: Me = QuaternionMatrixEJA.real_embed(M)
2039 sage: QuaternionMatrixEJA.real_unembed(Me) == M
2043 super(QuaternionMatrixEJA
,cls
).real_unembed(M
)
2045 d
= cls
.dimension_over_reals()
2047 # Use the base ring of the matrix to ensure that its entries can be
2048 # multiplied by elements of the quaternion algebra.
2049 field
= M
.base_ring()
2050 Q
= QuaternionAlgebra(field
,-1,-1)
2053 # Go top-left to bottom-right (reading order), converting every
2054 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
2057 for l
in range(n
/d
):
2058 for m
in range(n
/d
):
2059 submat
= ComplexMatrixEJA
.real_unembed(
2060 M
[d
*l
:d
*l
+d
,d
*m
:d
*m
+d
] )
2061 if submat
[0,0] != submat
[1,1].conjugate():
2062 raise ValueError('bad on-diagonal submatrix')
2063 if submat
[0,1] != -submat
[1,0].conjugate():
2064 raise ValueError('bad off-diagonal submatrix')
2065 z
= submat
[0,0].real()
2066 z
+= submat
[0,0].imag()*i
2067 z
+= submat
[0,1].real()*j
2068 z
+= submat
[0,1].imag()*k
2071 return matrix(Q
, n
/d
, elements
)
2074 class QuaternionHermitianEJA(ConcreteEJA
, QuaternionMatrixEJA
):
2076 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
2077 matrices, the usual symmetric Jordan product, and the
2078 real-part-of-trace inner product. It has dimension `2n^2 - n` over
2083 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2087 In theory, our "field" can be any subfield of the reals::
2089 sage: QuaternionHermitianEJA(2, field=RDF)
2090 Euclidean Jordan algebra of dimension 6 over Real Double Field
2091 sage: QuaternionHermitianEJA(2, field=RR)
2092 Euclidean Jordan algebra of dimension 6 over Real Field with
2093 53 bits of precision
2097 The dimension of this algebra is `2*n^2 - n`::
2099 sage: set_random_seed()
2100 sage: n_max = QuaternionHermitianEJA._max_random_instance_size()
2101 sage: n = ZZ.random_element(1, n_max)
2102 sage: J = QuaternionHermitianEJA(n)
2103 sage: J.dimension() == 2*(n^2) - n
2106 The Jordan multiplication is what we think it is::
2108 sage: set_random_seed()
2109 sage: J = QuaternionHermitianEJA.random_instance()
2110 sage: x,y = J.random_elements(2)
2111 sage: actual = (x*y).to_matrix()
2112 sage: X = x.to_matrix()
2113 sage: Y = y.to_matrix()
2114 sage: expected = (X*Y + Y*X)/2
2115 sage: actual == expected
2117 sage: J(expected) == x*y
2120 We can change the generator prefix::
2122 sage: QuaternionHermitianEJA(2, prefix='a').gens()
2123 (a0, a1, a2, a3, a4, a5)
2125 We can construct the (trivial) algebra of rank zero::
2127 sage: QuaternionHermitianEJA(0)
2128 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2132 def _denormalized_basis(cls
, n
):
2134 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
2136 Why do we embed these? Basically, because all of numerical
2137 linear algebra assumes that you're working with vectors consisting
2138 of `n` entries from a field and scalars from the same field. There's
2139 no way to tell SageMath that (for example) the vectors contain
2140 complex numbers, while the scalar field is real.
2144 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2148 sage: set_random_seed()
2149 sage: n = ZZ.random_element(1,5)
2150 sage: B = QuaternionHermitianEJA._denormalized_basis(n)
2151 sage: all( M.is_symmetric() for M in B )
2156 Q
= QuaternionAlgebra(QQ
,-1,-1)
2159 # This is like the symmetric case, but we need to be careful:
2161 # * We want conjugate-symmetry, not just symmetry.
2162 # * The diagonal will (as a result) be real.
2166 for j
in range(i
+1):
2167 Eij
= matrix(Q
, n
, lambda k
,l
: k
==i
and l
==j
)
2169 Sij
= cls
.real_embed(Eij
)
2172 # The second, third, and fourth ones have a minus
2173 # because they're conjugated.
2174 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
2176 Sij_I
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
2178 Sij_J
= cls
.real_embed(J
*Eij
- J
*Eij
.transpose())
2180 Sij_K
= cls
.real_embed(K
*Eij
- K
*Eij
.transpose())
2183 # Since we embedded these, we can drop back to the "field" that we
2184 # started with instead of the quaternion algebra "Q".
2185 return tuple( s
.change_ring(field
) for s
in S
)
2188 def __init__(self
, n
, **kwargs
):
2189 # We know this is a valid EJA, but will double-check
2190 # if the user passes check_axioms=True.
2191 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2193 super(QuaternionHermitianEJA
, self
).__init
__(self
._denormalized
_basis
(n
),
2194 self
.jordan_product
,
2195 self
.trace_inner_product
,
2197 # TODO: this could be factored out somehow, but is left here
2198 # because the MatrixEJA is not presently a subclass of the
2199 # FDEJA class that defines rank() and one().
2200 self
.rank
.set_cache(n
)
2201 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
2202 self
.one
.set_cache(self(idV
))
2206 def _max_random_instance_size():
2208 The maximum rank of a random QuaternionHermitianEJA.
2210 return 2 # Dimension 6
2213 def random_instance(cls
, **kwargs
):
2215 Return a random instance of this type of algebra.
2217 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2218 return cls(n
, **kwargs
)
2221 class HadamardEJA(ConcreteEJA
):
2223 Return the Euclidean Jordan Algebra corresponding to the set
2224 `R^n` under the Hadamard product.
2226 Note: this is nothing more than the Cartesian product of ``n``
2227 copies of the spin algebra. Once Cartesian product algebras
2228 are implemented, this can go.
2232 sage: from mjo.eja.eja_algebra import HadamardEJA
2236 This multiplication table can be verified by hand::
2238 sage: J = HadamardEJA(3)
2239 sage: e0,e1,e2 = J.gens()
2255 We can change the generator prefix::
2257 sage: HadamardEJA(3, prefix='r').gens()
2261 def __init__(self
, n
, **kwargs
):
2263 jordan_product
= lambda x
,y
: x
2264 inner_product
= lambda x
,y
: x
2266 def jordan_product(x
,y
):
2268 return P( xi
*yi
for (xi
,yi
) in zip(x
,y
) )
2270 def inner_product(x
,y
):
2273 # New defaults for keyword arguments. Don't orthonormalize
2274 # because our basis is already orthonormal with respect to our
2275 # inner-product. Don't check the axioms, because we know this
2276 # is a valid EJA... but do double-check if the user passes
2277 # check_axioms=True. Note: we DON'T override the "check_field"
2278 # default here, because the user can pass in a field!
2279 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2280 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2282 column_basis
= tuple( b
.column() for b
in FreeModule(ZZ
, n
).basis() )
2283 super().__init
__(column_basis
, jordan_product
, inner_product
, **kwargs
)
2284 self
.rank
.set_cache(n
)
2287 self
.one
.set_cache( self
.zero() )
2289 self
.one
.set_cache( sum(self
.gens()) )
2292 def _max_random_instance_size():
2294 The maximum dimension of a random HadamardEJA.
2299 def random_instance(cls
, **kwargs
):
2301 Return a random instance of this type of algebra.
2303 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2304 return cls(n
, **kwargs
)
2307 class BilinearFormEJA(ConcreteEJA
):
2309 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2310 with the half-trace inner product and jordan product ``x*y =
2311 (<Bx,y>,y_bar>, x0*y_bar + y0*x_bar)`` where `B = 1 \times B22` is
2312 a symmetric positive-definite "bilinear form" matrix. Its
2313 dimension is the size of `B`, and it has rank two in dimensions
2314 larger than two. It reduces to the ``JordanSpinEJA`` when `B` is
2315 the identity matrix of order ``n``.
2317 We insist that the one-by-one upper-left identity block of `B` be
2318 passed in as well so that we can be passed a matrix of size zero
2319 to construct a trivial algebra.
2323 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
2324 ....: JordanSpinEJA)
2328 When no bilinear form is specified, the identity matrix is used,
2329 and the resulting algebra is the Jordan spin algebra::
2331 sage: B = matrix.identity(AA,3)
2332 sage: J0 = BilinearFormEJA(B)
2333 sage: J1 = JordanSpinEJA(3)
2334 sage: J0.multiplication_table() == J0.multiplication_table()
2337 An error is raised if the matrix `B` does not correspond to a
2338 positive-definite bilinear form::
2340 sage: B = matrix.random(QQ,2,3)
2341 sage: J = BilinearFormEJA(B)
2342 Traceback (most recent call last):
2344 ValueError: bilinear form is not positive-definite
2345 sage: B = matrix.zero(QQ,3)
2346 sage: J = BilinearFormEJA(B)
2347 Traceback (most recent call last):
2349 ValueError: bilinear form is not positive-definite
2353 We can create a zero-dimensional algebra::
2355 sage: B = matrix.identity(AA,0)
2356 sage: J = BilinearFormEJA(B)
2360 We can check the multiplication condition given in the Jordan, von
2361 Neumann, and Wigner paper (and also discussed on my "On the
2362 symmetry..." paper). Note that this relies heavily on the standard
2363 choice of basis, as does anything utilizing the bilinear form
2364 matrix. We opt not to orthonormalize the basis, because if we
2365 did, we would have to normalize the `s_{i}` in a similar manner::
2367 sage: set_random_seed()
2368 sage: n = ZZ.random_element(5)
2369 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2370 sage: B11 = matrix.identity(QQ,1)
2371 sage: B22 = M.transpose()*M
2372 sage: B = block_matrix(2,2,[ [B11,0 ],
2374 sage: J = BilinearFormEJA(B, orthonormalize=False)
2375 sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
2376 sage: V = J.vector_space()
2377 sage: sis = [ J( V([0] + (M.inverse()*ei).list()).column() )
2378 ....: for ei in eis ]
2379 sage: actual = [ sis[i]*sis[j]
2380 ....: for i in range(n-1)
2381 ....: for j in range(n-1) ]
2382 sage: expected = [ J.one() if i == j else J.zero()
2383 ....: for i in range(n-1)
2384 ....: for j in range(n-1) ]
2385 sage: actual == expected
2389 def __init__(self
, B
, **kwargs
):
2390 # The matrix "B" is supplied by the user in most cases,
2391 # so it makes sense to check whether or not its positive-
2392 # definite unless we are specifically asked not to...
2393 if ("check_axioms" not in kwargs
) or kwargs
["check_axioms"]:
2394 if not B
.is_positive_definite():
2395 raise ValueError("bilinear form is not positive-definite")
2397 # However, all of the other data for this EJA is computed
2398 # by us in manner that guarantees the axioms are
2399 # satisfied. So, again, unless we are specifically asked to
2400 # verify things, we'll skip the rest of the checks.
2401 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2403 def inner_product(x
,y
):
2404 return (y
.T
*B
*x
)[0,0]
2406 def jordan_product(x
,y
):
2412 z0
= inner_product(y
,x
)
2413 zbar
= y0
*xbar
+ x0
*ybar
2414 return P([z0
] + zbar
.list())
2417 column_basis
= tuple( b
.column() for b
in FreeModule(ZZ
, n
).basis() )
2418 super(BilinearFormEJA
, self
).__init
__(column_basis
,
2423 # The rank of this algebra is two, unless we're in a
2424 # one-dimensional ambient space (because the rank is bounded
2425 # by the ambient dimension).
2426 self
.rank
.set_cache(min(n
,2))
2429 self
.one
.set_cache( self
.zero() )
2431 self
.one
.set_cache( self
.monomial(0) )
2434 def _max_random_instance_size():
2436 The maximum dimension of a random BilinearFormEJA.
2441 def random_instance(cls
, **kwargs
):
2443 Return a random instance of this algebra.
2445 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2447 B
= matrix
.identity(ZZ
, n
)
2448 return cls(B
, **kwargs
)
2450 B11
= matrix
.identity(ZZ
, 1)
2451 M
= matrix
.random(ZZ
, n
-1)
2452 I
= matrix
.identity(ZZ
, n
-1)
2454 while alpha
.is_zero():
2455 alpha
= ZZ
.random_element().abs()
2456 B22
= M
.transpose()*M
+ alpha
*I
2458 from sage
.matrix
.special
import block_matrix
2459 B
= block_matrix(2,2, [ [B11
, ZZ(0) ],
2462 return cls(B
, **kwargs
)
2465 class JordanSpinEJA(BilinearFormEJA
):
2467 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2468 with the usual inner product and jordan product ``x*y =
2469 (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2474 sage: from mjo.eja.eja_algebra import JordanSpinEJA
2478 This multiplication table can be verified by hand::
2480 sage: J = JordanSpinEJA(4)
2481 sage: e0,e1,e2,e3 = J.gens()
2497 We can change the generator prefix::
2499 sage: JordanSpinEJA(2, prefix='B').gens()
2504 Ensure that we have the usual inner product on `R^n`::
2506 sage: set_random_seed()
2507 sage: J = JordanSpinEJA.random_instance()
2508 sage: x,y = J.random_elements(2)
2509 sage: actual = x.inner_product(y)
2510 sage: expected = x.to_vector().inner_product(y.to_vector())
2511 sage: actual == expected
2515 def __init__(self
, n
, **kwargs
):
2516 # This is a special case of the BilinearFormEJA with the
2517 # identity matrix as its bilinear form.
2518 B
= matrix
.identity(ZZ
, n
)
2520 # Don't orthonormalize because our basis is already
2521 # orthonormal with respect to our inner-product.
2522 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2524 # But also don't pass check_field=False here, because the user
2525 # can pass in a field!
2526 super(JordanSpinEJA
, self
).__init
__(B
, **kwargs
)
2529 def _max_random_instance_size():
2531 The maximum dimension of a random JordanSpinEJA.
2536 def random_instance(cls
, **kwargs
):
2538 Return a random instance of this type of algebra.
2540 Needed here to override the implementation for ``BilinearFormEJA``.
2542 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2543 return cls(n
, **kwargs
)
2546 class TrivialEJA(ConcreteEJA
):
2548 The trivial Euclidean Jordan algebra consisting of only a zero element.
2552 sage: from mjo.eja.eja_algebra import TrivialEJA
2556 sage: J = TrivialEJA()
2563 sage: 7*J.one()*12*J.one()
2565 sage: J.one().inner_product(J.one())
2567 sage: J.one().norm()
2569 sage: J.one().subalgebra_generated_by()
2570 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2575 def __init__(self
, **kwargs
):
2576 jordan_product
= lambda x
,y
: x
2577 inner_product
= lambda x
,y
: 0
2580 # New defaults for keyword arguments
2581 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2582 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2584 super(TrivialEJA
, self
).__init
__(basis
,
2588 # The rank is zero using my definition, namely the dimension of the
2589 # largest subalgebra generated by any element.
2590 self
.rank
.set_cache(0)
2591 self
.one
.set_cache( self
.zero() )
2594 def random_instance(cls
, **kwargs
):
2595 # We don't take a "size" argument so the superclass method is
2596 # inappropriate for us.
2597 return cls(**kwargs
)
2599 # class DirectSumEJA(ConcreteEJA):
2601 # The external (orthogonal) direct sum of two other Euclidean Jordan
2602 # algebras. Essentially the Cartesian product of its two factors.
2603 # Every Euclidean Jordan algebra decomposes into an orthogonal
2604 # direct sum of simple Euclidean Jordan algebras, so no generality
2605 # is lost by providing only this construction.
2609 # sage: from mjo.eja.eja_algebra import (random_eja,
2610 # ....: HadamardEJA,
2611 # ....: RealSymmetricEJA,
2612 # ....: DirectSumEJA)
2616 # sage: J1 = HadamardEJA(2)
2617 # sage: J2 = RealSymmetricEJA(3)
2618 # sage: J = DirectSumEJA(J1,J2)
2619 # sage: J.dimension()
2626 # The external direct sum construction is only valid when the two factors
2627 # have the same base ring; an error is raised otherwise::
2629 # sage: set_random_seed()
2630 # sage: J1 = random_eja(field=AA)
2631 # sage: J2 = random_eja(field=QQ,orthonormalize=False)
2632 # sage: J = DirectSumEJA(J1,J2)
2633 # Traceback (most recent call last):
2635 # ValueError: algebras must share the same base field
2638 # def __init__(self, J1, J2, **kwargs):
2639 # if J1.base_ring() != J2.base_ring():
2640 # raise ValueError("algebras must share the same base field")
2641 # field = J1.base_ring()
2643 # self._factors = (J1, J2)
2644 # n1 = J1.dimension()
2645 # n2 = J2.dimension()
2647 # V = VectorSpace(field, n)
2648 # mult_table = [ [ V.zero() for j in range(i+1) ]
2649 # for i in range(n) ]
2650 # for i in range(n1):
2651 # for j in range(i+1):
2652 # p = (J1.monomial(i)*J1.monomial(j)).to_vector()
2653 # mult_table[i][j] = V(p.list() + [field.zero()]*n2)
2655 # for i in range(n2):
2656 # for j in range(i+1):
2657 # p = (J2.monomial(i)*J2.monomial(j)).to_vector()
2658 # mult_table[n1+i][n1+j] = V([field.zero()]*n1 + p.list())
2660 # # TODO: build the IP table here from the two constituent IP
2661 # # matrices (it'll be block diagonal, I think).
2662 # ip_table = [ [ field.zero() for j in range(i+1) ]
2663 # for i in range(n) ]
2664 # super(DirectSumEJA, self).__init__(field,
2667 # check_axioms=False,
2669 # self.rank.set_cache(J1.rank() + J2.rank())
2672 # def factors(self):
2674 # Return the pair of this algebra's factors.
2678 # sage: from mjo.eja.eja_algebra import (HadamardEJA,
2679 # ....: JordanSpinEJA,
2680 # ....: DirectSumEJA)
2684 # sage: J1 = HadamardEJA(2, field=QQ)
2685 # sage: J2 = JordanSpinEJA(3, field=QQ)
2686 # sage: J = DirectSumEJA(J1,J2)
2688 # (Euclidean Jordan algebra of dimension 2 over Rational Field,
2689 # Euclidean Jordan algebra of dimension 3 over Rational Field)
2692 # return self._factors
2694 # def projections(self):
2696 # Return a pair of projections onto this algebra's factors.
2700 # sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
2701 # ....: ComplexHermitianEJA,
2702 # ....: DirectSumEJA)
2706 # sage: J1 = JordanSpinEJA(2)
2707 # sage: J2 = ComplexHermitianEJA(2)
2708 # sage: J = DirectSumEJA(J1,J2)
2709 # sage: (pi_left, pi_right) = J.projections()
2710 # sage: J.one().to_vector()
2711 # (1, 0, 1, 0, 0, 1)
2712 # sage: pi_left(J.one()).to_vector()
2714 # sage: pi_right(J.one()).to_vector()
2718 # (J1,J2) = self.factors()
2719 # m = J1.dimension()
2720 # n = J2.dimension()
2721 # V_basis = self.vector_space().basis()
2722 # # Need to specify the dimensions explicitly so that we don't
2723 # # wind up with a zero-by-zero matrix when we want e.g. a
2724 # # zero-by-two matrix (important for composing things).
2725 # P1 = matrix(self.base_ring(), m, m+n, V_basis[:m])
2726 # P2 = matrix(self.base_ring(), n, m+n, V_basis[m:])
2727 # pi_left = FiniteDimensionalEJAOperator(self,J1,P1)
2728 # pi_right = FiniteDimensionalEJAOperator(self,J2,P2)
2729 # return (pi_left, pi_right)
2731 # def inclusions(self):
2733 # Return the pair of inclusion maps from our factors into us.
2737 # sage: from mjo.eja.eja_algebra import (random_eja,
2738 # ....: JordanSpinEJA,
2739 # ....: RealSymmetricEJA,
2740 # ....: DirectSumEJA)
2744 # sage: J1 = JordanSpinEJA(3)
2745 # sage: J2 = RealSymmetricEJA(2)
2746 # sage: J = DirectSumEJA(J1,J2)
2747 # sage: (iota_left, iota_right) = J.inclusions()
2748 # sage: iota_left(J1.zero()) == J.zero()
2750 # sage: iota_right(J2.zero()) == J.zero()
2752 # sage: J1.one().to_vector()
2754 # sage: iota_left(J1.one()).to_vector()
2755 # (1, 0, 0, 0, 0, 0)
2756 # sage: J2.one().to_vector()
2758 # sage: iota_right(J2.one()).to_vector()
2759 # (0, 0, 0, 1, 0, 1)
2760 # sage: J.one().to_vector()
2761 # (1, 0, 0, 1, 0, 1)
2765 # Composing a projection with the corresponding inclusion should
2766 # produce the identity map, and mismatching them should produce
2769 # sage: set_random_seed()
2770 # sage: J1 = random_eja()
2771 # sage: J2 = random_eja()
2772 # sage: J = DirectSumEJA(J1,J2)
2773 # sage: (iota_left, iota_right) = J.inclusions()
2774 # sage: (pi_left, pi_right) = J.projections()
2775 # sage: pi_left*iota_left == J1.one().operator()
2777 # sage: pi_right*iota_right == J2.one().operator()
2779 # sage: (pi_left*iota_right).is_zero()
2781 # sage: (pi_right*iota_left).is_zero()
2785 # (J1,J2) = self.factors()
2786 # m = J1.dimension()
2787 # n = J2.dimension()
2788 # V_basis = self.vector_space().basis()
2789 # # Need to specify the dimensions explicitly so that we don't
2790 # # wind up with a zero-by-zero matrix when we want e.g. a
2791 # # two-by-zero matrix (important for composing things).
2792 # I1 = matrix.column(self.base_ring(), m, m+n, V_basis[:m])
2793 # I2 = matrix.column(self.base_ring(), n, m+n, V_basis[m:])
2794 # iota_left = FiniteDimensionalEJAOperator(J1,self,I1)
2795 # iota_right = FiniteDimensionalEJAOperator(J2,self,I2)
2796 # return (iota_left, iota_right)
2798 # def inner_product(self, x, y):
2800 # The standard Cartesian inner-product.
2802 # We project ``x`` and ``y`` onto our factors, and add up the
2803 # inner-products from the subalgebras.
2808 # sage: from mjo.eja.eja_algebra import (HadamardEJA,
2809 # ....: QuaternionHermitianEJA,
2810 # ....: DirectSumEJA)
2814 # sage: J1 = HadamardEJA(3,field=QQ)
2815 # sage: J2 = QuaternionHermitianEJA(2,field=QQ,orthonormalize=False)
2816 # sage: J = DirectSumEJA(J1,J2)
2817 # sage: x1 = J1.one()
2819 # sage: y1 = J2.one()
2821 # sage: x1.inner_product(x2)
2823 # sage: y1.inner_product(y2)
2825 # sage: J.one().inner_product(J.one())
2829 # (pi_left, pi_right) = self.projections()
2835 # return (x1.inner_product(y1) + x2.inner_product(y2))
2839 random_eja
= ConcreteEJA
.random_instance