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 FiniteDimensionalEuclideanJordanAlgebraElement
33 from mjo
.eja
.eja_operator
import FiniteDimensionalEuclideanJordanAlgebraOperator
34 from mjo
.eja
.eja_utils
import _mat2vec
36 class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule
):
38 def _coerce_map_from_base_ring(self
):
40 Disable the map from the base ring into the algebra.
42 Performing a nonsense conversion like this automatically
43 is counterpedagogical. The fallback is to try the usual
44 element constructor, which should also fail.
48 sage: from mjo.eja.eja_algebra import random_eja
52 sage: set_random_seed()
53 sage: J = random_eja()
55 Traceback (most recent call last):
57 ValueError: not an element of this algebra
73 sage: from mjo.eja.eja_algebra import (
74 ....: FiniteDimensionalEuclideanJordanAlgebra,
80 By definition, Jordan multiplication commutes::
82 sage: set_random_seed()
83 sage: J = random_eja()
84 sage: x,y = J.random_elements(2)
90 The ``field`` we're given must be real with ``check_field=True``::
92 sage: JordanSpinEJA(2,QQbar)
93 Traceback (most recent call last):
95 ValueError: scalar field is not real
97 The multiplication table must be square with ``check_axioms=True``::
99 sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,((),()))
100 Traceback (most recent call last):
102 ValueError: multiplication table is not square
106 if not field
.is_subring(RR
):
107 # Note: this does return true for the real algebraic
108 # field, the rationals, and any quadratic field where
109 # we've specified a real embedding.
110 raise ValueError("scalar field is not real")
112 # The multiplication table had better be square
115 if not all( len(l
) == n
for l
in mult_table
):
116 raise ValueError("multiplication table is not square")
118 self
._matrix
_basis
= matrix_basis
121 category
= MagmaticAlgebras(field
).FiniteDimensional()
122 category
= category
.WithBasis().Unital()
124 fda
= super(FiniteDimensionalEuclideanJordanAlgebra
, self
)
129 self
.print_options(bracket
='')
131 # The multiplication table we're given is necessarily in terms
132 # of vectors, because we don't have an algebra yet for
133 # anything to be an element of. However, it's faster in the
134 # long run to have the multiplication table be in terms of
135 # algebra elements. We do this after calling the superclass
136 # constructor so that from_vector() knows what to do.
137 self
._multiplication
_table
= [
138 list(map(lambda x
: self
.from_vector(x
), ls
))
143 if not self
._is
_commutative
():
144 raise ValueError("algebra is not commutative")
145 if not self
._is
_jordanian
():
146 raise ValueError("Jordan identity does not hold")
147 if not self
._inner
_product
_is
_associative
():
148 raise ValueError("inner product is not associative")
150 def _element_constructor_(self
, elt
):
152 Construct an element of this algebra from its vector or matrix
155 This gets called only after the parent element _call_ method
156 fails to find a coercion for the argument.
160 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
162 ....: RealSymmetricEJA)
166 The identity in `S^n` is converted to the identity in the EJA::
168 sage: J = RealSymmetricEJA(3)
169 sage: I = matrix.identity(QQ,3)
170 sage: J(I) == J.one()
173 This skew-symmetric matrix can't be represented in the EJA::
175 sage: J = RealSymmetricEJA(3)
176 sage: A = matrix(QQ,3, lambda i,j: i-j)
178 Traceback (most recent call last):
180 ValueError: not an element of this algebra
184 Ensure that we can convert any element of the two non-matrix
185 simple algebras (whose matrix representations are columns)
186 back and forth faithfully::
188 sage: set_random_seed()
189 sage: J = HadamardEJA.random_instance()
190 sage: x = J.random_element()
191 sage: J(x.to_vector().column()) == x
193 sage: J = JordanSpinEJA.random_instance()
194 sage: x = J.random_element()
195 sage: J(x.to_vector().column()) == x
198 msg
= "not an element of this algebra"
200 # The superclass implementation of random_element()
201 # needs to be able to coerce "0" into the algebra.
203 elif elt
in self
.base_ring():
204 # Ensure that no base ring -> algebra coercion is performed
205 # by this method. There's some stupidity in sage that would
206 # otherwise propagate to this method; for example, sage thinks
207 # that the integer 3 belongs to the space of 2-by-2 matrices.
208 raise ValueError(msg
)
210 if elt
not in self
.matrix_space():
211 raise ValueError(msg
)
213 # Thanks for nothing! Matrix spaces aren't vector spaces in
214 # Sage, so we have to figure out its matrix-basis coordinates
215 # ourselves. We use the basis space's ring instead of the
216 # element's ring because the basis space might be an algebraic
217 # closure whereas the base ring of the 3-by-3 identity matrix
218 # could be QQ instead of QQbar.
219 V
= VectorSpace(self
.base_ring(), elt
.nrows()*elt
.ncols())
220 W
= V
.span_of_basis( _mat2vec(s
) for s
in self
.matrix_basis() )
223 coords
= W
.coordinate_vector(_mat2vec(elt
))
224 except ArithmeticError: # vector is not in free module
225 raise ValueError(msg
)
227 return self
.from_vector(coords
)
231 Return a string representation of ``self``.
235 sage: from mjo.eja.eja_algebra import JordanSpinEJA
239 Ensure that it says what we think it says::
241 sage: JordanSpinEJA(2, field=AA)
242 Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
243 sage: JordanSpinEJA(3, field=RDF)
244 Euclidean Jordan algebra of dimension 3 over Real Double Field
247 fmt
= "Euclidean Jordan algebra of dimension {} over {}"
248 return fmt
.format(self
.dimension(), self
.base_ring())
250 def product_on_basis(self
, i
, j
):
251 return self
._multiplication
_table
[i
][j
]
253 def _is_commutative(self
):
255 Whether or not this algebra's multiplication table is commutative.
257 This method should of course always return ``True``, unless
258 this algebra was constructed with ``check_axioms=False`` and
259 passed an invalid multiplication table.
261 return all( self
.product_on_basis(i
,j
) == self
.product_on_basis(i
,j
)
262 for i
in range(self
.dimension())
263 for j
in range(self
.dimension()) )
265 def _is_jordanian(self
):
267 Whether or not this algebra's multiplication table respects the
268 Jordan identity `(x^{2})(xy) = x(x^{2}y)`.
270 We only check one arrangement of `x` and `y`, so for a
271 ``True`` result to be truly true, you should also check
272 :meth:`_is_commutative`. This method should of course always
273 return ``True``, unless this algebra was constructed with
274 ``check_axioms=False`` and passed an invalid multiplication table.
276 return all( (self
.monomial(i
)**2)*(self
.monomial(i
)*self
.monomial(j
))
278 (self
.monomial(i
))*((self
.monomial(i
)**2)*self
.monomial(j
))
279 for i
in range(self
.dimension())
280 for j
in range(self
.dimension()) )
282 def _inner_product_is_associative(self
):
284 Return whether or not this algebra's inner product `B` is
285 associative; that is, whether or not `B(xy,z) = B(x,yz)`.
287 This method should of course always return ``True``, unless
288 this algebra was constructed with ``check_axioms=False`` and
289 passed an invalid multiplication table.
292 # Used to check whether or not something is zero in an inexact
293 # ring. This number is sufficient to allow the construction of
294 # QuaternionHermitianEJA(2, RDF) with check_axioms=True.
297 for i
in range(self
.dimension()):
298 for j
in range(self
.dimension()):
299 for k
in range(self
.dimension()):
303 diff
= (x
*y
).inner_product(z
) - x
.inner_product(y
*z
)
305 if self
.base_ring().is_exact():
309 if diff
.abs() > epsilon
:
315 def characteristic_polynomial_of(self
):
317 Return the algebra's "characteristic polynomial of" function,
318 which is itself a multivariate polynomial that, when evaluated
319 at the coordinates of some algebra element, returns that
320 element's characteristic polynomial.
322 The resulting polynomial has `n+1` variables, where `n` is the
323 dimension of this algebra. The first `n` variables correspond to
324 the coordinates of an algebra element: when evaluated at the
325 coordinates of an algebra element with respect to a certain
326 basis, the result is a univariate polynomial (in the one
327 remaining variable ``t``), namely the characteristic polynomial
332 sage: from mjo.eja.eja_algebra import JordanSpinEJA, TrivialEJA
336 The characteristic polynomial in the spin algebra is given in
337 Alizadeh, Example 11.11::
339 sage: J = JordanSpinEJA(3)
340 sage: p = J.characteristic_polynomial_of(); p
341 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
342 sage: xvec = J.one().to_vector()
346 By definition, the characteristic polynomial is a monic
347 degree-zero polynomial in a rank-zero algebra. Note that
348 Cayley-Hamilton is indeed satisfied since the polynomial
349 ``1`` evaluates to the identity element of the algebra on
352 sage: J = TrivialEJA()
353 sage: J.characteristic_polynomial_of()
360 # The list of coefficient polynomials a_0, a_1, a_2, ..., a_(r-1).
361 a
= self
._charpoly
_coefficients
()
363 # We go to a bit of trouble here to reorder the
364 # indeterminates, so that it's easier to evaluate the
365 # characteristic polynomial at x's coordinates and get back
366 # something in terms of t, which is what we want.
367 S
= PolynomialRing(self
.base_ring(),'t')
371 S
= PolynomialRing(S
, R
.variable_names())
374 return (t
**r
+ sum( a
[k
]*(t
**k
) for k
in range(r
) ))
376 def coordinate_polynomial_ring(self
):
378 The multivariate polynomial ring in which this algebra's
379 :meth:`characteristic_polynomial_of` lives.
383 sage: from mjo.eja.eja_algebra import (HadamardEJA,
384 ....: RealSymmetricEJA)
388 sage: J = HadamardEJA(2)
389 sage: J.coordinate_polynomial_ring()
390 Multivariate Polynomial Ring in X1, X2...
391 sage: J = RealSymmetricEJA(3,QQ)
392 sage: J.coordinate_polynomial_ring()
393 Multivariate Polynomial Ring in X1, X2, X3, X4, X5, X6...
396 var_names
= tuple( "X%d" % z
for z
in range(1, self
.dimension()+1) )
397 return PolynomialRing(self
.base_ring(), var_names
)
399 def inner_product(self
, x
, y
):
401 The inner product associated with this Euclidean Jordan algebra.
403 Defaults to the trace inner product, but can be overridden by
404 subclasses if they are sure that the necessary properties are
409 sage: from mjo.eja.eja_algebra import (random_eja,
411 ....: BilinearFormEJA)
415 Our inner product is "associative," which means the following for
416 a symmetric bilinear form::
418 sage: set_random_seed()
419 sage: J = random_eja()
420 sage: x,y,z = J.random_elements(3)
421 sage: (x*y).inner_product(z) == y.inner_product(x*z)
426 Ensure that this is the usual inner product for the algebras
429 sage: set_random_seed()
430 sage: J = HadamardEJA.random_instance()
431 sage: x,y = J.random_elements(2)
432 sage: actual = x.inner_product(y)
433 sage: expected = x.to_vector().inner_product(y.to_vector())
434 sage: actual == expected
437 Ensure that this is one-half of the trace inner-product in a
438 BilinearFormEJA that isn't just the reals (when ``n`` isn't
439 one). This is in Faraut and Koranyi, and also my "On the
442 sage: set_random_seed()
443 sage: J = BilinearFormEJA.random_instance()
444 sage: n = J.dimension()
445 sage: x = J.random_element()
446 sage: y = J.random_element()
447 sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
450 B
= self
._inner
_product
_matrix
451 return (B
*x
.to_vector()).inner_product(y
.to_vector())
454 def is_trivial(self
):
456 Return whether or not this algebra is trivial.
458 A trivial algebra contains only the zero element.
462 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
467 sage: J = ComplexHermitianEJA(3)
473 sage: J = TrivialEJA()
478 return self
.dimension() == 0
481 def multiplication_table(self
):
483 Return a visual representation of this algebra's multiplication
484 table (on basis elements).
488 sage: from mjo.eja.eja_algebra import JordanSpinEJA
492 sage: J = JordanSpinEJA(4)
493 sage: J.multiplication_table()
494 +----++----+----+----+----+
495 | * || e0 | e1 | e2 | e3 |
496 +====++====+====+====+====+
497 | e0 || e0 | e1 | e2 | e3 |
498 +----++----+----+----+----+
499 | e1 || e1 | e0 | 0 | 0 |
500 +----++----+----+----+----+
501 | e2 || e2 | 0 | e0 | 0 |
502 +----++----+----+----+----+
503 | e3 || e3 | 0 | 0 | e0 |
504 +----++----+----+----+----+
507 M
= list(self
._multiplication
_table
) # copy
508 for i
in range(len(M
)):
509 # M had better be "square"
510 M
[i
] = [self
.monomial(i
)] + M
[i
]
511 M
= [["*"] + list(self
.gens())] + M
512 return table(M
, header_row
=True, header_column
=True, frame
=True)
515 def matrix_basis(self
):
517 Return an (often more natural) representation of this algebras
518 basis as an ordered tuple of matrices.
520 Every finite-dimensional Euclidean Jordan Algebra is a, up to
521 Jordan isomorphism, a direct sum of five simple
522 algebras---four of which comprise Hermitian matrices. And the
523 last type of algebra can of course be thought of as `n`-by-`1`
524 column matrices (ambiguusly called column vectors) to avoid
525 special cases. As a result, matrices (and column vectors) are
526 a natural representation format for Euclidean Jordan algebra
529 But, when we construct an algebra from a basis of matrices,
530 those matrix representations are lost in favor of coordinate
531 vectors *with respect to* that basis. We could eventually
532 convert back if we tried hard enough, but having the original
533 representations handy is valuable enough that we simply store
534 them and return them from this method.
536 Why implement this for non-matrix algebras? Avoiding special
537 cases for the :class:`BilinearFormEJA` pays with simplicity in
538 its own right. But mainly, we would like to be able to assume
539 that elements of a :class:`DirectSumEJA` can be displayed
540 nicely, without having to have special classes for direct sums
541 one of whose components was a matrix algebra.
545 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
546 ....: RealSymmetricEJA)
550 sage: J = RealSymmetricEJA(2)
552 Finite family {0: e0, 1: e1, 2: e2}
553 sage: J.matrix_basis()
555 [1 0] [ 0 0.7071067811865475?] [0 0]
556 [0 0], [0.7071067811865475? 0], [0 1]
561 sage: J = JordanSpinEJA(2)
563 Finite family {0: e0, 1: e1}
564 sage: J.matrix_basis()
570 if self
._matrix
_basis
is None:
571 M
= self
.matrix_space()
572 return tuple( M(b
.to_vector()) for b
in self
.basis() )
574 return self
._matrix
_basis
577 def matrix_space(self
):
579 Return the matrix space in which this algebra's elements live, if
580 we think of them as matrices (including column vectors of the
583 Generally this will be an `n`-by-`1` column-vector space,
584 except when the algebra is trivial. There it's `n`-by-`n`
585 (where `n` is zero), to ensure that two elements of the matrix
586 space (empty matrices) can be multiplied.
588 Matrix algebras override this with something more useful.
590 if self
.is_trivial():
591 return MatrixSpace(self
.base_ring(), 0)
592 elif self
._matrix
_basis
is None or len(self
._matrix
_basis
) == 0:
593 return MatrixSpace(self
.base_ring(), self
.dimension(), 1)
595 return self
._matrix
_basis
[0].matrix_space()
601 Return the unit element of this algebra.
605 sage: from mjo.eja.eja_algebra import (HadamardEJA,
610 sage: J = HadamardEJA(5)
612 e0 + e1 + e2 + e3 + e4
616 The identity element acts like the identity::
618 sage: set_random_seed()
619 sage: J = random_eja()
620 sage: x = J.random_element()
621 sage: J.one()*x == x and x*J.one() == x
624 The matrix of the unit element's operator is the identity::
626 sage: set_random_seed()
627 sage: J = random_eja()
628 sage: actual = J.one().operator().matrix()
629 sage: expected = matrix.identity(J.base_ring(), J.dimension())
630 sage: actual == expected
633 Ensure that the cached unit element (often precomputed by
634 hand) agrees with the computed one::
636 sage: set_random_seed()
637 sage: J = random_eja()
638 sage: cached = J.one()
639 sage: J.one.clear_cache()
640 sage: J.one() == cached
644 # We can brute-force compute the matrices of the operators
645 # that correspond to the basis elements of this algebra.
646 # If some linear combination of those basis elements is the
647 # algebra identity, then the same linear combination of
648 # their matrices has to be the identity matrix.
650 # Of course, matrices aren't vectors in sage, so we have to
651 # appeal to the "long vectors" isometry.
652 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
654 # Now we use basic linear algebra to find the coefficients,
655 # of the matrices-as-vectors-linear-combination, which should
656 # work for the original algebra basis too.
657 A
= matrix(self
.base_ring(), oper_vecs
)
659 # We used the isometry on the left-hand side already, but we
660 # still need to do it for the right-hand side. Recall that we
661 # wanted something that summed to the identity matrix.
662 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
664 # Now if there's an identity element in the algebra, this
665 # should work. We solve on the left to avoid having to
666 # transpose the matrix "A".
667 return self
.from_vector(A
.solve_left(b
))
670 def peirce_decomposition(self
, c
):
672 The Peirce decomposition of this algebra relative to the
675 In the future, this can be extended to a complete system of
676 orthogonal idempotents.
680 - ``c`` -- an idempotent of this algebra.
684 A triple (J0, J5, J1) containing two subalgebras and one subspace
687 - ``J0`` -- the algebra on the eigenspace of ``c.operator()``
688 corresponding to the eigenvalue zero.
690 - ``J5`` -- the eigenspace (NOT a subalgebra) of ``c.operator()``
691 corresponding to the eigenvalue one-half.
693 - ``J1`` -- the algebra on the eigenspace of ``c.operator()``
694 corresponding to the eigenvalue one.
696 These are the only possible eigenspaces for that operator, and this
697 algebra is a direct sum of them. The spaces ``J0`` and ``J1`` are
698 orthogonal, and are subalgebras of this algebra with the appropriate
703 sage: from mjo.eja.eja_algebra import random_eja, RealSymmetricEJA
707 The canonical example comes from the symmetric matrices, which
708 decompose into diagonal and off-diagonal parts::
710 sage: J = RealSymmetricEJA(3)
711 sage: C = matrix(QQ, [ [1,0,0],
715 sage: J0,J5,J1 = J.peirce_decomposition(c)
717 Euclidean Jordan algebra of dimension 1...
719 Vector space of degree 6 and dimension 2...
721 Euclidean Jordan algebra of dimension 3...
722 sage: J0.one().to_matrix()
726 sage: orig_df = AA.options.display_format
727 sage: AA.options.display_format = 'radical'
728 sage: J.from_vector(J5.basis()[0]).to_matrix()
732 sage: J.from_vector(J5.basis()[1]).to_matrix()
736 sage: AA.options.display_format = orig_df
737 sage: J1.one().to_matrix()
744 Every algebra decomposes trivially with respect to its identity
747 sage: set_random_seed()
748 sage: J = random_eja()
749 sage: J0,J5,J1 = J.peirce_decomposition(J.one())
750 sage: J0.dimension() == 0 and J5.dimension() == 0
752 sage: J1.superalgebra() == J and J1.dimension() == J.dimension()
755 The decomposition is into eigenspaces, and its components are
756 therefore necessarily orthogonal. Moreover, the identity
757 elements in the two subalgebras are the projections onto their
758 respective subspaces of the superalgebra's identity element::
760 sage: set_random_seed()
761 sage: J = random_eja()
762 sage: x = J.random_element()
763 sage: if not J.is_trivial():
764 ....: while x.is_nilpotent():
765 ....: x = J.random_element()
766 sage: c = x.subalgebra_idempotent()
767 sage: J0,J5,J1 = J.peirce_decomposition(c)
769 sage: for (w,y,z) in zip(J0.basis(), J5.basis(), J1.basis()):
770 ....: w = w.superalgebra_element()
771 ....: y = J.from_vector(y)
772 ....: z = z.superalgebra_element()
773 ....: ipsum += w.inner_product(y).abs()
774 ....: ipsum += w.inner_product(z).abs()
775 ....: ipsum += y.inner_product(z).abs()
778 sage: J1(c) == J1.one()
780 sage: J0(J.one() - c) == J0.one()
784 if not c
.is_idempotent():
785 raise ValueError("element is not idempotent: %s" % c
)
787 from mjo
.eja
.eja_subalgebra
import FiniteDimensionalEuclideanJordanSubalgebra
789 # Default these to what they should be if they turn out to be
790 # trivial, because eigenspaces_left() won't return eigenvalues
791 # corresponding to trivial spaces (e.g. it returns only the
792 # eigenspace corresponding to lambda=1 if you take the
793 # decomposition relative to the identity element).
794 trivial
= FiniteDimensionalEuclideanJordanSubalgebra(self
, ())
795 J0
= trivial
# eigenvalue zero
796 J5
= VectorSpace(self
.base_ring(), 0) # eigenvalue one-half
797 J1
= trivial
# eigenvalue one
799 for (eigval
, eigspace
) in c
.operator().matrix().right_eigenspaces():
800 if eigval
== ~
(self
.base_ring()(2)):
803 gens
= tuple( self
.from_vector(b
) for b
in eigspace
.basis() )
804 subalg
= FiniteDimensionalEuclideanJordanSubalgebra(self
,
812 raise ValueError("unexpected eigenvalue: %s" % eigval
)
817 def random_element(self
, thorough
=False):
819 Return a random element of this algebra.
821 Our algebra superclass method only returns a linear
822 combination of at most two basis elements. We instead
823 want the vector space "random element" method that
824 returns a more diverse selection.
828 - ``thorough`` -- (boolean; default False) whether or not we
829 should generate irrational coefficients for the random
830 element when our base ring is irrational; this slows the
831 algebra operations to a crawl, but any truly random method
835 # For a general base ring... maybe we can trust this to do the
836 # right thing? Unlikely, but.
837 V
= self
.vector_space()
838 v
= V
.random_element()
840 if self
.base_ring() is AA
:
841 # The "random element" method of the algebraic reals is
842 # stupid at the moment, and only returns integers between
843 # -2 and 2, inclusive:
845 # https://trac.sagemath.org/ticket/30875
847 # Instead, we implement our own "random vector" method,
848 # and then coerce that into the algebra. We use the vector
849 # space degree here instead of the dimension because a
850 # subalgebra could (for example) be spanned by only two
851 # vectors, each with five coordinates. We need to
852 # generate all five coordinates.
854 v
*= QQbar
.random_element().real()
856 v
*= QQ
.random_element()
858 return self
.from_vector(V
.coordinate_vector(v
))
860 def random_elements(self
, count
, thorough
=False):
862 Return ``count`` random elements as a tuple.
866 - ``thorough`` -- (boolean; default False) whether or not we
867 should generate irrational coefficients for the random
868 elements when our base ring is irrational; this slows the
869 algebra operations to a crawl, but any truly random method
874 sage: from mjo.eja.eja_algebra import JordanSpinEJA
878 sage: J = JordanSpinEJA(3)
879 sage: x,y,z = J.random_elements(3)
880 sage: all( [ x in J, y in J, z in J ])
882 sage: len( J.random_elements(10) ) == 10
886 return tuple( self
.random_element(thorough
)
887 for idx
in range(count
) )
891 def _charpoly_coefficients(self
):
893 The `r` polynomial coefficients of the "characteristic polynomial
897 R
= self
.coordinate_polynomial_ring()
899 F
= R
.fraction_field()
902 # From a result in my book, these are the entries of the
903 # basis representation of L_x.
904 return sum( vars[k
]*self
.monomial(k
).operator().matrix()[i
,j
]
907 L_x
= matrix(F
, n
, n
, L_x_i_j
)
910 if self
.rank
.is_in_cache():
912 # There's no need to pad the system with redundant
913 # columns if we *know* they'll be redundant.
916 # Compute an extra power in case the rank is equal to
917 # the dimension (otherwise, we would stop at x^(r-1)).
918 x_powers
= [ (L_x
**k
)*self
.one().to_vector()
919 for k
in range(n
+1) ]
920 A
= matrix
.column(F
, x_powers
[:n
])
921 AE
= A
.extended_echelon_form()
928 # The theory says that only the first "r" coefficients are
929 # nonzero, and they actually live in the original polynomial
930 # ring and not the fraction field. We negate them because
931 # in the actual characteristic polynomial, they get moved
932 # to the other side where x^r lives.
933 return -A_rref
.solve_right(E
*b
).change_ring(R
)[:r
]
938 Return the rank of this EJA.
940 This is a cached method because we know the rank a priori for
941 all of the algebras we can construct. Thus we can avoid the
942 expensive ``_charpoly_coefficients()`` call unless we truly
943 need to compute the whole characteristic polynomial.
947 sage: from mjo.eja.eja_algebra import (HadamardEJA,
949 ....: RealSymmetricEJA,
950 ....: ComplexHermitianEJA,
951 ....: QuaternionHermitianEJA,
956 The rank of the Jordan spin algebra is always two::
958 sage: JordanSpinEJA(2).rank()
960 sage: JordanSpinEJA(3).rank()
962 sage: JordanSpinEJA(4).rank()
965 The rank of the `n`-by-`n` Hermitian real, complex, or
966 quaternion matrices is `n`::
968 sage: RealSymmetricEJA(4).rank()
970 sage: ComplexHermitianEJA(3).rank()
972 sage: QuaternionHermitianEJA(2).rank()
977 Ensure that every EJA that we know how to construct has a
978 positive integer rank, unless the algebra is trivial in
979 which case its rank will be zero::
981 sage: set_random_seed()
982 sage: J = random_eja()
986 sage: r > 0 or (r == 0 and J.is_trivial())
989 Ensure that computing the rank actually works, since the ranks
990 of all simple algebras are known and will be cached by default::
992 sage: J = HadamardEJA(4)
993 sage: J.rank.clear_cache()
999 sage: J = JordanSpinEJA(4)
1000 sage: J.rank.clear_cache()
1006 sage: J = RealSymmetricEJA(3)
1007 sage: J.rank.clear_cache()
1013 sage: J = ComplexHermitianEJA(2)
1014 sage: J.rank.clear_cache()
1020 sage: J = QuaternionHermitianEJA(2)
1021 sage: J.rank.clear_cache()
1025 return len(self
._charpoly
_coefficients
())
1028 def vector_space(self
):
1030 Return the vector space that underlies this algebra.
1034 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1038 sage: J = RealSymmetricEJA(2)
1039 sage: J.vector_space()
1040 Vector space of dimension 3 over...
1043 return self
.zero().to_vector().parent().ambient_vector_space()
1046 Element
= FiniteDimensionalEuclideanJordanAlgebraElement
1049 class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra
):
1051 Algebras whose basis consists of vectors with rational
1052 entries. Equivalently, algebras whose multiplication tables
1053 contain only rational coefficients.
1055 When an EJA has a basis that can be made rational, we can speed up
1056 the computation of its characteristic polynomial by doing it over
1057 ``QQ``. All of the named EJA constructors that we provide fall
1061 def _charpoly_coefficients(self
):
1063 Override the parent method with something that tries to compute
1064 over a faster (non-extension) field.
1068 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1072 The base ring of the resulting polynomial coefficients is what
1073 it should be, and not the rationals (unless the algebra was
1074 already over the rationals)::
1076 sage: J = JordanSpinEJA(3)
1077 sage: J._charpoly_coefficients()
1078 (X1^2 - X2^2 - X3^2, -2*X1)
1079 sage: a0 = J._charpoly_coefficients()[0]
1081 Algebraic Real Field
1082 sage: a0.base_ring()
1083 Algebraic Real Field
1086 if self
.base_ring() is QQ
:
1087 # There's no need to construct *another* algebra over the
1088 # rationals if this one is already over the rationals.
1089 superclass
= super(RationalBasisEuclideanJordanAlgebra
, self
)
1090 return superclass
._charpoly
_coefficients
()
1093 map(lambda x
: x
.to_vector(), ls
)
1094 for ls
in self
._multiplication
_table
1097 # Do the computation over the rationals. The answer will be
1098 # the same, because our basis coordinates are (essentially)
1100 J
= FiniteDimensionalEuclideanJordanAlgebra(QQ
,
1104 a
= J
._charpoly
_coefficients
()
1105 return tuple(map(lambda x
: x
.change_ring(self
.base_ring()), a
))
1108 class ConcreteEuclideanJordanAlgebra
:
1110 A class for the Euclidean Jordan algebras that we know by name.
1112 These are the Jordan algebras whose basis, multiplication table,
1113 rank, and so on are known a priori. More to the point, they are
1114 the Euclidean Jordan algebras for which we are able to conjure up
1115 a "random instance."
1119 sage: from mjo.eja.eja_algebra import ConcreteEuclideanJordanAlgebra
1123 Our basis is normalized with respect to the algebra's inner
1124 product, unless we specify otherwise::
1126 sage: set_random_seed()
1127 sage: J = ConcreteEuclideanJordanAlgebra.random_instance()
1128 sage: all( b.norm() == 1 for b in J.gens() )
1131 Since our basis is orthonormal with respect to the algebra's inner
1132 product, and since we know that this algebra is an EJA, any
1133 left-multiplication operator's matrix will be symmetric because
1134 natural->EJA basis representation is an isometry and within the
1135 EJA the operator is self-adjoint by the Jordan axiom::
1137 sage: set_random_seed()
1138 sage: J = ConcreteEuclideanJordanAlgebra.random_instance()
1139 sage: x = J.random_element()
1140 sage: x.operator().is_self_adjoint()
1145 def _max_random_instance_size():
1147 Return an integer "size" that is an upper bound on the size of
1148 this algebra when it is used in a random test
1149 case. Unfortunately, the term "size" is ambiguous -- when
1150 dealing with `R^n` under either the Hadamard or Jordan spin
1151 product, the "size" refers to the dimension `n`. When dealing
1152 with a matrix algebra (real symmetric or complex/quaternion
1153 Hermitian), it refers to the size of the matrix, which is far
1154 less than the dimension of the underlying vector space.
1156 This method must be implemented in each subclass.
1158 raise NotImplementedError
1161 def random_instance(cls
, field
=AA
, **kwargs
):
1163 Return a random instance of this type of algebra.
1165 This method should be implemented in each subclass.
1167 from sage
.misc
.prandom
import choice
1168 eja_class
= choice(cls
.__subclasses
__())
1169 return eja_class
.random_instance(field
)
1172 class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra
):
1174 def __init__(self
, field
, basis
, normalize_basis
=True, **kwargs
):
1176 Compared to the superclass constructor, we take a basis instead of
1177 a multiplication table because the latter can be computed in terms
1178 of the former when the product is known (like it is here).
1180 # Used in this class's fast _charpoly_coefficients() override.
1181 self
._basis
_normalizers
= None
1183 # We're going to loop through this a few times, so now's a good
1184 # time to ensure that it isn't a generator expression.
1185 basis
= tuple(basis
)
1187 algebra_dim
= len(basis
)
1188 degree
= 0 # size of the matrices
1190 degree
= basis
[0].nrows()
1192 if algebra_dim
> 1 and normalize_basis
:
1193 # We'll need sqrt(2) to normalize the basis, and this
1194 # winds up in the multiplication table, so the whole
1195 # algebra needs to be over the field extension.
1196 R
= PolynomialRing(field
, 'z')
1199 if p
.is_irreducible():
1200 field
= field
.extension(p
, 'sqrt2', embedding
=RLF(2).sqrt())
1201 basis
= tuple( s
.change_ring(field
) for s
in basis
)
1202 self
._basis
_normalizers
= tuple(
1203 ~
(self
.matrix_inner_product(s
,s
).sqrt()) for s
in basis
)
1204 basis
= tuple(s
*c
for (s
,c
) in zip(basis
,self
._basis
_normalizers
))
1206 # Now compute the multiplication and inner product tables.
1207 # We have to do this *after* normalizing the basis, because
1208 # scaling affects the answers.
1209 V
= VectorSpace(field
, degree
**2)
1210 W
= V
.span_of_basis( _mat2vec(s
) for s
in basis
)
1211 mult_table
= [[W
.zero() for j
in range(algebra_dim
)]
1212 for i
in range(algebra_dim
)]
1213 ip_table
= [[W
.zero() for j
in range(algebra_dim
)]
1214 for i
in range(algebra_dim
)]
1215 for i
in range(algebra_dim
):
1216 for j
in range(algebra_dim
):
1217 mat_entry
= (basis
[i
]*basis
[j
] + basis
[j
]*basis
[i
])/2
1218 mult_table
[i
][j
] = W
.coordinate_vector(_mat2vec(mat_entry
))
1221 # HACK: ignore the error here if we don't need the
1222 # inner product (as is the case when we construct
1223 # a dummy QQ-algebra for fast charpoly coefficients.
1224 ip_table
[i
][j
] = self
.matrix_inner_product(basis
[i
],
1231 self
._inner
_product
_matrix
= matrix(field
,ip_table
)
1235 super(MatrixEuclideanJordanAlgebra
, self
).__init
__(field
,
1240 if algebra_dim
== 0:
1241 self
.one
.set_cache(self
.zero())
1243 n
= basis
[0].nrows()
1244 # The identity wrt (A,B) -> (AB + BA)/2 is independent of the
1245 # details of this algebra.
1246 self
.one
.set_cache(self(matrix
.identity(field
,n
)))
1250 def _charpoly_coefficients(self
):
1252 Override the parent method with something that tries to compute
1253 over a faster (non-extension) field.
1255 if self
._basis
_normalizers
is None or self
.base_ring() is QQ
:
1256 # We didn't normalize, or the basis we started with had
1257 # entries in a nice field already. Just compute the thing.
1258 return super(MatrixEuclideanJordanAlgebra
, self
)._charpoly
_coefficients
()
1260 basis
= ( (b
/n
) for (b
,n
) in zip(self
.matrix_basis(),
1261 self
._basis
_normalizers
) )
1263 # Do this over the rationals and convert back at the end.
1264 # Only works because we know the entries of the basis are
1265 # integers. The argument ``check_axioms=False`` is required
1266 # because the trace inner-product method for this
1267 # class is a stub and can't actually be checked.
1268 J
= MatrixEuclideanJordanAlgebra(QQ
,
1270 normalize_basis
=False,
1273 a
= J
._charpoly
_coefficients
()
1275 # Unfortunately, changing the basis does change the
1276 # coefficients of the characteristic polynomial, but since
1277 # these are really the coefficients of the "characteristic
1278 # polynomial of" function, everything is still nice and
1279 # unevaluated. It's therefore "obvious" how scaling the
1280 # basis affects the coordinate variables X1, X2, et
1281 # cetera. Scaling the first basis vector up by "n" adds a
1282 # factor of 1/n into every "X1" term, for example. So here
1283 # we simply undo the basis_normalizer scaling that we
1284 # performed earlier.
1286 # The a[0] access here is safe because trivial algebras
1287 # won't have any basis normalizers and therefore won't
1288 # make it to this "else" branch.
1289 XS
= a
[0].parent().gens()
1290 subs_dict
= { XS
[i
]: self
._basis
_normalizers
[i
]*XS
[i
]
1291 for i
in range(len(XS
)) }
1292 return tuple( a_i
.subs(subs_dict
) for a_i
in a
)
1298 Embed the matrix ``M`` into a space of real matrices.
1300 The matrix ``M`` can have entries in any field at the moment:
1301 the real numbers, complex numbers, or quaternions. And although
1302 they are not a field, we can probably support octonions at some
1303 point, too. This function returns a real matrix that "acts like"
1304 the original with respect to matrix multiplication; i.e.
1306 real_embed(M*N) = real_embed(M)*real_embed(N)
1309 raise NotImplementedError
1313 def real_unembed(M
):
1315 The inverse of :meth:`real_embed`.
1317 raise NotImplementedError
1320 def matrix_inner_product(cls
,X
,Y
):
1321 Xu
= cls
.real_unembed(X
)
1322 Yu
= cls
.real_unembed(Y
)
1323 tr
= (Xu
*Yu
).trace()
1326 # Works in QQ, AA, RDF, et cetera.
1328 except AttributeError:
1329 # A quaternion doesn't have a real() method, but does
1330 # have coefficient_tuple() method that returns the
1331 # coefficients of 1, i, j, and k -- in that order.
1332 return tr
.coefficient_tuple()[0]
1335 class RealMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1339 The identity function, for embedding real matrices into real
1345 def real_unembed(M
):
1347 The identity function, for unembedding real matrices from real
1353 class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra
,
1354 ConcreteEuclideanJordanAlgebra
):
1356 The rank-n simple EJA consisting of real symmetric n-by-n
1357 matrices, the usual symmetric Jordan product, and the trace inner
1358 product. It has dimension `(n^2 + n)/2` over the reals.
1362 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1366 sage: J = RealSymmetricEJA(2)
1367 sage: e0, e1, e2 = J.gens()
1375 In theory, our "field" can be any subfield of the reals::
1377 sage: RealSymmetricEJA(2, RDF)
1378 Euclidean Jordan algebra of dimension 3 over Real Double Field
1379 sage: RealSymmetricEJA(2, RR)
1380 Euclidean Jordan algebra of dimension 3 over Real Field with
1381 53 bits of precision
1385 The dimension of this algebra is `(n^2 + n) / 2`::
1387 sage: set_random_seed()
1388 sage: n_max = RealSymmetricEJA._max_random_instance_size()
1389 sage: n = ZZ.random_element(1, n_max)
1390 sage: J = RealSymmetricEJA(n)
1391 sage: J.dimension() == (n^2 + n)/2
1394 The Jordan multiplication is what we think it is::
1396 sage: set_random_seed()
1397 sage: J = RealSymmetricEJA.random_instance()
1398 sage: x,y = J.random_elements(2)
1399 sage: actual = (x*y).to_matrix()
1400 sage: X = x.to_matrix()
1401 sage: Y = y.to_matrix()
1402 sage: expected = (X*Y + Y*X)/2
1403 sage: actual == expected
1405 sage: J(expected) == x*y
1408 We can change the generator prefix::
1410 sage: RealSymmetricEJA(3, prefix='q').gens()
1411 (q0, q1, q2, q3, q4, q5)
1413 We can construct the (trivial) algebra of rank zero::
1415 sage: RealSymmetricEJA(0)
1416 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1420 def _denormalized_basis(cls
, n
, field
):
1422 Return a basis for the space of real symmetric n-by-n matrices.
1426 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1430 sage: set_random_seed()
1431 sage: n = ZZ.random_element(1,5)
1432 sage: B = RealSymmetricEJA._denormalized_basis(n,QQ)
1433 sage: all( M.is_symmetric() for M in B)
1437 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1441 for j
in range(i
+1):
1442 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
1446 Sij
= Eij
+ Eij
.transpose()
1452 def _max_random_instance_size():
1453 return 4 # Dimension 10
1456 def random_instance(cls
, field
=AA
, **kwargs
):
1458 Return a random instance of this type of algebra.
1460 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1461 return cls(n
, field
, **kwargs
)
1463 def __init__(self
, n
, field
=AA
, **kwargs
):
1464 basis
= self
._denormalized
_basis
(n
, field
)
1465 super(RealSymmetricEJA
, self
).__init
__(field
,
1469 self
.rank
.set_cache(n
)
1472 class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1476 Embed the n-by-n complex matrix ``M`` into the space of real
1477 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1478 bi` to the block matrix ``[[a,b],[-b,a]]``.
1482 sage: from mjo.eja.eja_algebra import \
1483 ....: ComplexMatrixEuclideanJordanAlgebra
1487 sage: F = QuadraticField(-1, 'I')
1488 sage: x1 = F(4 - 2*i)
1489 sage: x2 = F(1 + 2*i)
1492 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1493 sage: ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1502 Embedding is a homomorphism (isomorphism, in fact)::
1504 sage: set_random_seed()
1505 sage: n = ZZ.random_element(3)
1506 sage: F = QuadraticField(-1, 'I')
1507 sage: X = random_matrix(F, n)
1508 sage: Y = random_matrix(F, n)
1509 sage: Xe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X)
1510 sage: Ye = ComplexMatrixEuclideanJordanAlgebra.real_embed(Y)
1511 sage: XYe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1518 raise ValueError("the matrix 'M' must be square")
1520 # We don't need any adjoined elements...
1521 field
= M
.base_ring().base_ring()
1525 a
= z
.list()[0] # real part, I guess
1526 b
= z
.list()[1] # imag part, I guess
1527 blocks
.append(matrix(field
, 2, [[a
,b
],[-b
,a
]]))
1529 return matrix
.block(field
, n
, blocks
)
1533 def real_unembed(M
):
1535 The inverse of _embed_complex_matrix().
1539 sage: from mjo.eja.eja_algebra import \
1540 ....: ComplexMatrixEuclideanJordanAlgebra
1544 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1545 ....: [-2, 1, -4, 3],
1546 ....: [ 9, 10, 11, 12],
1547 ....: [-10, 9, -12, 11] ])
1548 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(A)
1550 [ 10*I + 9 12*I + 11]
1554 Unembedding is the inverse of embedding::
1556 sage: set_random_seed()
1557 sage: F = QuadraticField(-1, 'I')
1558 sage: M = random_matrix(F, 3)
1559 sage: Me = ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1560 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1566 raise ValueError("the matrix 'M' must be square")
1567 if not n
.mod(2).is_zero():
1568 raise ValueError("the matrix 'M' must be a complex embedding")
1570 # If "M" was normalized, its base ring might have roots
1571 # adjoined and they can stick around after unembedding.
1572 field
= M
.base_ring()
1573 R
= PolynomialRing(field
, 'z')
1576 # Sage doesn't know how to embed AA into QQbar, i.e. how
1577 # to adjoin sqrt(-1) to AA.
1580 F
= field
.extension(z
**2 + 1, 'I', embedding
=CLF(-1).sqrt())
1583 # Go top-left to bottom-right (reading order), converting every
1584 # 2-by-2 block we see to a single complex element.
1586 for k
in range(n
/2):
1587 for j
in range(n
/2):
1588 submat
= M
[2*k
:2*k
+2,2*j
:2*j
+2]
1589 if submat
[0,0] != submat
[1,1]:
1590 raise ValueError('bad on-diagonal submatrix')
1591 if submat
[0,1] != -submat
[1,0]:
1592 raise ValueError('bad off-diagonal submatrix')
1593 z
= submat
[0,0] + submat
[0,1]*i
1596 return matrix(F
, n
/2, elements
)
1600 def matrix_inner_product(cls
,X
,Y
):
1602 Compute a matrix inner product in this algebra directly from
1607 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1611 This gives the same answer as the slow, default method implemented
1612 in :class:`MatrixEuclideanJordanAlgebra`::
1614 sage: set_random_seed()
1615 sage: J = ComplexHermitianEJA.random_instance()
1616 sage: x,y = J.random_elements(2)
1617 sage: Xe = x.to_matrix()
1618 sage: Ye = y.to_matrix()
1619 sage: X = ComplexHermitianEJA.real_unembed(Xe)
1620 sage: Y = ComplexHermitianEJA.real_unembed(Ye)
1621 sage: expected = (X*Y).trace().real()
1622 sage: actual = ComplexHermitianEJA.matrix_inner_product(Xe,Ye)
1623 sage: actual == expected
1627 return RealMatrixEuclideanJordanAlgebra
.matrix_inner_product(X
,Y
)/2
1630 class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra
,
1631 ConcreteEuclideanJordanAlgebra
):
1633 The rank-n simple EJA consisting of complex Hermitian n-by-n
1634 matrices over the real numbers, the usual symmetric Jordan product,
1635 and the real-part-of-trace inner product. It has dimension `n^2` over
1640 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1644 In theory, our "field" can be any subfield of the reals::
1646 sage: ComplexHermitianEJA(2, RDF)
1647 Euclidean Jordan algebra of dimension 4 over Real Double Field
1648 sage: ComplexHermitianEJA(2, RR)
1649 Euclidean Jordan algebra of dimension 4 over Real Field with
1650 53 bits of precision
1654 The dimension of this algebra is `n^2`::
1656 sage: set_random_seed()
1657 sage: n_max = ComplexHermitianEJA._max_random_instance_size()
1658 sage: n = ZZ.random_element(1, n_max)
1659 sage: J = ComplexHermitianEJA(n)
1660 sage: J.dimension() == n^2
1663 The Jordan multiplication is what we think it is::
1665 sage: set_random_seed()
1666 sage: J = ComplexHermitianEJA.random_instance()
1667 sage: x,y = J.random_elements(2)
1668 sage: actual = (x*y).to_matrix()
1669 sage: X = x.to_matrix()
1670 sage: Y = y.to_matrix()
1671 sage: expected = (X*Y + Y*X)/2
1672 sage: actual == expected
1674 sage: J(expected) == x*y
1677 We can change the generator prefix::
1679 sage: ComplexHermitianEJA(2, prefix='z').gens()
1682 We can construct the (trivial) algebra of rank zero::
1684 sage: ComplexHermitianEJA(0)
1685 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1690 def _denormalized_basis(cls
, n
, field
):
1692 Returns a basis for the space of complex Hermitian n-by-n matrices.
1694 Why do we embed these? Basically, because all of numerical linear
1695 algebra assumes that you're working with vectors consisting of `n`
1696 entries from a field and scalars from the same field. There's no way
1697 to tell SageMath that (for example) the vectors contain complex
1698 numbers, while the scalar field is real.
1702 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1706 sage: set_random_seed()
1707 sage: n = ZZ.random_element(1,5)
1708 sage: field = QuadraticField(2, 'sqrt2')
1709 sage: B = ComplexHermitianEJA._denormalized_basis(n, field)
1710 sage: all( M.is_symmetric() for M in B)
1714 R
= PolynomialRing(field
, 'z')
1716 F
= field
.extension(z
**2 + 1, 'I')
1719 # This is like the symmetric case, but we need to be careful:
1721 # * We want conjugate-symmetry, not just symmetry.
1722 # * The diagonal will (as a result) be real.
1726 for j
in range(i
+1):
1727 Eij
= matrix(F
, n
, lambda k
,l
: k
==i
and l
==j
)
1729 Sij
= cls
.real_embed(Eij
)
1732 # The second one has a minus because it's conjugated.
1733 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
1735 Sij_imag
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
1738 # Since we embedded these, we can drop back to the "field" that we
1739 # started with instead of the complex extension "F".
1740 return ( s
.change_ring(field
) for s
in S
)
1743 def __init__(self
, n
, field
=AA
, **kwargs
):
1744 basis
= self
._denormalized
_basis
(n
,field
)
1745 super(ComplexHermitianEJA
,self
).__init
__(field
,
1749 self
.rank
.set_cache(n
)
1752 def _max_random_instance_size():
1753 return 3 # Dimension 9
1756 def random_instance(cls
, field
=AA
, **kwargs
):
1758 Return a random instance of this type of algebra.
1760 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1761 return cls(n
, field
, **kwargs
)
1763 class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1767 Embed the n-by-n quaternion matrix ``M`` into the space of real
1768 matrices of size 4n-by-4n by first sending each quaternion entry `z
1769 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
1770 c+di],[-c + di, a-bi]]`, and then embedding those into a real
1775 sage: from mjo.eja.eja_algebra import \
1776 ....: QuaternionMatrixEuclideanJordanAlgebra
1780 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1781 sage: i,j,k = Q.gens()
1782 sage: x = 1 + 2*i + 3*j + 4*k
1783 sage: M = matrix(Q, 1, [[x]])
1784 sage: QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1790 Embedding is a homomorphism (isomorphism, in fact)::
1792 sage: set_random_seed()
1793 sage: n = ZZ.random_element(2)
1794 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1795 sage: X = random_matrix(Q, n)
1796 sage: Y = random_matrix(Q, n)
1797 sage: Xe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X)
1798 sage: Ye = QuaternionMatrixEuclideanJordanAlgebra.real_embed(Y)
1799 sage: XYe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1804 quaternions
= M
.base_ring()
1807 raise ValueError("the matrix 'M' must be square")
1809 F
= QuadraticField(-1, 'I')
1814 t
= z
.coefficient_tuple()
1819 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
1820 [-c
+ d
*i
, a
- b
*i
]])
1821 realM
= ComplexMatrixEuclideanJordanAlgebra
.real_embed(cplxM
)
1822 blocks
.append(realM
)
1824 # We should have real entries by now, so use the realest field
1825 # we've got for the return value.
1826 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
1831 def real_unembed(M
):
1833 The inverse of _embed_quaternion_matrix().
1837 sage: from mjo.eja.eja_algebra import \
1838 ....: QuaternionMatrixEuclideanJordanAlgebra
1842 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
1843 ....: [-2, 1, -4, 3],
1844 ....: [-3, 4, 1, -2],
1845 ....: [-4, -3, 2, 1]])
1846 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(M)
1847 [1 + 2*i + 3*j + 4*k]
1851 Unembedding is the inverse of embedding::
1853 sage: set_random_seed()
1854 sage: Q = QuaternionAlgebra(QQ, -1, -1)
1855 sage: M = random_matrix(Q, 3)
1856 sage: Me = QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1857 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1863 raise ValueError("the matrix 'M' must be square")
1864 if not n
.mod(4).is_zero():
1865 raise ValueError("the matrix 'M' must be a quaternion embedding")
1867 # Use the base ring of the matrix to ensure that its entries can be
1868 # multiplied by elements of the quaternion algebra.
1869 field
= M
.base_ring()
1870 Q
= QuaternionAlgebra(field
,-1,-1)
1873 # Go top-left to bottom-right (reading order), converting every
1874 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
1877 for l
in range(n
/4):
1878 for m
in range(n
/4):
1879 submat
= ComplexMatrixEuclideanJordanAlgebra
.real_unembed(
1880 M
[4*l
:4*l
+4,4*m
:4*m
+4] )
1881 if submat
[0,0] != submat
[1,1].conjugate():
1882 raise ValueError('bad on-diagonal submatrix')
1883 if submat
[0,1] != -submat
[1,0].conjugate():
1884 raise ValueError('bad off-diagonal submatrix')
1885 z
= submat
[0,0].real()
1886 z
+= submat
[0,0].imag()*i
1887 z
+= submat
[0,1].real()*j
1888 z
+= submat
[0,1].imag()*k
1891 return matrix(Q
, n
/4, elements
)
1895 def matrix_inner_product(cls
,X
,Y
):
1897 Compute a matrix inner product in this algebra directly from
1902 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1906 This gives the same answer as the slow, default method implemented
1907 in :class:`MatrixEuclideanJordanAlgebra`::
1909 sage: set_random_seed()
1910 sage: J = QuaternionHermitianEJA.random_instance()
1911 sage: x,y = J.random_elements(2)
1912 sage: Xe = x.to_matrix()
1913 sage: Ye = y.to_matrix()
1914 sage: X = QuaternionHermitianEJA.real_unembed(Xe)
1915 sage: Y = QuaternionHermitianEJA.real_unembed(Ye)
1916 sage: expected = (X*Y).trace().coefficient_tuple()[0]
1917 sage: actual = QuaternionHermitianEJA.matrix_inner_product(Xe,Ye)
1918 sage: actual == expected
1922 return RealMatrixEuclideanJordanAlgebra
.matrix_inner_product(X
,Y
)/4
1925 class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra
,
1926 ConcreteEuclideanJordanAlgebra
):
1928 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
1929 matrices, the usual symmetric Jordan product, and the
1930 real-part-of-trace inner product. It has dimension `2n^2 - n` over
1935 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1939 In theory, our "field" can be any subfield of the reals::
1941 sage: QuaternionHermitianEJA(2, RDF)
1942 Euclidean Jordan algebra of dimension 6 over Real Double Field
1943 sage: QuaternionHermitianEJA(2, RR)
1944 Euclidean Jordan algebra of dimension 6 over Real Field with
1945 53 bits of precision
1949 The dimension of this algebra is `2*n^2 - n`::
1951 sage: set_random_seed()
1952 sage: n_max = QuaternionHermitianEJA._max_random_instance_size()
1953 sage: n = ZZ.random_element(1, n_max)
1954 sage: J = QuaternionHermitianEJA(n)
1955 sage: J.dimension() == 2*(n^2) - n
1958 The Jordan multiplication is what we think it is::
1960 sage: set_random_seed()
1961 sage: J = QuaternionHermitianEJA.random_instance()
1962 sage: x,y = J.random_elements(2)
1963 sage: actual = (x*y).to_matrix()
1964 sage: X = x.to_matrix()
1965 sage: Y = y.to_matrix()
1966 sage: expected = (X*Y + Y*X)/2
1967 sage: actual == expected
1969 sage: J(expected) == x*y
1972 We can change the generator prefix::
1974 sage: QuaternionHermitianEJA(2, prefix='a').gens()
1975 (a0, a1, a2, a3, a4, a5)
1977 We can construct the (trivial) algebra of rank zero::
1979 sage: QuaternionHermitianEJA(0)
1980 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1984 def _denormalized_basis(cls
, n
, field
):
1986 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
1988 Why do we embed these? Basically, because all of numerical
1989 linear algebra assumes that you're working with vectors consisting
1990 of `n` entries from a field and scalars from the same field. There's
1991 no way to tell SageMath that (for example) the vectors contain
1992 complex numbers, while the scalar field is real.
1996 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2000 sage: set_random_seed()
2001 sage: n = ZZ.random_element(1,5)
2002 sage: B = QuaternionHermitianEJA._denormalized_basis(n,QQ)
2003 sage: all( M.is_symmetric() for M in B )
2007 Q
= QuaternionAlgebra(QQ
,-1,-1)
2010 # This is like the symmetric case, but we need to be careful:
2012 # * We want conjugate-symmetry, not just symmetry.
2013 # * The diagonal will (as a result) be real.
2017 for j
in range(i
+1):
2018 Eij
= matrix(Q
, n
, lambda k
,l
: k
==i
and l
==j
)
2020 Sij
= cls
.real_embed(Eij
)
2023 # The second, third, and fourth ones have a minus
2024 # because they're conjugated.
2025 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
2027 Sij_I
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
2029 Sij_J
= cls
.real_embed(J
*Eij
- J
*Eij
.transpose())
2031 Sij_K
= cls
.real_embed(K
*Eij
- K
*Eij
.transpose())
2034 # Since we embedded these, we can drop back to the "field" that we
2035 # started with instead of the quaternion algebra "Q".
2036 return ( s
.change_ring(field
) for s
in S
)
2039 def __init__(self
, n
, field
=AA
, **kwargs
):
2040 basis
= self
._denormalized
_basis
(n
,field
)
2041 super(QuaternionHermitianEJA
,self
).__init
__(field
,
2045 self
.rank
.set_cache(n
)
2048 def _max_random_instance_size():
2050 The maximum rank of a random QuaternionHermitianEJA.
2052 return 2 # Dimension 6
2055 def random_instance(cls
, field
=AA
, **kwargs
):
2057 Return a random instance of this type of algebra.
2059 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2060 return cls(n
, field
, **kwargs
)
2063 class HadamardEJA(RationalBasisEuclideanJordanAlgebra
,
2064 ConcreteEuclideanJordanAlgebra
):
2066 Return the Euclidean Jordan Algebra corresponding to the set
2067 `R^n` under the Hadamard product.
2069 Note: this is nothing more than the Cartesian product of ``n``
2070 copies of the spin algebra. Once Cartesian product algebras
2071 are implemented, this can go.
2075 sage: from mjo.eja.eja_algebra import HadamardEJA
2079 This multiplication table can be verified by hand::
2081 sage: J = HadamardEJA(3)
2082 sage: e0,e1,e2 = J.gens()
2098 We can change the generator prefix::
2100 sage: HadamardEJA(3, prefix='r').gens()
2104 def __init__(self
, n
, field
=AA
, **kwargs
):
2105 V
= VectorSpace(field
, n
)
2106 mult_table
= [ [ V
.gen(i
)*(i
== j
) for j
in range(n
) ]
2109 # Inner products are real numbers and not algebra
2110 # elements, so once we turn the algebra element
2111 # into a vector in inner_product(), we never go
2112 # back. As a result -- contrary to what we do with
2113 # self._multiplication_table -- we store the inner
2114 # product table as a plain old matrix and not as
2115 # an algebra operator.
2116 ip_table
= matrix
.identity(field
,n
)
2117 self
._inner
_product
_matrix
= ip_table
2119 super(HadamardEJA
, self
).__init
__(field
,
2123 self
.rank
.set_cache(n
)
2126 self
.one
.set_cache( self
.zero() )
2128 self
.one
.set_cache( sum(self
.gens()) )
2131 def _max_random_instance_size():
2133 The maximum dimension of a random HadamardEJA.
2138 def random_instance(cls
, field
=AA
, **kwargs
):
2140 Return a random instance of this type of algebra.
2142 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2143 return cls(n
, field
, **kwargs
)
2146 class BilinearFormEJA(RationalBasisEuclideanJordanAlgebra
,
2147 ConcreteEuclideanJordanAlgebra
):
2149 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2150 with the half-trace inner product and jordan product ``x*y =
2151 (<Bx,y>,y_bar>, x0*y_bar + y0*x_bar)`` where `B = 1 \times B22` is
2152 a symmetric positive-definite "bilinear form" matrix. Its
2153 dimension is the size of `B`, and it has rank two in dimensions
2154 larger than two. It reduces to the ``JordanSpinEJA`` when `B` is
2155 the identity matrix of order ``n``.
2157 We insist that the one-by-one upper-left identity block of `B` be
2158 passed in as well so that we can be passed a matrix of size zero
2159 to construct a trivial algebra.
2163 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
2164 ....: JordanSpinEJA)
2168 When no bilinear form is specified, the identity matrix is used,
2169 and the resulting algebra is the Jordan spin algebra::
2171 sage: B = matrix.identity(AA,3)
2172 sage: J0 = BilinearFormEJA(B)
2173 sage: J1 = JordanSpinEJA(3)
2174 sage: J0.multiplication_table() == J0.multiplication_table()
2177 An error is raised if the matrix `B` does not correspond to a
2178 positive-definite bilinear form::
2180 sage: B = matrix.random(QQ,2,3)
2181 sage: J = BilinearFormEJA(B)
2182 Traceback (most recent call last):
2184 ValueError: bilinear form is not positive-definite
2185 sage: B = matrix.zero(QQ,3)
2186 sage: J = BilinearFormEJA(B)
2187 Traceback (most recent call last):
2189 ValueError: bilinear form is not positive-definite
2193 We can create a zero-dimensional algebra::
2195 sage: B = matrix.identity(AA,0)
2196 sage: J = BilinearFormEJA(B)
2200 We can check the multiplication condition given in the Jordan, von
2201 Neumann, and Wigner paper (and also discussed on my "On the
2202 symmetry..." paper). Note that this relies heavily on the standard
2203 choice of basis, as does anything utilizing the bilinear form matrix::
2205 sage: set_random_seed()
2206 sage: n = ZZ.random_element(5)
2207 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2208 sage: B11 = matrix.identity(QQ,1)
2209 sage: B22 = M.transpose()*M
2210 sage: B = block_matrix(2,2,[ [B11,0 ],
2212 sage: J = BilinearFormEJA(B)
2213 sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
2214 sage: V = J.vector_space()
2215 sage: sis = [ J.from_vector(V([0] + (M.inverse()*ei).list()))
2216 ....: for ei in eis ]
2217 sage: actual = [ sis[i]*sis[j]
2218 ....: for i in range(n-1)
2219 ....: for j in range(n-1) ]
2220 sage: expected = [ J.one() if i == j else J.zero()
2221 ....: for i in range(n-1)
2222 ....: for j in range(n-1) ]
2223 sage: actual == expected
2226 def __init__(self
, B
, field
=AA
, **kwargs
):
2229 if not B
.is_positive_definite():
2230 raise ValueError("bilinear form is not positive-definite")
2232 V
= VectorSpace(field
, n
)
2233 mult_table
= [[V
.zero() for j
in range(n
)] for i
in range(n
)]
2242 z0
= (B
*x
).inner_product(y
)
2243 zbar
= y0
*xbar
+ x0
*ybar
2244 z
= V([z0
] + zbar
.list())
2245 mult_table
[i
][j
] = z
2247 # Inner products are real numbers and not algebra
2248 # elements, so once we turn the algebra element
2249 # into a vector in inner_product(), we never go
2250 # back. As a result -- contrary to what we do with
2251 # self._multiplication_table -- we store the inner
2252 # product table as a plain old matrix and not as
2253 # an algebra operator.
2255 self
._inner
_product
_matrix
= ip_table
2257 super(BilinearFormEJA
, self
).__init
__(field
,
2262 # The rank of this algebra is two, unless we're in a
2263 # one-dimensional ambient space (because the rank is bounded
2264 # by the ambient dimension).
2265 self
.rank
.set_cache(min(n
,2))
2268 self
.one
.set_cache( self
.zero() )
2270 self
.one
.set_cache( self
.monomial(0) )
2273 def _max_random_instance_size():
2275 The maximum dimension of a random BilinearFormEJA.
2280 def random_instance(cls
, field
=AA
, **kwargs
):
2282 Return a random instance of this algebra.
2284 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2286 B
= matrix
.identity(field
, n
)
2287 return cls(B
, field
, **kwargs
)
2289 B11
= matrix
.identity(field
,1)
2290 M
= matrix
.random(field
, n
-1)
2291 I
= matrix
.identity(field
, n
-1)
2292 alpha
= field
.zero()
2293 while alpha
.is_zero():
2294 alpha
= field
.random_element().abs()
2295 B22
= M
.transpose()*M
+ alpha
*I
2297 from sage
.matrix
.special
import block_matrix
2298 B
= block_matrix(2,2, [ [B11
, ZZ(0) ],
2301 return cls(B
, field
, **kwargs
)
2304 class JordanSpinEJA(BilinearFormEJA
):
2306 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2307 with the usual inner product and jordan product ``x*y =
2308 (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2313 sage: from mjo.eja.eja_algebra import JordanSpinEJA
2317 This multiplication table can be verified by hand::
2319 sage: J = JordanSpinEJA(4)
2320 sage: e0,e1,e2,e3 = J.gens()
2336 We can change the generator prefix::
2338 sage: JordanSpinEJA(2, prefix='B').gens()
2343 Ensure that we have the usual inner product on `R^n`::
2345 sage: set_random_seed()
2346 sage: J = JordanSpinEJA.random_instance()
2347 sage: x,y = J.random_elements(2)
2348 sage: actual = x.inner_product(y)
2349 sage: expected = x.to_vector().inner_product(y.to_vector())
2350 sage: actual == expected
2354 def __init__(self
, n
, field
=AA
, **kwargs
):
2355 # This is a special case of the BilinearFormEJA with the identity
2356 # matrix as its bilinear form.
2357 B
= matrix
.identity(field
, n
)
2358 super(JordanSpinEJA
, self
).__init
__(B
, field
, **kwargs
)
2361 def _max_random_instance_size():
2363 The maximum dimension of a random JordanSpinEJA.
2368 def random_instance(cls
, field
=AA
, **kwargs
):
2370 Return a random instance of this type of algebra.
2372 Needed here to override the implementation for ``BilinearFormEJA``.
2374 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2375 return cls(n
, field
, **kwargs
)
2378 class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra
,
2379 ConcreteEuclideanJordanAlgebra
):
2381 The trivial Euclidean Jordan algebra consisting of only a zero element.
2385 sage: from mjo.eja.eja_algebra import TrivialEJA
2389 sage: J = TrivialEJA()
2396 sage: 7*J.one()*12*J.one()
2398 sage: J.one().inner_product(J.one())
2400 sage: J.one().norm()
2402 sage: J.one().subalgebra_generated_by()
2403 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2408 def __init__(self
, field
=AA
, **kwargs
):
2410 self
._inner
_product
_matrix
= matrix(field
,0)
2411 super(TrivialEJA
, self
).__init
__(field
,
2415 # The rank is zero using my definition, namely the dimension of the
2416 # largest subalgebra generated by any element.
2417 self
.rank
.set_cache(0)
2418 self
.one
.set_cache( self
.zero() )
2421 def random_instance(cls
, field
=AA
, **kwargs
):
2422 # We don't take a "size" argument so the superclass method is
2423 # inappropriate for us.
2424 return cls(field
, **kwargs
)
2426 class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra
):
2428 The external (orthogonal) direct sum of two other Euclidean Jordan
2429 algebras. Essentially the Cartesian product of its two factors.
2430 Every Euclidean Jordan algebra decomposes into an orthogonal
2431 direct sum of simple Euclidean Jordan algebras, so no generality
2432 is lost by providing only this construction.
2436 sage: from mjo.eja.eja_algebra import (random_eja,
2438 ....: RealSymmetricEJA,
2443 sage: J1 = HadamardEJA(2)
2444 sage: J2 = RealSymmetricEJA(3)
2445 sage: J = DirectSumEJA(J1,J2)
2453 The external direct sum construction is only valid when the two factors
2454 have the same base ring; an error is raised otherwise::
2456 sage: set_random_seed()
2457 sage: J1 = random_eja(AA)
2458 sage: J2 = random_eja(QQ)
2459 sage: J = DirectSumEJA(J1,J2)
2460 Traceback (most recent call last):
2462 ValueError: algebras must share the same base field
2465 def __init__(self
, J1
, J2
, **kwargs
):
2466 if J1
.base_ring() != J2
.base_ring():
2467 raise ValueError("algebras must share the same base field")
2468 field
= J1
.base_ring()
2470 self
._factors
= (J1
, J2
)
2474 V
= VectorSpace(field
, n
)
2475 mult_table
= [ [ V
.zero() for j
in range(n
) ]
2479 p
= (J1
.monomial(i
)*J1
.monomial(j
)).to_vector()
2480 mult_table
[i
][j
] = V(p
.list() + [field
.zero()]*n2
)
2484 p
= (J2
.monomial(i
)*J2
.monomial(j
)).to_vector()
2485 mult_table
[n1
+i
][n1
+j
] = V([field
.zero()]*n1
+ p
.list())
2487 super(DirectSumEJA
, self
).__init
__(field
,
2491 self
.rank
.set_cache(J1
.rank() + J2
.rank())
2496 Return the pair of this algebra's factors.
2500 sage: from mjo.eja.eja_algebra import (HadamardEJA,
2501 ....: JordanSpinEJA,
2506 sage: J1 = HadamardEJA(2,QQ)
2507 sage: J2 = JordanSpinEJA(3,QQ)
2508 sage: J = DirectSumEJA(J1,J2)
2510 (Euclidean Jordan algebra of dimension 2 over Rational Field,
2511 Euclidean Jordan algebra of dimension 3 over Rational Field)
2514 return self
._factors
2516 def projections(self
):
2518 Return a pair of projections onto this algebra's factors.
2522 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
2523 ....: ComplexHermitianEJA,
2528 sage: J1 = JordanSpinEJA(2)
2529 sage: J2 = ComplexHermitianEJA(2)
2530 sage: J = DirectSumEJA(J1,J2)
2531 sage: (pi_left, pi_right) = J.projections()
2532 sage: J.one().to_vector()
2534 sage: pi_left(J.one()).to_vector()
2536 sage: pi_right(J.one()).to_vector()
2540 (J1
,J2
) = self
.factors()
2543 V_basis
= self
.vector_space().basis()
2544 # Need to specify the dimensions explicitly so that we don't
2545 # wind up with a zero-by-zero matrix when we want e.g. a
2546 # zero-by-two matrix (important for composing things).
2547 P1
= matrix(self
.base_ring(), m
, m
+n
, V_basis
[:m
])
2548 P2
= matrix(self
.base_ring(), n
, m
+n
, V_basis
[m
:])
2549 pi_left
= FiniteDimensionalEuclideanJordanAlgebraOperator(self
,J1
,P1
)
2550 pi_right
= FiniteDimensionalEuclideanJordanAlgebraOperator(self
,J2
,P2
)
2551 return (pi_left
, pi_right
)
2553 def inclusions(self
):
2555 Return the pair of inclusion maps from our factors into us.
2559 sage: from mjo.eja.eja_algebra import (random_eja,
2560 ....: JordanSpinEJA,
2561 ....: RealSymmetricEJA,
2566 sage: J1 = JordanSpinEJA(3)
2567 sage: J2 = RealSymmetricEJA(2)
2568 sage: J = DirectSumEJA(J1,J2)
2569 sage: (iota_left, iota_right) = J.inclusions()
2570 sage: iota_left(J1.zero()) == J.zero()
2572 sage: iota_right(J2.zero()) == J.zero()
2574 sage: J1.one().to_vector()
2576 sage: iota_left(J1.one()).to_vector()
2578 sage: J2.one().to_vector()
2580 sage: iota_right(J2.one()).to_vector()
2582 sage: J.one().to_vector()
2587 Composing a projection with the corresponding inclusion should
2588 produce the identity map, and mismatching them should produce
2591 sage: set_random_seed()
2592 sage: J1 = random_eja()
2593 sage: J2 = random_eja()
2594 sage: J = DirectSumEJA(J1,J2)
2595 sage: (iota_left, iota_right) = J.inclusions()
2596 sage: (pi_left, pi_right) = J.projections()
2597 sage: pi_left*iota_left == J1.one().operator()
2599 sage: pi_right*iota_right == J2.one().operator()
2601 sage: (pi_left*iota_right).is_zero()
2603 sage: (pi_right*iota_left).is_zero()
2607 (J1
,J2
) = self
.factors()
2610 V_basis
= self
.vector_space().basis()
2611 # Need to specify the dimensions explicitly so that we don't
2612 # wind up with a zero-by-zero matrix when we want e.g. a
2613 # two-by-zero matrix (important for composing things).
2614 I1
= matrix
.column(self
.base_ring(), m
, m
+n
, V_basis
[:m
])
2615 I2
= matrix
.column(self
.base_ring(), n
, m
+n
, V_basis
[m
:])
2616 iota_left
= FiniteDimensionalEuclideanJordanAlgebraOperator(J1
,self
,I1
)
2617 iota_right
= FiniteDimensionalEuclideanJordanAlgebraOperator(J2
,self
,I2
)
2618 return (iota_left
, iota_right
)
2620 def inner_product(self
, x
, y
):
2622 The standard Cartesian inner-product.
2624 We project ``x`` and ``y`` onto our factors, and add up the
2625 inner-products from the subalgebras.
2630 sage: from mjo.eja.eja_algebra import (HadamardEJA,
2631 ....: QuaternionHermitianEJA,
2636 sage: J1 = HadamardEJA(3,QQ)
2637 sage: J2 = QuaternionHermitianEJA(2,QQ,normalize_basis=False)
2638 sage: J = DirectSumEJA(J1,J2)
2643 sage: x1.inner_product(x2)
2645 sage: y1.inner_product(y2)
2647 sage: J.one().inner_product(J.one())
2651 (pi_left
, pi_right
) = self
.projections()
2657 return (x1
.inner_product(y1
) + x2
.inner_product(y2
))
2661 random_eja
= ConcreteEuclideanJordanAlgebra
.random_instance