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
= [ [ self
.vector_space().zero()
140 # take advantage of symmetry
143 elt
= self
.from_vector(mult_table
[i
][j
])
144 self
._multiplication
_table
[i
][j
] = elt
145 self
._multiplication
_table
[j
][i
] = elt
148 if not self
._is
_commutative
():
149 raise ValueError("algebra is not commutative")
150 if not self
._is
_jordanian
():
151 raise ValueError("Jordan identity does not hold")
152 if not self
._inner
_product
_is
_associative
():
153 raise ValueError("inner product is not associative")
155 def _element_constructor_(self
, elt
):
157 Construct an element of this algebra from its vector or matrix
160 This gets called only after the parent element _call_ method
161 fails to find a coercion for the argument.
165 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
167 ....: RealSymmetricEJA)
171 The identity in `S^n` is converted to the identity in the EJA::
173 sage: J = RealSymmetricEJA(3)
174 sage: I = matrix.identity(QQ,3)
175 sage: J(I) == J.one()
178 This skew-symmetric matrix can't be represented in the EJA::
180 sage: J = RealSymmetricEJA(3)
181 sage: A = matrix(QQ,3, lambda i,j: i-j)
183 Traceback (most recent call last):
185 ValueError: not an element of this algebra
189 Ensure that we can convert any element of the two non-matrix
190 simple algebras (whose matrix representations are columns)
191 back and forth faithfully::
193 sage: set_random_seed()
194 sage: J = HadamardEJA.random_instance()
195 sage: x = J.random_element()
196 sage: J(x.to_vector().column()) == x
198 sage: J = JordanSpinEJA.random_instance()
199 sage: x = J.random_element()
200 sage: J(x.to_vector().column()) == x
203 msg
= "not an element of this algebra"
205 # The superclass implementation of random_element()
206 # needs to be able to coerce "0" into the algebra.
208 elif elt
in self
.base_ring():
209 # Ensure that no base ring -> algebra coercion is performed
210 # by this method. There's some stupidity in sage that would
211 # otherwise propagate to this method; for example, sage thinks
212 # that the integer 3 belongs to the space of 2-by-2 matrices.
213 raise ValueError(msg
)
215 if elt
not in self
.matrix_space():
216 raise ValueError(msg
)
218 # Thanks for nothing! Matrix spaces aren't vector spaces in
219 # Sage, so we have to figure out its matrix-basis coordinates
220 # ourselves. We use the basis space's ring instead of the
221 # element's ring because the basis space might be an algebraic
222 # closure whereas the base ring of the 3-by-3 identity matrix
223 # could be QQ instead of QQbar.
224 V
= VectorSpace(self
.base_ring(), elt
.nrows()*elt
.ncols())
225 W
= V
.span_of_basis( _mat2vec(s
) for s
in self
.matrix_basis() )
228 coords
= W
.coordinate_vector(_mat2vec(elt
))
229 except ArithmeticError: # vector is not in free module
230 raise ValueError(msg
)
232 return self
.from_vector(coords
)
236 Return a string representation of ``self``.
240 sage: from mjo.eja.eja_algebra import JordanSpinEJA
244 Ensure that it says what we think it says::
246 sage: JordanSpinEJA(2, field=AA)
247 Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
248 sage: JordanSpinEJA(3, field=RDF)
249 Euclidean Jordan algebra of dimension 3 over Real Double Field
252 fmt
= "Euclidean Jordan algebra of dimension {} over {}"
253 return fmt
.format(self
.dimension(), self
.base_ring())
255 def product_on_basis(self
, i
, j
):
256 return self
._multiplication
_table
[i
][j
]
258 def _is_commutative(self
):
260 Whether or not this algebra's multiplication table is commutative.
262 This method should of course always return ``True``, unless
263 this algebra was constructed with ``check_axioms=False`` and
264 passed an invalid multiplication table.
266 return all( self
.product_on_basis(i
,j
) == self
.product_on_basis(i
,j
)
267 for i
in range(self
.dimension())
268 for j
in range(self
.dimension()) )
270 def _is_jordanian(self
):
272 Whether or not this algebra's multiplication table respects the
273 Jordan identity `(x^{2})(xy) = x(x^{2}y)`.
275 We only check one arrangement of `x` and `y`, so for a
276 ``True`` result to be truly true, you should also check
277 :meth:`_is_commutative`. This method should of course always
278 return ``True``, unless this algebra was constructed with
279 ``check_axioms=False`` and passed an invalid multiplication table.
281 return all( (self
.monomial(i
)**2)*(self
.monomial(i
)*self
.monomial(j
))
283 (self
.monomial(i
))*((self
.monomial(i
)**2)*self
.monomial(j
))
284 for i
in range(self
.dimension())
285 for j
in range(self
.dimension()) )
287 def _inner_product_is_associative(self
):
289 Return whether or not this algebra's inner product `B` is
290 associative; that is, whether or not `B(xy,z) = B(x,yz)`.
292 This method should of course always return ``True``, unless
293 this algebra was constructed with ``check_axioms=False`` and
294 passed an invalid multiplication table.
297 # Used to check whether or not something is zero in an inexact
298 # ring. This number is sufficient to allow the construction of
299 # QuaternionHermitianEJA(2, RDF) with check_axioms=True.
302 for i
in range(self
.dimension()):
303 for j
in range(self
.dimension()):
304 for k
in range(self
.dimension()):
308 diff
= (x
*y
).inner_product(z
) - x
.inner_product(y
*z
)
310 if self
.base_ring().is_exact():
314 if diff
.abs() > epsilon
:
320 def characteristic_polynomial_of(self
):
322 Return the algebra's "characteristic polynomial of" function,
323 which is itself a multivariate polynomial that, when evaluated
324 at the coordinates of some algebra element, returns that
325 element's characteristic polynomial.
327 The resulting polynomial has `n+1` variables, where `n` is the
328 dimension of this algebra. The first `n` variables correspond to
329 the coordinates of an algebra element: when evaluated at the
330 coordinates of an algebra element with respect to a certain
331 basis, the result is a univariate polynomial (in the one
332 remaining variable ``t``), namely the characteristic polynomial
337 sage: from mjo.eja.eja_algebra import JordanSpinEJA, TrivialEJA
341 The characteristic polynomial in the spin algebra is given in
342 Alizadeh, Example 11.11::
344 sage: J = JordanSpinEJA(3)
345 sage: p = J.characteristic_polynomial_of(); p
346 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
347 sage: xvec = J.one().to_vector()
351 By definition, the characteristic polynomial is a monic
352 degree-zero polynomial in a rank-zero algebra. Note that
353 Cayley-Hamilton is indeed satisfied since the polynomial
354 ``1`` evaluates to the identity element of the algebra on
357 sage: J = TrivialEJA()
358 sage: J.characteristic_polynomial_of()
365 # The list of coefficient polynomials a_0, a_1, a_2, ..., a_(r-1).
366 a
= self
._charpoly
_coefficients
()
368 # We go to a bit of trouble here to reorder the
369 # indeterminates, so that it's easier to evaluate the
370 # characteristic polynomial at x's coordinates and get back
371 # something in terms of t, which is what we want.
372 S
= PolynomialRing(self
.base_ring(),'t')
376 S
= PolynomialRing(S
, R
.variable_names())
379 return (t
**r
+ sum( a
[k
]*(t
**k
) for k
in range(r
) ))
381 def coordinate_polynomial_ring(self
):
383 The multivariate polynomial ring in which this algebra's
384 :meth:`characteristic_polynomial_of` lives.
388 sage: from mjo.eja.eja_algebra import (HadamardEJA,
389 ....: RealSymmetricEJA)
393 sage: J = HadamardEJA(2)
394 sage: J.coordinate_polynomial_ring()
395 Multivariate Polynomial Ring in X1, X2...
396 sage: J = RealSymmetricEJA(3,QQ)
397 sage: J.coordinate_polynomial_ring()
398 Multivariate Polynomial Ring in X1, X2, X3, X4, X5, X6...
401 var_names
= tuple( "X%d" % z
for z
in range(1, self
.dimension()+1) )
402 return PolynomialRing(self
.base_ring(), var_names
)
404 def inner_product(self
, x
, y
):
406 The inner product associated with this Euclidean Jordan algebra.
408 Defaults to the trace inner product, but can be overridden by
409 subclasses if they are sure that the necessary properties are
414 sage: from mjo.eja.eja_algebra import (random_eja,
416 ....: BilinearFormEJA)
420 Our inner product is "associative," which means the following for
421 a symmetric bilinear form::
423 sage: set_random_seed()
424 sage: J = random_eja()
425 sage: x,y,z = J.random_elements(3)
426 sage: (x*y).inner_product(z) == y.inner_product(x*z)
431 Ensure that this is the usual inner product for the algebras
434 sage: set_random_seed()
435 sage: J = HadamardEJA.random_instance()
436 sage: x,y = J.random_elements(2)
437 sage: actual = x.inner_product(y)
438 sage: expected = x.to_vector().inner_product(y.to_vector())
439 sage: actual == expected
442 Ensure that this is one-half of the trace inner-product in a
443 BilinearFormEJA that isn't just the reals (when ``n`` isn't
444 one). This is in Faraut and Koranyi, and also my "On the
447 sage: set_random_seed()
448 sage: J = BilinearFormEJA.random_instance()
449 sage: n = J.dimension()
450 sage: x = J.random_element()
451 sage: y = J.random_element()
452 sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
455 B
= self
._inner
_product
_matrix
456 return (B
*x
.to_vector()).inner_product(y
.to_vector())
459 def is_trivial(self
):
461 Return whether or not this algebra is trivial.
463 A trivial algebra contains only the zero element.
467 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
472 sage: J = ComplexHermitianEJA(3)
478 sage: J = TrivialEJA()
483 return self
.dimension() == 0
486 def multiplication_table(self
):
488 Return a visual representation of this algebra's multiplication
489 table (on basis elements).
493 sage: from mjo.eja.eja_algebra import JordanSpinEJA
497 sage: J = JordanSpinEJA(4)
498 sage: J.multiplication_table()
499 +----++----+----+----+----+
500 | * || e0 | e1 | e2 | e3 |
501 +====++====+====+====+====+
502 | e0 || e0 | e1 | e2 | e3 |
503 +----++----+----+----+----+
504 | e1 || e1 | e0 | 0 | 0 |
505 +----++----+----+----+----+
506 | e2 || e2 | 0 | e0 | 0 |
507 +----++----+----+----+----+
508 | e3 || e3 | 0 | 0 | e0 |
509 +----++----+----+----+----+
512 M
= list(self
._multiplication
_table
) # copy
513 for i
in range(len(M
)):
514 # M had better be "square"
515 M
[i
] = [self
.monomial(i
)] + M
[i
]
516 M
= [["*"] + list(self
.gens())] + M
517 return table(M
, header_row
=True, header_column
=True, frame
=True)
520 def matrix_basis(self
):
522 Return an (often more natural) representation of this algebras
523 basis as an ordered tuple of matrices.
525 Every finite-dimensional Euclidean Jordan Algebra is a, up to
526 Jordan isomorphism, a direct sum of five simple
527 algebras---four of which comprise Hermitian matrices. And the
528 last type of algebra can of course be thought of as `n`-by-`1`
529 column matrices (ambiguusly called column vectors) to avoid
530 special cases. As a result, matrices (and column vectors) are
531 a natural representation format for Euclidean Jordan algebra
534 But, when we construct an algebra from a basis of matrices,
535 those matrix representations are lost in favor of coordinate
536 vectors *with respect to* that basis. We could eventually
537 convert back if we tried hard enough, but having the original
538 representations handy is valuable enough that we simply store
539 them and return them from this method.
541 Why implement this for non-matrix algebras? Avoiding special
542 cases for the :class:`BilinearFormEJA` pays with simplicity in
543 its own right. But mainly, we would like to be able to assume
544 that elements of a :class:`DirectSumEJA` can be displayed
545 nicely, without having to have special classes for direct sums
546 one of whose components was a matrix algebra.
550 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
551 ....: RealSymmetricEJA)
555 sage: J = RealSymmetricEJA(2)
557 Finite family {0: e0, 1: e1, 2: e2}
558 sage: J.matrix_basis()
560 [1 0] [ 0 0.7071067811865475?] [0 0]
561 [0 0], [0.7071067811865475? 0], [0 1]
566 sage: J = JordanSpinEJA(2)
568 Finite family {0: e0, 1: e1}
569 sage: J.matrix_basis()
575 if self
._matrix
_basis
is None:
576 M
= self
.matrix_space()
577 return tuple( M(b
.to_vector()) for b
in self
.basis() )
579 return self
._matrix
_basis
582 def matrix_space(self
):
584 Return the matrix space in which this algebra's elements live, if
585 we think of them as matrices (including column vectors of the
588 Generally this will be an `n`-by-`1` column-vector space,
589 except when the algebra is trivial. There it's `n`-by-`n`
590 (where `n` is zero), to ensure that two elements of the matrix
591 space (empty matrices) can be multiplied.
593 Matrix algebras override this with something more useful.
595 if self
.is_trivial():
596 return MatrixSpace(self
.base_ring(), 0)
597 elif self
._matrix
_basis
is None or len(self
._matrix
_basis
) == 0:
598 return MatrixSpace(self
.base_ring(), self
.dimension(), 1)
600 return self
._matrix
_basis
[0].matrix_space()
606 Return the unit element of this algebra.
610 sage: from mjo.eja.eja_algebra import (HadamardEJA,
615 sage: J = HadamardEJA(5)
617 e0 + e1 + e2 + e3 + e4
621 The identity element acts like the identity::
623 sage: set_random_seed()
624 sage: J = random_eja()
625 sage: x = J.random_element()
626 sage: J.one()*x == x and x*J.one() == x
629 The matrix of the unit element's operator is the identity::
631 sage: set_random_seed()
632 sage: J = random_eja()
633 sage: actual = J.one().operator().matrix()
634 sage: expected = matrix.identity(J.base_ring(), J.dimension())
635 sage: actual == expected
638 Ensure that the cached unit element (often precomputed by
639 hand) agrees with the computed one::
641 sage: set_random_seed()
642 sage: J = random_eja()
643 sage: cached = J.one()
644 sage: J.one.clear_cache()
645 sage: J.one() == cached
649 # We can brute-force compute the matrices of the operators
650 # that correspond to the basis elements of this algebra.
651 # If some linear combination of those basis elements is the
652 # algebra identity, then the same linear combination of
653 # their matrices has to be the identity matrix.
655 # Of course, matrices aren't vectors in sage, so we have to
656 # appeal to the "long vectors" isometry.
657 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
659 # Now we use basic linear algebra to find the coefficients,
660 # of the matrices-as-vectors-linear-combination, which should
661 # work for the original algebra basis too.
662 A
= matrix(self
.base_ring(), oper_vecs
)
664 # We used the isometry on the left-hand side already, but we
665 # still need to do it for the right-hand side. Recall that we
666 # wanted something that summed to the identity matrix.
667 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
669 # Now if there's an identity element in the algebra, this
670 # should work. We solve on the left to avoid having to
671 # transpose the matrix "A".
672 return self
.from_vector(A
.solve_left(b
))
675 def peirce_decomposition(self
, c
):
677 The Peirce decomposition of this algebra relative to the
680 In the future, this can be extended to a complete system of
681 orthogonal idempotents.
685 - ``c`` -- an idempotent of this algebra.
689 A triple (J0, J5, J1) containing two subalgebras and one subspace
692 - ``J0`` -- the algebra on the eigenspace of ``c.operator()``
693 corresponding to the eigenvalue zero.
695 - ``J5`` -- the eigenspace (NOT a subalgebra) of ``c.operator()``
696 corresponding to the eigenvalue one-half.
698 - ``J1`` -- the algebra on the eigenspace of ``c.operator()``
699 corresponding to the eigenvalue one.
701 These are the only possible eigenspaces for that operator, and this
702 algebra is a direct sum of them. The spaces ``J0`` and ``J1`` are
703 orthogonal, and are subalgebras of this algebra with the appropriate
708 sage: from mjo.eja.eja_algebra import random_eja, RealSymmetricEJA
712 The canonical example comes from the symmetric matrices, which
713 decompose into diagonal and off-diagonal parts::
715 sage: J = RealSymmetricEJA(3)
716 sage: C = matrix(QQ, [ [1,0,0],
720 sage: J0,J5,J1 = J.peirce_decomposition(c)
722 Euclidean Jordan algebra of dimension 1...
724 Vector space of degree 6 and dimension 2...
726 Euclidean Jordan algebra of dimension 3...
727 sage: J0.one().to_matrix()
731 sage: orig_df = AA.options.display_format
732 sage: AA.options.display_format = 'radical'
733 sage: J.from_vector(J5.basis()[0]).to_matrix()
737 sage: J.from_vector(J5.basis()[1]).to_matrix()
741 sage: AA.options.display_format = orig_df
742 sage: J1.one().to_matrix()
749 Every algebra decomposes trivially with respect to its identity
752 sage: set_random_seed()
753 sage: J = random_eja()
754 sage: J0,J5,J1 = J.peirce_decomposition(J.one())
755 sage: J0.dimension() == 0 and J5.dimension() == 0
757 sage: J1.superalgebra() == J and J1.dimension() == J.dimension()
760 The decomposition is into eigenspaces, and its components are
761 therefore necessarily orthogonal. Moreover, the identity
762 elements in the two subalgebras are the projections onto their
763 respective subspaces of the superalgebra's identity element::
765 sage: set_random_seed()
766 sage: J = random_eja()
767 sage: x = J.random_element()
768 sage: if not J.is_trivial():
769 ....: while x.is_nilpotent():
770 ....: x = J.random_element()
771 sage: c = x.subalgebra_idempotent()
772 sage: J0,J5,J1 = J.peirce_decomposition(c)
774 sage: for (w,y,z) in zip(J0.basis(), J5.basis(), J1.basis()):
775 ....: w = w.superalgebra_element()
776 ....: y = J.from_vector(y)
777 ....: z = z.superalgebra_element()
778 ....: ipsum += w.inner_product(y).abs()
779 ....: ipsum += w.inner_product(z).abs()
780 ....: ipsum += y.inner_product(z).abs()
783 sage: J1(c) == J1.one()
785 sage: J0(J.one() - c) == J0.one()
789 if not c
.is_idempotent():
790 raise ValueError("element is not idempotent: %s" % c
)
792 from mjo
.eja
.eja_subalgebra
import FiniteDimensionalEuclideanJordanSubalgebra
794 # Default these to what they should be if they turn out to be
795 # trivial, because eigenspaces_left() won't return eigenvalues
796 # corresponding to trivial spaces (e.g. it returns only the
797 # eigenspace corresponding to lambda=1 if you take the
798 # decomposition relative to the identity element).
799 trivial
= FiniteDimensionalEuclideanJordanSubalgebra(self
, ())
800 J0
= trivial
# eigenvalue zero
801 J5
= VectorSpace(self
.base_ring(), 0) # eigenvalue one-half
802 J1
= trivial
# eigenvalue one
804 for (eigval
, eigspace
) in c
.operator().matrix().right_eigenspaces():
805 if eigval
== ~
(self
.base_ring()(2)):
808 gens
= tuple( self
.from_vector(b
) for b
in eigspace
.basis() )
809 subalg
= FiniteDimensionalEuclideanJordanSubalgebra(self
,
817 raise ValueError("unexpected eigenvalue: %s" % eigval
)
822 def random_element(self
, thorough
=False):
824 Return a random element of this algebra.
826 Our algebra superclass method only returns a linear
827 combination of at most two basis elements. We instead
828 want the vector space "random element" method that
829 returns a more diverse selection.
833 - ``thorough`` -- (boolean; default False) whether or not we
834 should generate irrational coefficients for the random
835 element when our base ring is irrational; this slows the
836 algebra operations to a crawl, but any truly random method
840 # For a general base ring... maybe we can trust this to do the
841 # right thing? Unlikely, but.
842 V
= self
.vector_space()
843 v
= V
.random_element()
845 if self
.base_ring() is AA
:
846 # The "random element" method of the algebraic reals is
847 # stupid at the moment, and only returns integers between
848 # -2 and 2, inclusive:
850 # https://trac.sagemath.org/ticket/30875
852 # Instead, we implement our own "random vector" method,
853 # and then coerce that into the algebra. We use the vector
854 # space degree here instead of the dimension because a
855 # subalgebra could (for example) be spanned by only two
856 # vectors, each with five coordinates. We need to
857 # generate all five coordinates.
859 v
*= QQbar
.random_element().real()
861 v
*= QQ
.random_element()
863 return self
.from_vector(V
.coordinate_vector(v
))
865 def random_elements(self
, count
, thorough
=False):
867 Return ``count`` random elements as a tuple.
871 - ``thorough`` -- (boolean; default False) whether or not we
872 should generate irrational coefficients for the random
873 elements when our base ring is irrational; this slows the
874 algebra operations to a crawl, but any truly random method
879 sage: from mjo.eja.eja_algebra import JordanSpinEJA
883 sage: J = JordanSpinEJA(3)
884 sage: x,y,z = J.random_elements(3)
885 sage: all( [ x in J, y in J, z in J ])
887 sage: len( J.random_elements(10) ) == 10
891 return tuple( self
.random_element(thorough
)
892 for idx
in range(count
) )
896 def _charpoly_coefficients(self
):
898 The `r` polynomial coefficients of the "characteristic polynomial
902 R
= self
.coordinate_polynomial_ring()
904 F
= R
.fraction_field()
907 # From a result in my book, these are the entries of the
908 # basis representation of L_x.
909 return sum( vars[k
]*self
.monomial(k
).operator().matrix()[i
,j
]
912 L_x
= matrix(F
, n
, n
, L_x_i_j
)
915 if self
.rank
.is_in_cache():
917 # There's no need to pad the system with redundant
918 # columns if we *know* they'll be redundant.
921 # Compute an extra power in case the rank is equal to
922 # the dimension (otherwise, we would stop at x^(r-1)).
923 x_powers
= [ (L_x
**k
)*self
.one().to_vector()
924 for k
in range(n
+1) ]
925 A
= matrix
.column(F
, x_powers
[:n
])
926 AE
= A
.extended_echelon_form()
933 # The theory says that only the first "r" coefficients are
934 # nonzero, and they actually live in the original polynomial
935 # ring and not the fraction field. We negate them because
936 # in the actual characteristic polynomial, they get moved
937 # to the other side where x^r lives.
938 return -A_rref
.solve_right(E
*b
).change_ring(R
)[:r
]
943 Return the rank of this EJA.
945 This is a cached method because we know the rank a priori for
946 all of the algebras we can construct. Thus we can avoid the
947 expensive ``_charpoly_coefficients()`` call unless we truly
948 need to compute the whole characteristic polynomial.
952 sage: from mjo.eja.eja_algebra import (HadamardEJA,
954 ....: RealSymmetricEJA,
955 ....: ComplexHermitianEJA,
956 ....: QuaternionHermitianEJA,
961 The rank of the Jordan spin algebra is always two::
963 sage: JordanSpinEJA(2).rank()
965 sage: JordanSpinEJA(3).rank()
967 sage: JordanSpinEJA(4).rank()
970 The rank of the `n`-by-`n` Hermitian real, complex, or
971 quaternion matrices is `n`::
973 sage: RealSymmetricEJA(4).rank()
975 sage: ComplexHermitianEJA(3).rank()
977 sage: QuaternionHermitianEJA(2).rank()
982 Ensure that every EJA that we know how to construct has a
983 positive integer rank, unless the algebra is trivial in
984 which case its rank will be zero::
986 sage: set_random_seed()
987 sage: J = random_eja()
991 sage: r > 0 or (r == 0 and J.is_trivial())
994 Ensure that computing the rank actually works, since the ranks
995 of all simple algebras are known and will be cached by default::
997 sage: set_random_seed() # long time
998 sage: J = random_eja() # long time
999 sage: caches = J.rank() # long time
1000 sage: J.rank.clear_cache() # long time
1001 sage: J.rank() == cached # long time
1005 return len(self
._charpoly
_coefficients
())
1008 def vector_space(self
):
1010 Return the vector space that underlies this algebra.
1014 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1018 sage: J = RealSymmetricEJA(2)
1019 sage: J.vector_space()
1020 Vector space of dimension 3 over...
1023 return self
.zero().to_vector().parent().ambient_vector_space()
1026 Element
= FiniteDimensionalEuclideanJordanAlgebraElement
1028 class RationalBasisEuclideanJordanAlgebraNg(FiniteDimensionalEuclideanJordanAlgebra
):
1034 orthonormalize
=True,
1041 vector_basis
= basis
1043 from sage
.structure
.element
import is_Matrix
1044 basis_is_matrices
= False
1048 if is_Matrix(basis
[0]):
1049 basis_is_matrices
= True
1050 vector_basis
= tuple( map(_mat2vec
,basis
) )
1051 degree
= basis
[0].nrows()**2
1053 degree
= basis
[0].degree()
1055 V
= VectorSpace(field
, degree
)
1057 # Compute this from "Q" (obtained from Gram-Schmidt) below as
1058 # R = Q.solve_right(A), where the rows of "Q" are the
1059 # orthonormalized vector_basis and and the rows of "A" are the
1060 # original vector_basis.
1061 self
._deorthonormalization
_matrix
= None
1064 from mjo
.eja
.eja_utils
import gram_schmidt
1065 vector_basis
= gram_schmidt(vector_basis
, inner_product
)
1066 W
= V
.span_of_basis( vector_basis
)
1067 if basis_is_matrices
:
1068 from mjo
.eja
.eja_utils
import _vec2mat
1069 basis
= tuple( map(_vec2mat
,vector_basis
) )
1071 W
= V
.span_of_basis( vector_basis
)
1073 mult_table
= [ [0 for i
in range(n
)] for j
in range(n
) ]
1074 ip_table
= [ [0 for i
in range(n
)] for j
in range(n
) ]
1077 for j
in range(i
+1):
1078 # do another mat2vec because the multiplication
1079 # table is in terms of vectors
1080 elt
= _mat2vec(jordan_product(basis
[i
],basis
[j
]))
1081 elt
= W
.coordinate_vector(elt
)
1082 mult_table
[i
][j
] = elt
1083 mult_table
[j
][i
] = elt
1084 ip
= inner_product(basis
[i
],basis
[j
])
1088 self
._inner
_product
_matrix
= matrix(field
,ip_table
)
1090 if basis_is_matrices
:
1094 basis
= tuple( x
.column() for x
in basis
)
1096 super().__init
__(field
,
1100 basis
, # matrix basis
1104 class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra
):
1106 Algebras whose basis consists of vectors with rational
1107 entries. Equivalently, algebras whose multiplication tables
1108 contain only rational coefficients.
1110 When an EJA has a basis that can be made rational, we can speed up
1111 the computation of its characteristic polynomial by doing it over
1112 ``QQ``. All of the named EJA constructors that we provide fall
1116 def _charpoly_coefficients(self
):
1118 Override the parent method with something that tries to compute
1119 over a faster (non-extension) field.
1123 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1127 The base ring of the resulting polynomial coefficients is what
1128 it should be, and not the rationals (unless the algebra was
1129 already over the rationals)::
1131 sage: J = JordanSpinEJA(3)
1132 sage: J._charpoly_coefficients()
1133 (X1^2 - X2^2 - X3^2, -2*X1)
1134 sage: a0 = J._charpoly_coefficients()[0]
1136 Algebraic Real Field
1137 sage: a0.base_ring()
1138 Algebraic Real Field
1141 if self
.base_ring() is QQ
:
1142 # There's no need to construct *another* algebra over the
1143 # rationals if this one is already over the rationals.
1144 superclass
= super(RationalBasisEuclideanJordanAlgebra
, self
)
1145 return superclass
._charpoly
_coefficients
()
1148 tuple(map(lambda x
: x
.to_vector(), ls
))
1149 for ls
in self
._multiplication
_table
1152 # Do the computation over the rationals. The answer will be
1153 # the same, because our basis coordinates are (essentially)
1155 J
= FiniteDimensionalEuclideanJordanAlgebra(QQ
,
1159 a
= J
._charpoly
_coefficients
()
1160 return tuple(map(lambda x
: x
.change_ring(self
.base_ring()), a
))
1163 class ConcreteEuclideanJordanAlgebra
:
1165 A class for the Euclidean Jordan algebras that we know by name.
1167 These are the Jordan algebras whose basis, multiplication table,
1168 rank, and so on are known a priori. More to the point, they are
1169 the Euclidean Jordan algebras for which we are able to conjure up
1170 a "random instance."
1174 sage: from mjo.eja.eja_algebra import ConcreteEuclideanJordanAlgebra
1178 Our basis is normalized with respect to the algebra's inner
1179 product, unless we specify otherwise::
1181 sage: set_random_seed()
1182 sage: J = ConcreteEuclideanJordanAlgebra.random_instance()
1183 sage: all( b.norm() == 1 for b in J.gens() )
1186 Since our basis is orthonormal with respect to the algebra's inner
1187 product, and since we know that this algebra is an EJA, any
1188 left-multiplication operator's matrix will be symmetric because
1189 natural->EJA basis representation is an isometry and within the
1190 EJA the operator is self-adjoint by the Jordan axiom::
1192 sage: set_random_seed()
1193 sage: J = ConcreteEuclideanJordanAlgebra.random_instance()
1194 sage: x = J.random_element()
1195 sage: x.operator().is_self_adjoint()
1200 def _max_random_instance_size():
1202 Return an integer "size" that is an upper bound on the size of
1203 this algebra when it is used in a random test
1204 case. Unfortunately, the term "size" is ambiguous -- when
1205 dealing with `R^n` under either the Hadamard or Jordan spin
1206 product, the "size" refers to the dimension `n`. When dealing
1207 with a matrix algebra (real symmetric or complex/quaternion
1208 Hermitian), it refers to the size of the matrix, which is far
1209 less than the dimension of the underlying vector space.
1211 This method must be implemented in each subclass.
1213 raise NotImplementedError
1216 def random_instance(cls
, field
=AA
, **kwargs
):
1218 Return a random instance of this type of algebra.
1220 This method should be implemented in each subclass.
1222 from sage
.misc
.prandom
import choice
1223 eja_class
= choice(cls
.__subclasses
__())
1224 return eja_class
.random_instance(field
)
1227 class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra
):
1229 def __init__(self
, field
, basis
, normalize_basis
=True, **kwargs
):
1231 Compared to the superclass constructor, we take a basis instead of
1232 a multiplication table because the latter can be computed in terms
1233 of the former when the product is known (like it is here).
1235 # Used in this class's fast _charpoly_coefficients() override.
1236 self
._basis
_normalizers
= None
1238 # We're going to loop through this a few times, so now's a good
1239 # time to ensure that it isn't a generator expression.
1240 basis
= tuple(basis
)
1242 algebra_dim
= len(basis
)
1243 degree
= 0 # size of the matrices
1245 degree
= basis
[0].nrows()
1247 if algebra_dim
> 1 and normalize_basis
:
1248 # We'll need sqrt(2) to normalize the basis, and this
1249 # winds up in the multiplication table, so the whole
1250 # algebra needs to be over the field extension.
1251 R
= PolynomialRing(field
, 'z')
1254 if p
.is_irreducible():
1255 field
= field
.extension(p
, 'sqrt2', embedding
=RLF(2).sqrt())
1256 basis
= tuple( s
.change_ring(field
) for s
in basis
)
1257 self
._basis
_normalizers
= tuple(
1258 ~
(self
.matrix_inner_product(s
,s
).sqrt()) for s
in basis
)
1259 basis
= tuple(s
*c
for (s
,c
) in zip(basis
,self
._basis
_normalizers
))
1261 # Now compute the multiplication and inner product tables.
1262 # We have to do this *after* normalizing the basis, because
1263 # scaling affects the answers.
1264 V
= VectorSpace(field
, degree
**2)
1265 W
= V
.span_of_basis( _mat2vec(s
) for s
in basis
)
1266 mult_table
= [[W
.zero() for j
in range(algebra_dim
)]
1267 for i
in range(algebra_dim
)]
1268 ip_table
= [[W
.zero() for j
in range(algebra_dim
)]
1269 for i
in range(algebra_dim
)]
1270 for i
in range(algebra_dim
):
1271 for j
in range(algebra_dim
):
1272 mat_entry
= (basis
[i
]*basis
[j
] + basis
[j
]*basis
[i
])/2
1273 mult_table
[i
][j
] = W
.coordinate_vector(_mat2vec(mat_entry
))
1276 # HACK: ignore the error here if we don't need the
1277 # inner product (as is the case when we construct
1278 # a dummy QQ-algebra for fast charpoly coefficients.
1279 ip_table
[i
][j
] = self
.matrix_inner_product(basis
[i
],
1286 self
._inner
_product
_matrix
= matrix(field
,ip_table
)
1290 super(MatrixEuclideanJordanAlgebra
, self
).__init
__(field
,
1295 if algebra_dim
== 0:
1296 self
.one
.set_cache(self
.zero())
1298 n
= basis
[0].nrows()
1299 # The identity wrt (A,B) -> (AB + BA)/2 is independent of the
1300 # details of this algebra.
1301 self
.one
.set_cache(self(matrix
.identity(field
,n
)))
1305 def _charpoly_coefficients(self
):
1307 Override the parent method with something that tries to compute
1308 over a faster (non-extension) field.
1310 if self
._basis
_normalizers
is None or self
.base_ring() is QQ
:
1311 # We didn't normalize, or the basis we started with had
1312 # entries in a nice field already. Just compute the thing.
1313 return super(MatrixEuclideanJordanAlgebra
, self
)._charpoly
_coefficients
()
1315 basis
= ( (b
/n
) for (b
,n
) in zip(self
.matrix_basis(),
1316 self
._basis
_normalizers
) )
1318 # Do this over the rationals and convert back at the end.
1319 # Only works because we know the entries of the basis are
1320 # integers. The argument ``check_axioms=False`` is required
1321 # because the trace inner-product method for this
1322 # class is a stub and can't actually be checked.
1323 J
= MatrixEuclideanJordanAlgebra(QQ
,
1325 normalize_basis
=False,
1328 a
= J
._charpoly
_coefficients
()
1330 # Unfortunately, changing the basis does change the
1331 # coefficients of the characteristic polynomial, but since
1332 # these are really the coefficients of the "characteristic
1333 # polynomial of" function, everything is still nice and
1334 # unevaluated. It's therefore "obvious" how scaling the
1335 # basis affects the coordinate variables X1, X2, et
1336 # cetera. Scaling the first basis vector up by "n" adds a
1337 # factor of 1/n into every "X1" term, for example. So here
1338 # we simply undo the basis_normalizer scaling that we
1339 # performed earlier.
1341 # The a[0] access here is safe because trivial algebras
1342 # won't have any basis normalizers and therefore won't
1343 # make it to this "else" branch.
1344 XS
= a
[0].parent().gens()
1345 subs_dict
= { XS
[i
]: self
._basis
_normalizers
[i
]*XS
[i
]
1346 for i
in range(len(XS
)) }
1347 return tuple( a_i
.subs(subs_dict
) for a_i
in a
)
1353 Embed the matrix ``M`` into a space of real matrices.
1355 The matrix ``M`` can have entries in any field at the moment:
1356 the real numbers, complex numbers, or quaternions. And although
1357 they are not a field, we can probably support octonions at some
1358 point, too. This function returns a real matrix that "acts like"
1359 the original with respect to matrix multiplication; i.e.
1361 real_embed(M*N) = real_embed(M)*real_embed(N)
1364 raise NotImplementedError
1368 def real_unembed(M
):
1370 The inverse of :meth:`real_embed`.
1372 raise NotImplementedError
1375 def matrix_inner_product(cls
,X
,Y
):
1376 Xu
= cls
.real_unembed(X
)
1377 Yu
= cls
.real_unembed(Y
)
1378 tr
= (Xu
*Yu
).trace()
1381 # Works in QQ, AA, RDF, et cetera.
1383 except AttributeError:
1384 # A quaternion doesn't have a real() method, but does
1385 # have coefficient_tuple() method that returns the
1386 # coefficients of 1, i, j, and k -- in that order.
1387 return tr
.coefficient_tuple()[0]
1390 class RealMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1394 The identity function, for embedding real matrices into real
1400 def real_unembed(M
):
1402 The identity function, for unembedding real matrices from real
1408 class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra
,
1409 ConcreteEuclideanJordanAlgebra
):
1411 The rank-n simple EJA consisting of real symmetric n-by-n
1412 matrices, the usual symmetric Jordan product, and the trace inner
1413 product. It has dimension `(n^2 + n)/2` over the reals.
1417 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1421 sage: J = RealSymmetricEJA(2)
1422 sage: e0, e1, e2 = J.gens()
1430 In theory, our "field" can be any subfield of the reals::
1432 sage: RealSymmetricEJA(2, RDF)
1433 Euclidean Jordan algebra of dimension 3 over Real Double Field
1434 sage: RealSymmetricEJA(2, RR)
1435 Euclidean Jordan algebra of dimension 3 over Real Field with
1436 53 bits of precision
1440 The dimension of this algebra is `(n^2 + n) / 2`::
1442 sage: set_random_seed()
1443 sage: n_max = RealSymmetricEJA._max_random_instance_size()
1444 sage: n = ZZ.random_element(1, n_max)
1445 sage: J = RealSymmetricEJA(n)
1446 sage: J.dimension() == (n^2 + n)/2
1449 The Jordan multiplication is what we think it is::
1451 sage: set_random_seed()
1452 sage: J = RealSymmetricEJA.random_instance()
1453 sage: x,y = J.random_elements(2)
1454 sage: actual = (x*y).to_matrix()
1455 sage: X = x.to_matrix()
1456 sage: Y = y.to_matrix()
1457 sage: expected = (X*Y + Y*X)/2
1458 sage: actual == expected
1460 sage: J(expected) == x*y
1463 We can change the generator prefix::
1465 sage: RealSymmetricEJA(3, prefix='q').gens()
1466 (q0, q1, q2, q3, q4, q5)
1468 We can construct the (trivial) algebra of rank zero::
1470 sage: RealSymmetricEJA(0)
1471 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1475 def _denormalized_basis(cls
, n
, field
):
1477 Return a basis for the space of real symmetric n-by-n matrices.
1481 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1485 sage: set_random_seed()
1486 sage: n = ZZ.random_element(1,5)
1487 sage: B = RealSymmetricEJA._denormalized_basis(n,QQ)
1488 sage: all( M.is_symmetric() for M in B)
1492 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1496 for j
in range(i
+1):
1497 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
1501 Sij
= Eij
+ Eij
.transpose()
1507 def _max_random_instance_size():
1508 return 4 # Dimension 10
1511 def random_instance(cls
, field
=AA
, **kwargs
):
1513 Return a random instance of this type of algebra.
1515 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1516 return cls(n
, field
, **kwargs
)
1518 def __init__(self
, n
, field
=AA
, **kwargs
):
1519 basis
= self
._denormalized
_basis
(n
, field
)
1520 super(RealSymmetricEJA
, self
).__init
__(field
,
1524 self
.rank
.set_cache(n
)
1527 class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1531 Embed the n-by-n complex matrix ``M`` into the space of real
1532 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1533 bi` to the block matrix ``[[a,b],[-b,a]]``.
1537 sage: from mjo.eja.eja_algebra import \
1538 ....: ComplexMatrixEuclideanJordanAlgebra
1542 sage: F = QuadraticField(-1, 'I')
1543 sage: x1 = F(4 - 2*i)
1544 sage: x2 = F(1 + 2*i)
1547 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1548 sage: ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1557 Embedding is a homomorphism (isomorphism, in fact)::
1559 sage: set_random_seed()
1560 sage: n = ZZ.random_element(3)
1561 sage: F = QuadraticField(-1, 'I')
1562 sage: X = random_matrix(F, n)
1563 sage: Y = random_matrix(F, n)
1564 sage: Xe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X)
1565 sage: Ye = ComplexMatrixEuclideanJordanAlgebra.real_embed(Y)
1566 sage: XYe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1573 raise ValueError("the matrix 'M' must be square")
1575 # We don't need any adjoined elements...
1576 field
= M
.base_ring().base_ring()
1580 a
= z
.list()[0] # real part, I guess
1581 b
= z
.list()[1] # imag part, I guess
1582 blocks
.append(matrix(field
, 2, [[a
,b
],[-b
,a
]]))
1584 return matrix
.block(field
, n
, blocks
)
1588 def real_unembed(M
):
1590 The inverse of _embed_complex_matrix().
1594 sage: from mjo.eja.eja_algebra import \
1595 ....: ComplexMatrixEuclideanJordanAlgebra
1599 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1600 ....: [-2, 1, -4, 3],
1601 ....: [ 9, 10, 11, 12],
1602 ....: [-10, 9, -12, 11] ])
1603 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(A)
1605 [ 10*I + 9 12*I + 11]
1609 Unembedding is the inverse of embedding::
1611 sage: set_random_seed()
1612 sage: F = QuadraticField(-1, 'I')
1613 sage: M = random_matrix(F, 3)
1614 sage: Me = ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1615 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1621 raise ValueError("the matrix 'M' must be square")
1622 if not n
.mod(2).is_zero():
1623 raise ValueError("the matrix 'M' must be a complex embedding")
1625 # If "M" was normalized, its base ring might have roots
1626 # adjoined and they can stick around after unembedding.
1627 field
= M
.base_ring()
1628 R
= PolynomialRing(field
, 'z')
1631 # Sage doesn't know how to embed AA into QQbar, i.e. how
1632 # to adjoin sqrt(-1) to AA.
1635 F
= field
.extension(z
**2 + 1, 'I', embedding
=CLF(-1).sqrt())
1638 # Go top-left to bottom-right (reading order), converting every
1639 # 2-by-2 block we see to a single complex element.
1641 for k
in range(n
/2):
1642 for j
in range(n
/2):
1643 submat
= M
[2*k
:2*k
+2,2*j
:2*j
+2]
1644 if submat
[0,0] != submat
[1,1]:
1645 raise ValueError('bad on-diagonal submatrix')
1646 if submat
[0,1] != -submat
[1,0]:
1647 raise ValueError('bad off-diagonal submatrix')
1648 z
= submat
[0,0] + submat
[0,1]*i
1651 return matrix(F
, n
/2, elements
)
1655 def matrix_inner_product(cls
,X
,Y
):
1657 Compute a matrix inner product in this algebra directly from
1662 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1666 This gives the same answer as the slow, default method implemented
1667 in :class:`MatrixEuclideanJordanAlgebra`::
1669 sage: set_random_seed()
1670 sage: J = ComplexHermitianEJA.random_instance()
1671 sage: x,y = J.random_elements(2)
1672 sage: Xe = x.to_matrix()
1673 sage: Ye = y.to_matrix()
1674 sage: X = ComplexHermitianEJA.real_unembed(Xe)
1675 sage: Y = ComplexHermitianEJA.real_unembed(Ye)
1676 sage: expected = (X*Y).trace().real()
1677 sage: actual = ComplexHermitianEJA.matrix_inner_product(Xe,Ye)
1678 sage: actual == expected
1682 return RealMatrixEuclideanJordanAlgebra
.matrix_inner_product(X
,Y
)/2
1685 class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra
,
1686 ConcreteEuclideanJordanAlgebra
):
1688 The rank-n simple EJA consisting of complex Hermitian n-by-n
1689 matrices over the real numbers, the usual symmetric Jordan product,
1690 and the real-part-of-trace inner product. It has dimension `n^2` over
1695 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1699 In theory, our "field" can be any subfield of the reals::
1701 sage: ComplexHermitianEJA(2, RDF)
1702 Euclidean Jordan algebra of dimension 4 over Real Double Field
1703 sage: ComplexHermitianEJA(2, RR)
1704 Euclidean Jordan algebra of dimension 4 over Real Field with
1705 53 bits of precision
1709 The dimension of this algebra is `n^2`::
1711 sage: set_random_seed()
1712 sage: n_max = ComplexHermitianEJA._max_random_instance_size()
1713 sage: n = ZZ.random_element(1, n_max)
1714 sage: J = ComplexHermitianEJA(n)
1715 sage: J.dimension() == n^2
1718 The Jordan multiplication is what we think it is::
1720 sage: set_random_seed()
1721 sage: J = ComplexHermitianEJA.random_instance()
1722 sage: x,y = J.random_elements(2)
1723 sage: actual = (x*y).to_matrix()
1724 sage: X = x.to_matrix()
1725 sage: Y = y.to_matrix()
1726 sage: expected = (X*Y + Y*X)/2
1727 sage: actual == expected
1729 sage: J(expected) == x*y
1732 We can change the generator prefix::
1734 sage: ComplexHermitianEJA(2, prefix='z').gens()
1737 We can construct the (trivial) algebra of rank zero::
1739 sage: ComplexHermitianEJA(0)
1740 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1745 def _denormalized_basis(cls
, n
, field
):
1747 Returns a basis for the space of complex Hermitian n-by-n matrices.
1749 Why do we embed these? Basically, because all of numerical linear
1750 algebra assumes that you're working with vectors consisting of `n`
1751 entries from a field and scalars from the same field. There's no way
1752 to tell SageMath that (for example) the vectors contain complex
1753 numbers, while the scalar field is real.
1757 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1761 sage: set_random_seed()
1762 sage: n = ZZ.random_element(1,5)
1763 sage: field = QuadraticField(2, 'sqrt2')
1764 sage: B = ComplexHermitianEJA._denormalized_basis(n, field)
1765 sage: all( M.is_symmetric() for M in B)
1769 R
= PolynomialRing(field
, 'z')
1771 F
= field
.extension(z
**2 + 1, 'I')
1774 # This is like the symmetric case, but we need to be careful:
1776 # * We want conjugate-symmetry, not just symmetry.
1777 # * The diagonal will (as a result) be real.
1781 for j
in range(i
+1):
1782 Eij
= matrix(F
, n
, lambda k
,l
: k
==i
and l
==j
)
1784 Sij
= cls
.real_embed(Eij
)
1787 # The second one has a minus because it's conjugated.
1788 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
1790 Sij_imag
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
1793 # Since we embedded these, we can drop back to the "field" that we
1794 # started with instead of the complex extension "F".
1795 return ( s
.change_ring(field
) for s
in S
)
1798 def __init__(self
, n
, field
=AA
, **kwargs
):
1799 basis
= self
._denormalized
_basis
(n
,field
)
1800 super(ComplexHermitianEJA
,self
).__init
__(field
,
1804 self
.rank
.set_cache(n
)
1807 def _max_random_instance_size():
1808 return 3 # Dimension 9
1811 def random_instance(cls
, field
=AA
, **kwargs
):
1813 Return a random instance of this type of algebra.
1815 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1816 return cls(n
, field
, **kwargs
)
1818 class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1822 Embed the n-by-n quaternion matrix ``M`` into the space of real
1823 matrices of size 4n-by-4n by first sending each quaternion entry `z
1824 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
1825 c+di],[-c + di, a-bi]]`, and then embedding those into a real
1830 sage: from mjo.eja.eja_algebra import \
1831 ....: QuaternionMatrixEuclideanJordanAlgebra
1835 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1836 sage: i,j,k = Q.gens()
1837 sage: x = 1 + 2*i + 3*j + 4*k
1838 sage: M = matrix(Q, 1, [[x]])
1839 sage: QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1845 Embedding is a homomorphism (isomorphism, in fact)::
1847 sage: set_random_seed()
1848 sage: n = ZZ.random_element(2)
1849 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1850 sage: X = random_matrix(Q, n)
1851 sage: Y = random_matrix(Q, n)
1852 sage: Xe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X)
1853 sage: Ye = QuaternionMatrixEuclideanJordanAlgebra.real_embed(Y)
1854 sage: XYe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1859 quaternions
= M
.base_ring()
1862 raise ValueError("the matrix 'M' must be square")
1864 F
= QuadraticField(-1, 'I')
1869 t
= z
.coefficient_tuple()
1874 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
1875 [-c
+ d
*i
, a
- b
*i
]])
1876 realM
= ComplexMatrixEuclideanJordanAlgebra
.real_embed(cplxM
)
1877 blocks
.append(realM
)
1879 # We should have real entries by now, so use the realest field
1880 # we've got for the return value.
1881 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
1886 def real_unembed(M
):
1888 The inverse of _embed_quaternion_matrix().
1892 sage: from mjo.eja.eja_algebra import \
1893 ....: QuaternionMatrixEuclideanJordanAlgebra
1897 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
1898 ....: [-2, 1, -4, 3],
1899 ....: [-3, 4, 1, -2],
1900 ....: [-4, -3, 2, 1]])
1901 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(M)
1902 [1 + 2*i + 3*j + 4*k]
1906 Unembedding is the inverse of embedding::
1908 sage: set_random_seed()
1909 sage: Q = QuaternionAlgebra(QQ, -1, -1)
1910 sage: M = random_matrix(Q, 3)
1911 sage: Me = QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1912 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1918 raise ValueError("the matrix 'M' must be square")
1919 if not n
.mod(4).is_zero():
1920 raise ValueError("the matrix 'M' must be a quaternion embedding")
1922 # Use the base ring of the matrix to ensure that its entries can be
1923 # multiplied by elements of the quaternion algebra.
1924 field
= M
.base_ring()
1925 Q
= QuaternionAlgebra(field
,-1,-1)
1928 # Go top-left to bottom-right (reading order), converting every
1929 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
1932 for l
in range(n
/4):
1933 for m
in range(n
/4):
1934 submat
= ComplexMatrixEuclideanJordanAlgebra
.real_unembed(
1935 M
[4*l
:4*l
+4,4*m
:4*m
+4] )
1936 if submat
[0,0] != submat
[1,1].conjugate():
1937 raise ValueError('bad on-diagonal submatrix')
1938 if submat
[0,1] != -submat
[1,0].conjugate():
1939 raise ValueError('bad off-diagonal submatrix')
1940 z
= submat
[0,0].real()
1941 z
+= submat
[0,0].imag()*i
1942 z
+= submat
[0,1].real()*j
1943 z
+= submat
[0,1].imag()*k
1946 return matrix(Q
, n
/4, elements
)
1950 def matrix_inner_product(cls
,X
,Y
):
1952 Compute a matrix inner product in this algebra directly from
1957 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1961 This gives the same answer as the slow, default method implemented
1962 in :class:`MatrixEuclideanJordanAlgebra`::
1964 sage: set_random_seed()
1965 sage: J = QuaternionHermitianEJA.random_instance()
1966 sage: x,y = J.random_elements(2)
1967 sage: Xe = x.to_matrix()
1968 sage: Ye = y.to_matrix()
1969 sage: X = QuaternionHermitianEJA.real_unembed(Xe)
1970 sage: Y = QuaternionHermitianEJA.real_unembed(Ye)
1971 sage: expected = (X*Y).trace().coefficient_tuple()[0]
1972 sage: actual = QuaternionHermitianEJA.matrix_inner_product(Xe,Ye)
1973 sage: actual == expected
1977 return RealMatrixEuclideanJordanAlgebra
.matrix_inner_product(X
,Y
)/4
1980 class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra
,
1981 ConcreteEuclideanJordanAlgebra
):
1983 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
1984 matrices, the usual symmetric Jordan product, and the
1985 real-part-of-trace inner product. It has dimension `2n^2 - n` over
1990 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1994 In theory, our "field" can be any subfield of the reals::
1996 sage: QuaternionHermitianEJA(2, RDF)
1997 Euclidean Jordan algebra of dimension 6 over Real Double Field
1998 sage: QuaternionHermitianEJA(2, RR)
1999 Euclidean Jordan algebra of dimension 6 over Real Field with
2000 53 bits of precision
2004 The dimension of this algebra is `2*n^2 - n`::
2006 sage: set_random_seed()
2007 sage: n_max = QuaternionHermitianEJA._max_random_instance_size()
2008 sage: n = ZZ.random_element(1, n_max)
2009 sage: J = QuaternionHermitianEJA(n)
2010 sage: J.dimension() == 2*(n^2) - n
2013 The Jordan multiplication is what we think it is::
2015 sage: set_random_seed()
2016 sage: J = QuaternionHermitianEJA.random_instance()
2017 sage: x,y = J.random_elements(2)
2018 sage: actual = (x*y).to_matrix()
2019 sage: X = x.to_matrix()
2020 sage: Y = y.to_matrix()
2021 sage: expected = (X*Y + Y*X)/2
2022 sage: actual == expected
2024 sage: J(expected) == x*y
2027 We can change the generator prefix::
2029 sage: QuaternionHermitianEJA(2, prefix='a').gens()
2030 (a0, a1, a2, a3, a4, a5)
2032 We can construct the (trivial) algebra of rank zero::
2034 sage: QuaternionHermitianEJA(0)
2035 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2039 def _denormalized_basis(cls
, n
, field
):
2041 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
2043 Why do we embed these? Basically, because all of numerical
2044 linear algebra assumes that you're working with vectors consisting
2045 of `n` entries from a field and scalars from the same field. There's
2046 no way to tell SageMath that (for example) the vectors contain
2047 complex numbers, while the scalar field is real.
2051 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2055 sage: set_random_seed()
2056 sage: n = ZZ.random_element(1,5)
2057 sage: B = QuaternionHermitianEJA._denormalized_basis(n,QQ)
2058 sage: all( M.is_symmetric() for M in B )
2062 Q
= QuaternionAlgebra(QQ
,-1,-1)
2065 # This is like the symmetric case, but we need to be careful:
2067 # * We want conjugate-symmetry, not just symmetry.
2068 # * The diagonal will (as a result) be real.
2072 for j
in range(i
+1):
2073 Eij
= matrix(Q
, n
, lambda k
,l
: k
==i
and l
==j
)
2075 Sij
= cls
.real_embed(Eij
)
2078 # The second, third, and fourth ones have a minus
2079 # because they're conjugated.
2080 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
2082 Sij_I
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
2084 Sij_J
= cls
.real_embed(J
*Eij
- J
*Eij
.transpose())
2086 Sij_K
= cls
.real_embed(K
*Eij
- K
*Eij
.transpose())
2089 # Since we embedded these, we can drop back to the "field" that we
2090 # started with instead of the quaternion algebra "Q".
2091 return ( s
.change_ring(field
) for s
in S
)
2094 def __init__(self
, n
, field
=AA
, **kwargs
):
2095 basis
= self
._denormalized
_basis
(n
,field
)
2096 super(QuaternionHermitianEJA
,self
).__init
__(field
,
2100 self
.rank
.set_cache(n
)
2103 def _max_random_instance_size():
2105 The maximum rank of a random QuaternionHermitianEJA.
2107 return 2 # Dimension 6
2110 def random_instance(cls
, field
=AA
, **kwargs
):
2112 Return a random instance of this type of algebra.
2114 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2115 return cls(n
, field
, **kwargs
)
2118 class HadamardEJA(RationalBasisEuclideanJordanAlgebraNg
,
2119 ConcreteEuclideanJordanAlgebra
):
2121 Return the Euclidean Jordan Algebra corresponding to the set
2122 `R^n` under the Hadamard product.
2124 Note: this is nothing more than the Cartesian product of ``n``
2125 copies of the spin algebra. Once Cartesian product algebras
2126 are implemented, this can go.
2130 sage: from mjo.eja.eja_algebra import HadamardEJA
2134 This multiplication table can be verified by hand::
2136 sage: J = HadamardEJA(3)
2137 sage: e0,e1,e2 = J.gens()
2153 We can change the generator prefix::
2155 sage: HadamardEJA(3, prefix='r').gens()
2159 def __init__(self
, n
, field
=AA
, **kwargs
):
2160 V
= VectorSpace(field
, n
)
2163 def jordan_product(x
,y
):
2164 return V([ xi
*yi
for (xi
,yi
) in zip(x
,y
) ])
2165 def inner_product(x
,y
):
2166 return x
.inner_product(y
)
2168 super(HadamardEJA
, self
).__init
__(field
,
2173 self
.rank
.set_cache(n
)
2176 self
.one
.set_cache( self
.zero() )
2178 self
.one
.set_cache( sum(self
.gens()) )
2181 def _max_random_instance_size():
2183 The maximum dimension of a random HadamardEJA.
2188 def random_instance(cls
, field
=AA
, **kwargs
):
2190 Return a random instance of this type of algebra.
2192 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2193 return cls(n
, field
, **kwargs
)
2196 class BilinearFormEJA(RationalBasisEuclideanJordanAlgebra
,
2197 ConcreteEuclideanJordanAlgebra
):
2199 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2200 with the half-trace inner product and jordan product ``x*y =
2201 (<Bx,y>,y_bar>, x0*y_bar + y0*x_bar)`` where `B = 1 \times B22` is
2202 a symmetric positive-definite "bilinear form" matrix. Its
2203 dimension is the size of `B`, and it has rank two in dimensions
2204 larger than two. It reduces to the ``JordanSpinEJA`` when `B` is
2205 the identity matrix of order ``n``.
2207 We insist that the one-by-one upper-left identity block of `B` be
2208 passed in as well so that we can be passed a matrix of size zero
2209 to construct a trivial algebra.
2213 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
2214 ....: JordanSpinEJA)
2218 When no bilinear form is specified, the identity matrix is used,
2219 and the resulting algebra is the Jordan spin algebra::
2221 sage: B = matrix.identity(AA,3)
2222 sage: J0 = BilinearFormEJA(B)
2223 sage: J1 = JordanSpinEJA(3)
2224 sage: J0.multiplication_table() == J0.multiplication_table()
2227 An error is raised if the matrix `B` does not correspond to a
2228 positive-definite bilinear form::
2230 sage: B = matrix.random(QQ,2,3)
2231 sage: J = BilinearFormEJA(B)
2232 Traceback (most recent call last):
2234 ValueError: bilinear form is not positive-definite
2235 sage: B = matrix.zero(QQ,3)
2236 sage: J = BilinearFormEJA(B)
2237 Traceback (most recent call last):
2239 ValueError: bilinear form is not positive-definite
2243 We can create a zero-dimensional algebra::
2245 sage: B = matrix.identity(AA,0)
2246 sage: J = BilinearFormEJA(B)
2250 We can check the multiplication condition given in the Jordan, von
2251 Neumann, and Wigner paper (and also discussed on my "On the
2252 symmetry..." paper). Note that this relies heavily on the standard
2253 choice of basis, as does anything utilizing the bilinear form matrix::
2255 sage: set_random_seed()
2256 sage: n = ZZ.random_element(5)
2257 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2258 sage: B11 = matrix.identity(QQ,1)
2259 sage: B22 = M.transpose()*M
2260 sage: B = block_matrix(2,2,[ [B11,0 ],
2262 sage: J = BilinearFormEJA(B)
2263 sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
2264 sage: V = J.vector_space()
2265 sage: sis = [ J.from_vector(V([0] + (M.inverse()*ei).list()))
2266 ....: for ei in eis ]
2267 sage: actual = [ sis[i]*sis[j]
2268 ....: for i in range(n-1)
2269 ....: for j in range(n-1) ]
2270 sage: expected = [ J.one() if i == j else J.zero()
2271 ....: for i in range(n-1)
2272 ....: for j in range(n-1) ]
2273 sage: actual == expected
2276 def __init__(self
, B
, field
=AA
, **kwargs
):
2279 if not B
.is_positive_definite():
2280 raise ValueError("bilinear form is not positive-definite")
2282 V
= VectorSpace(field
, n
)
2283 mult_table
= [[V
.zero() for j
in range(n
)] for i
in range(n
)]
2292 z0
= (B
*x
).inner_product(y
)
2293 zbar
= y0
*xbar
+ x0
*ybar
2294 z
= V([z0
] + zbar
.list())
2295 mult_table
[i
][j
] = z
2297 # Inner products are real numbers and not algebra
2298 # elements, so once we turn the algebra element
2299 # into a vector in inner_product(), we never go
2300 # back. As a result -- contrary to what we do with
2301 # self._multiplication_table -- we store the inner
2302 # product table as a plain old matrix and not as
2303 # an algebra operator.
2305 self
._inner
_product
_matrix
= ip_table
2307 super(BilinearFormEJA
, self
).__init
__(field
,
2312 # The rank of this algebra is two, unless we're in a
2313 # one-dimensional ambient space (because the rank is bounded
2314 # by the ambient dimension).
2315 self
.rank
.set_cache(min(n
,2))
2318 self
.one
.set_cache( self
.zero() )
2320 self
.one
.set_cache( self
.monomial(0) )
2323 def _max_random_instance_size():
2325 The maximum dimension of a random BilinearFormEJA.
2330 def random_instance(cls
, field
=AA
, **kwargs
):
2332 Return a random instance of this algebra.
2334 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2336 B
= matrix
.identity(field
, n
)
2337 return cls(B
, field
, **kwargs
)
2339 B11
= matrix
.identity(field
,1)
2340 M
= matrix
.random(field
, n
-1)
2341 I
= matrix
.identity(field
, n
-1)
2342 alpha
= field
.zero()
2343 while alpha
.is_zero():
2344 alpha
= field
.random_element().abs()
2345 B22
= M
.transpose()*M
+ alpha
*I
2347 from sage
.matrix
.special
import block_matrix
2348 B
= block_matrix(2,2, [ [B11
, ZZ(0) ],
2351 return cls(B
, field
, **kwargs
)
2354 class JordanSpinEJA(BilinearFormEJA
):
2356 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2357 with the usual inner product and jordan product ``x*y =
2358 (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2363 sage: from mjo.eja.eja_algebra import JordanSpinEJA
2367 This multiplication table can be verified by hand::
2369 sage: J = JordanSpinEJA(4)
2370 sage: e0,e1,e2,e3 = J.gens()
2386 We can change the generator prefix::
2388 sage: JordanSpinEJA(2, prefix='B').gens()
2393 Ensure that we have the usual inner product on `R^n`::
2395 sage: set_random_seed()
2396 sage: J = JordanSpinEJA.random_instance()
2397 sage: x,y = J.random_elements(2)
2398 sage: actual = x.inner_product(y)
2399 sage: expected = x.to_vector().inner_product(y.to_vector())
2400 sage: actual == expected
2404 def __init__(self
, n
, field
=AA
, **kwargs
):
2405 # This is a special case of the BilinearFormEJA with the identity
2406 # matrix as its bilinear form.
2407 B
= matrix
.identity(field
, n
)
2408 super(JordanSpinEJA
, self
).__init
__(B
, field
, **kwargs
)
2411 def _max_random_instance_size():
2413 The maximum dimension of a random JordanSpinEJA.
2418 def random_instance(cls
, field
=AA
, **kwargs
):
2420 Return a random instance of this type of algebra.
2422 Needed here to override the implementation for ``BilinearFormEJA``.
2424 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2425 return cls(n
, field
, **kwargs
)
2428 class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra
,
2429 ConcreteEuclideanJordanAlgebra
):
2431 The trivial Euclidean Jordan algebra consisting of only a zero element.
2435 sage: from mjo.eja.eja_algebra import TrivialEJA
2439 sage: J = TrivialEJA()
2446 sage: 7*J.one()*12*J.one()
2448 sage: J.one().inner_product(J.one())
2450 sage: J.one().norm()
2452 sage: J.one().subalgebra_generated_by()
2453 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2458 def __init__(self
, field
=AA
, **kwargs
):
2460 self
._inner
_product
_matrix
= matrix(field
,0)
2461 super(TrivialEJA
, self
).__init
__(field
,
2465 # The rank is zero using my definition, namely the dimension of the
2466 # largest subalgebra generated by any element.
2467 self
.rank
.set_cache(0)
2468 self
.one
.set_cache( self
.zero() )
2471 def random_instance(cls
, field
=AA
, **kwargs
):
2472 # We don't take a "size" argument so the superclass method is
2473 # inappropriate for us.
2474 return cls(field
, **kwargs
)
2476 class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra
):
2478 The external (orthogonal) direct sum of two other Euclidean Jordan
2479 algebras. Essentially the Cartesian product of its two factors.
2480 Every Euclidean Jordan algebra decomposes into an orthogonal
2481 direct sum of simple Euclidean Jordan algebras, so no generality
2482 is lost by providing only this construction.
2486 sage: from mjo.eja.eja_algebra import (random_eja,
2488 ....: RealSymmetricEJA,
2493 sage: J1 = HadamardEJA(2)
2494 sage: J2 = RealSymmetricEJA(3)
2495 sage: J = DirectSumEJA(J1,J2)
2503 The external direct sum construction is only valid when the two factors
2504 have the same base ring; an error is raised otherwise::
2506 sage: set_random_seed()
2507 sage: J1 = random_eja(AA)
2508 sage: J2 = random_eja(QQ)
2509 sage: J = DirectSumEJA(J1,J2)
2510 Traceback (most recent call last):
2512 ValueError: algebras must share the same base field
2515 def __init__(self
, J1
, J2
, **kwargs
):
2516 if J1
.base_ring() != J2
.base_ring():
2517 raise ValueError("algebras must share the same base field")
2518 field
= J1
.base_ring()
2520 self
._factors
= (J1
, J2
)
2524 V
= VectorSpace(field
, n
)
2525 mult_table
= [ [ V
.zero() for j
in range(n
) ]
2529 p
= (J1
.monomial(i
)*J1
.monomial(j
)).to_vector()
2530 mult_table
[i
][j
] = V(p
.list() + [field
.zero()]*n2
)
2534 p
= (J2
.monomial(i
)*J2
.monomial(j
)).to_vector()
2535 mult_table
[n1
+i
][n1
+j
] = V([field
.zero()]*n1
+ p
.list())
2537 super(DirectSumEJA
, self
).__init
__(field
,
2541 self
.rank
.set_cache(J1
.rank() + J2
.rank())
2546 Return the pair of this algebra's factors.
2550 sage: from mjo.eja.eja_algebra import (HadamardEJA,
2551 ....: JordanSpinEJA,
2556 sage: J1 = HadamardEJA(2,QQ)
2557 sage: J2 = JordanSpinEJA(3,QQ)
2558 sage: J = DirectSumEJA(J1,J2)
2560 (Euclidean Jordan algebra of dimension 2 over Rational Field,
2561 Euclidean Jordan algebra of dimension 3 over Rational Field)
2564 return self
._factors
2566 def projections(self
):
2568 Return a pair of projections onto this algebra's factors.
2572 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
2573 ....: ComplexHermitianEJA,
2578 sage: J1 = JordanSpinEJA(2)
2579 sage: J2 = ComplexHermitianEJA(2)
2580 sage: J = DirectSumEJA(J1,J2)
2581 sage: (pi_left, pi_right) = J.projections()
2582 sage: J.one().to_vector()
2584 sage: pi_left(J.one()).to_vector()
2586 sage: pi_right(J.one()).to_vector()
2590 (J1
,J2
) = self
.factors()
2593 V_basis
= self
.vector_space().basis()
2594 # Need to specify the dimensions explicitly so that we don't
2595 # wind up with a zero-by-zero matrix when we want e.g. a
2596 # zero-by-two matrix (important for composing things).
2597 P1
= matrix(self
.base_ring(), m
, m
+n
, V_basis
[:m
])
2598 P2
= matrix(self
.base_ring(), n
, m
+n
, V_basis
[m
:])
2599 pi_left
= FiniteDimensionalEuclideanJordanAlgebraOperator(self
,J1
,P1
)
2600 pi_right
= FiniteDimensionalEuclideanJordanAlgebraOperator(self
,J2
,P2
)
2601 return (pi_left
, pi_right
)
2603 def inclusions(self
):
2605 Return the pair of inclusion maps from our factors into us.
2609 sage: from mjo.eja.eja_algebra import (random_eja,
2610 ....: JordanSpinEJA,
2611 ....: RealSymmetricEJA,
2616 sage: J1 = JordanSpinEJA(3)
2617 sage: J2 = RealSymmetricEJA(2)
2618 sage: J = DirectSumEJA(J1,J2)
2619 sage: (iota_left, iota_right) = J.inclusions()
2620 sage: iota_left(J1.zero()) == J.zero()
2622 sage: iota_right(J2.zero()) == J.zero()
2624 sage: J1.one().to_vector()
2626 sage: iota_left(J1.one()).to_vector()
2628 sage: J2.one().to_vector()
2630 sage: iota_right(J2.one()).to_vector()
2632 sage: J.one().to_vector()
2637 Composing a projection with the corresponding inclusion should
2638 produce the identity map, and mismatching them should produce
2641 sage: set_random_seed()
2642 sage: J1 = random_eja()
2643 sage: J2 = random_eja()
2644 sage: J = DirectSumEJA(J1,J2)
2645 sage: (iota_left, iota_right) = J.inclusions()
2646 sage: (pi_left, pi_right) = J.projections()
2647 sage: pi_left*iota_left == J1.one().operator()
2649 sage: pi_right*iota_right == J2.one().operator()
2651 sage: (pi_left*iota_right).is_zero()
2653 sage: (pi_right*iota_left).is_zero()
2657 (J1
,J2
) = self
.factors()
2660 V_basis
= self
.vector_space().basis()
2661 # Need to specify the dimensions explicitly so that we don't
2662 # wind up with a zero-by-zero matrix when we want e.g. a
2663 # two-by-zero matrix (important for composing things).
2664 I1
= matrix
.column(self
.base_ring(), m
, m
+n
, V_basis
[:m
])
2665 I2
= matrix
.column(self
.base_ring(), n
, m
+n
, V_basis
[m
:])
2666 iota_left
= FiniteDimensionalEuclideanJordanAlgebraOperator(J1
,self
,I1
)
2667 iota_right
= FiniteDimensionalEuclideanJordanAlgebraOperator(J2
,self
,I2
)
2668 return (iota_left
, iota_right
)
2670 def inner_product(self
, x
, y
):
2672 The standard Cartesian inner-product.
2674 We project ``x`` and ``y`` onto our factors, and add up the
2675 inner-products from the subalgebras.
2680 sage: from mjo.eja.eja_algebra import (HadamardEJA,
2681 ....: QuaternionHermitianEJA,
2686 sage: J1 = HadamardEJA(3,QQ)
2687 sage: J2 = QuaternionHermitianEJA(2,QQ,normalize_basis=False)
2688 sage: J = DirectSumEJA(J1,J2)
2693 sage: x1.inner_product(x2)
2695 sage: y1.inner_product(y2)
2697 sage: J.one().inner_product(J.one())
2701 (pi_left
, pi_right
) = self
.projections()
2707 return (x1
.inner_product(y1
) + x2
.inner_product(y2
))
2711 random_eja
= ConcreteEuclideanJordanAlgebra
.random_instance