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
.matrix
.matrix
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 self
._deorthonormalization
_matrix
= matrix
.identity(field
,n
)
1059 A
= matrix(field
, vector_basis
)
1060 # uh oh, this is only the "usual" inner product
1061 Q
,R
= A
.gram_schmidt(orthonormal
=True)
1062 self
._deorthonormalization
_matrix
= R
.inverse().transpose()
1063 vector_basis
= Q
.rows()
1064 W
= V
.span_of_basis( vector_basis
)
1065 if basis_is_matrices
:
1066 from mjo
.eja
.eja_utils
import _vec2mat
1067 basis
= tuple( map(_vec2mat
,vector_basis
) )
1069 mult_table
= [ [0 for i
in range(n
)] for j
in range(n
) ]
1070 ip_table
= [ [0 for i
in range(n
)] for j
in range(n
) ]
1073 for j
in range(i
+1):
1074 # do another mat2vec because the multiplication
1075 # table is in terms of vectors
1076 elt
= _mat2vec(jordan_product(basis
[i
],basis
[j
]))
1077 elt
= W
.coordinate_vector(elt
)
1078 mult_table
[i
][j
] = elt
1079 mult_table
[j
][i
] = elt
1080 ip
= inner_product(basis
[i
],basis
[j
])
1084 self
._inner
_product
_matrix
= matrix(field
,ip_table
)
1086 if basis_is_matrices
:
1090 super().__init
__(field
,
1098 class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra
):
1100 Algebras whose basis consists of vectors with rational
1101 entries. Equivalently, algebras whose multiplication tables
1102 contain only rational coefficients.
1104 When an EJA has a basis that can be made rational, we can speed up
1105 the computation of its characteristic polynomial by doing it over
1106 ``QQ``. All of the named EJA constructors that we provide fall
1110 def _charpoly_coefficients(self
):
1112 Override the parent method with something that tries to compute
1113 over a faster (non-extension) field.
1117 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1121 The base ring of the resulting polynomial coefficients is what
1122 it should be, and not the rationals (unless the algebra was
1123 already over the rationals)::
1125 sage: J = JordanSpinEJA(3)
1126 sage: J._charpoly_coefficients()
1127 (X1^2 - X2^2 - X3^2, -2*X1)
1128 sage: a0 = J._charpoly_coefficients()[0]
1130 Algebraic Real Field
1131 sage: a0.base_ring()
1132 Algebraic Real Field
1135 if self
.base_ring() is QQ
:
1136 # There's no need to construct *another* algebra over the
1137 # rationals if this one is already over the rationals.
1138 superclass
= super(RationalBasisEuclideanJordanAlgebra
, self
)
1139 return superclass
._charpoly
_coefficients
()
1142 tuple(map(lambda x
: x
.to_vector(), ls
))
1143 for ls
in self
._multiplication
_table
1146 # Do the computation over the rationals. The answer will be
1147 # the same, because our basis coordinates are (essentially)
1149 J
= FiniteDimensionalEuclideanJordanAlgebra(QQ
,
1153 a
= J
._charpoly
_coefficients
()
1154 return tuple(map(lambda x
: x
.change_ring(self
.base_ring()), a
))
1157 class ConcreteEuclideanJordanAlgebra
:
1159 A class for the Euclidean Jordan algebras that we know by name.
1161 These are the Jordan algebras whose basis, multiplication table,
1162 rank, and so on are known a priori. More to the point, they are
1163 the Euclidean Jordan algebras for which we are able to conjure up
1164 a "random instance."
1168 sage: from mjo.eja.eja_algebra import ConcreteEuclideanJordanAlgebra
1172 Our basis is normalized with respect to the algebra's inner
1173 product, unless we specify otherwise::
1175 sage: set_random_seed()
1176 sage: J = ConcreteEuclideanJordanAlgebra.random_instance()
1177 sage: all( b.norm() == 1 for b in J.gens() )
1180 Since our basis is orthonormal with respect to the algebra's inner
1181 product, and since we know that this algebra is an EJA, any
1182 left-multiplication operator's matrix will be symmetric because
1183 natural->EJA basis representation is an isometry and within the
1184 EJA the operator is self-adjoint by the Jordan axiom::
1186 sage: set_random_seed()
1187 sage: J = ConcreteEuclideanJordanAlgebra.random_instance()
1188 sage: x = J.random_element()
1189 sage: x.operator().is_self_adjoint()
1194 def _max_random_instance_size():
1196 Return an integer "size" that is an upper bound on the size of
1197 this algebra when it is used in a random test
1198 case. Unfortunately, the term "size" is ambiguous -- when
1199 dealing with `R^n` under either the Hadamard or Jordan spin
1200 product, the "size" refers to the dimension `n`. When dealing
1201 with a matrix algebra (real symmetric or complex/quaternion
1202 Hermitian), it refers to the size of the matrix, which is far
1203 less than the dimension of the underlying vector space.
1205 This method must be implemented in each subclass.
1207 raise NotImplementedError
1210 def random_instance(cls
, field
=AA
, **kwargs
):
1212 Return a random instance of this type of algebra.
1214 This method should be implemented in each subclass.
1216 from sage
.misc
.prandom
import choice
1217 eja_class
= choice(cls
.__subclasses
__())
1218 return eja_class
.random_instance(field
)
1221 class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra
):
1223 def __init__(self
, field
, basis
, normalize_basis
=True, **kwargs
):
1225 Compared to the superclass constructor, we take a basis instead of
1226 a multiplication table because the latter can be computed in terms
1227 of the former when the product is known (like it is here).
1229 # Used in this class's fast _charpoly_coefficients() override.
1230 self
._basis
_normalizers
= None
1232 # We're going to loop through this a few times, so now's a good
1233 # time to ensure that it isn't a generator expression.
1234 basis
= tuple(basis
)
1236 algebra_dim
= len(basis
)
1237 degree
= 0 # size of the matrices
1239 degree
= basis
[0].nrows()
1241 if algebra_dim
> 1 and normalize_basis
:
1242 # We'll need sqrt(2) to normalize the basis, and this
1243 # winds up in the multiplication table, so the whole
1244 # algebra needs to be over the field extension.
1245 R
= PolynomialRing(field
, 'z')
1248 if p
.is_irreducible():
1249 field
= field
.extension(p
, 'sqrt2', embedding
=RLF(2).sqrt())
1250 basis
= tuple( s
.change_ring(field
) for s
in basis
)
1251 self
._basis
_normalizers
= tuple(
1252 ~
(self
.matrix_inner_product(s
,s
).sqrt()) for s
in basis
)
1253 basis
= tuple(s
*c
for (s
,c
) in zip(basis
,self
._basis
_normalizers
))
1255 # Now compute the multiplication and inner product tables.
1256 # We have to do this *after* normalizing the basis, because
1257 # scaling affects the answers.
1258 V
= VectorSpace(field
, degree
**2)
1259 W
= V
.span_of_basis( _mat2vec(s
) for s
in basis
)
1260 mult_table
= [[W
.zero() for j
in range(algebra_dim
)]
1261 for i
in range(algebra_dim
)]
1262 ip_table
= [[W
.zero() for j
in range(algebra_dim
)]
1263 for i
in range(algebra_dim
)]
1264 for i
in range(algebra_dim
):
1265 for j
in range(algebra_dim
):
1266 mat_entry
= (basis
[i
]*basis
[j
] + basis
[j
]*basis
[i
])/2
1267 mult_table
[i
][j
] = W
.coordinate_vector(_mat2vec(mat_entry
))
1270 # HACK: ignore the error here if we don't need the
1271 # inner product (as is the case when we construct
1272 # a dummy QQ-algebra for fast charpoly coefficients.
1273 ip_table
[i
][j
] = self
.matrix_inner_product(basis
[i
],
1280 self
._inner
_product
_matrix
= matrix(field
,ip_table
)
1284 super(MatrixEuclideanJordanAlgebra
, self
).__init
__(field
,
1289 if algebra_dim
== 0:
1290 self
.one
.set_cache(self
.zero())
1292 n
= basis
[0].nrows()
1293 # The identity wrt (A,B) -> (AB + BA)/2 is independent of the
1294 # details of this algebra.
1295 self
.one
.set_cache(self(matrix
.identity(field
,n
)))
1299 def _charpoly_coefficients(self
):
1301 Override the parent method with something that tries to compute
1302 over a faster (non-extension) field.
1304 if self
._basis
_normalizers
is None or self
.base_ring() is QQ
:
1305 # We didn't normalize, or the basis we started with had
1306 # entries in a nice field already. Just compute the thing.
1307 return super(MatrixEuclideanJordanAlgebra
, self
)._charpoly
_coefficients
()
1309 basis
= ( (b
/n
) for (b
,n
) in zip(self
.matrix_basis(),
1310 self
._basis
_normalizers
) )
1312 # Do this over the rationals and convert back at the end.
1313 # Only works because we know the entries of the basis are
1314 # integers. The argument ``check_axioms=False`` is required
1315 # because the trace inner-product method for this
1316 # class is a stub and can't actually be checked.
1317 J
= MatrixEuclideanJordanAlgebra(QQ
,
1319 normalize_basis
=False,
1322 a
= J
._charpoly
_coefficients
()
1324 # Unfortunately, changing the basis does change the
1325 # coefficients of the characteristic polynomial, but since
1326 # these are really the coefficients of the "characteristic
1327 # polynomial of" function, everything is still nice and
1328 # unevaluated. It's therefore "obvious" how scaling the
1329 # basis affects the coordinate variables X1, X2, et
1330 # cetera. Scaling the first basis vector up by "n" adds a
1331 # factor of 1/n into every "X1" term, for example. So here
1332 # we simply undo the basis_normalizer scaling that we
1333 # performed earlier.
1335 # The a[0] access here is safe because trivial algebras
1336 # won't have any basis normalizers and therefore won't
1337 # make it to this "else" branch.
1338 XS
= a
[0].parent().gens()
1339 subs_dict
= { XS
[i
]: self
._basis
_normalizers
[i
]*XS
[i
]
1340 for i
in range(len(XS
)) }
1341 return tuple( a_i
.subs(subs_dict
) for a_i
in a
)
1347 Embed the matrix ``M`` into a space of real matrices.
1349 The matrix ``M`` can have entries in any field at the moment:
1350 the real numbers, complex numbers, or quaternions. And although
1351 they are not a field, we can probably support octonions at some
1352 point, too. This function returns a real matrix that "acts like"
1353 the original with respect to matrix multiplication; i.e.
1355 real_embed(M*N) = real_embed(M)*real_embed(N)
1358 raise NotImplementedError
1362 def real_unembed(M
):
1364 The inverse of :meth:`real_embed`.
1366 raise NotImplementedError
1369 def matrix_inner_product(cls
,X
,Y
):
1370 Xu
= cls
.real_unembed(X
)
1371 Yu
= cls
.real_unembed(Y
)
1372 tr
= (Xu
*Yu
).trace()
1375 # Works in QQ, AA, RDF, et cetera.
1377 except AttributeError:
1378 # A quaternion doesn't have a real() method, but does
1379 # have coefficient_tuple() method that returns the
1380 # coefficients of 1, i, j, and k -- in that order.
1381 return tr
.coefficient_tuple()[0]
1384 class RealMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1388 The identity function, for embedding real matrices into real
1394 def real_unembed(M
):
1396 The identity function, for unembedding real matrices from real
1402 class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra
,
1403 ConcreteEuclideanJordanAlgebra
):
1405 The rank-n simple EJA consisting of real symmetric n-by-n
1406 matrices, the usual symmetric Jordan product, and the trace inner
1407 product. It has dimension `(n^2 + n)/2` over the reals.
1411 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1415 sage: J = RealSymmetricEJA(2)
1416 sage: e0, e1, e2 = J.gens()
1424 In theory, our "field" can be any subfield of the reals::
1426 sage: RealSymmetricEJA(2, RDF)
1427 Euclidean Jordan algebra of dimension 3 over Real Double Field
1428 sage: RealSymmetricEJA(2, RR)
1429 Euclidean Jordan algebra of dimension 3 over Real Field with
1430 53 bits of precision
1434 The dimension of this algebra is `(n^2 + n) / 2`::
1436 sage: set_random_seed()
1437 sage: n_max = RealSymmetricEJA._max_random_instance_size()
1438 sage: n = ZZ.random_element(1, n_max)
1439 sage: J = RealSymmetricEJA(n)
1440 sage: J.dimension() == (n^2 + n)/2
1443 The Jordan multiplication is what we think it is::
1445 sage: set_random_seed()
1446 sage: J = RealSymmetricEJA.random_instance()
1447 sage: x,y = J.random_elements(2)
1448 sage: actual = (x*y).to_matrix()
1449 sage: X = x.to_matrix()
1450 sage: Y = y.to_matrix()
1451 sage: expected = (X*Y + Y*X)/2
1452 sage: actual == expected
1454 sage: J(expected) == x*y
1457 We can change the generator prefix::
1459 sage: RealSymmetricEJA(3, prefix='q').gens()
1460 (q0, q1, q2, q3, q4, q5)
1462 We can construct the (trivial) algebra of rank zero::
1464 sage: RealSymmetricEJA(0)
1465 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1469 def _denormalized_basis(cls
, n
, field
):
1471 Return a basis for the space of real symmetric n-by-n matrices.
1475 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1479 sage: set_random_seed()
1480 sage: n = ZZ.random_element(1,5)
1481 sage: B = RealSymmetricEJA._denormalized_basis(n,QQ)
1482 sage: all( M.is_symmetric() for M in B)
1486 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1490 for j
in range(i
+1):
1491 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
1495 Sij
= Eij
+ Eij
.transpose()
1501 def _max_random_instance_size():
1502 return 4 # Dimension 10
1505 def random_instance(cls
, field
=AA
, **kwargs
):
1507 Return a random instance of this type of algebra.
1509 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1510 return cls(n
, field
, **kwargs
)
1512 def __init__(self
, n
, field
=AA
, **kwargs
):
1513 basis
= self
._denormalized
_basis
(n
, field
)
1514 super(RealSymmetricEJA
, self
).__init
__(field
,
1518 self
.rank
.set_cache(n
)
1521 class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1525 Embed the n-by-n complex matrix ``M`` into the space of real
1526 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1527 bi` to the block matrix ``[[a,b],[-b,a]]``.
1531 sage: from mjo.eja.eja_algebra import \
1532 ....: ComplexMatrixEuclideanJordanAlgebra
1536 sage: F = QuadraticField(-1, 'I')
1537 sage: x1 = F(4 - 2*i)
1538 sage: x2 = F(1 + 2*i)
1541 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1542 sage: ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1551 Embedding is a homomorphism (isomorphism, in fact)::
1553 sage: set_random_seed()
1554 sage: n = ZZ.random_element(3)
1555 sage: F = QuadraticField(-1, 'I')
1556 sage: X = random_matrix(F, n)
1557 sage: Y = random_matrix(F, n)
1558 sage: Xe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X)
1559 sage: Ye = ComplexMatrixEuclideanJordanAlgebra.real_embed(Y)
1560 sage: XYe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1567 raise ValueError("the matrix 'M' must be square")
1569 # We don't need any adjoined elements...
1570 field
= M
.base_ring().base_ring()
1574 a
= z
.list()[0] # real part, I guess
1575 b
= z
.list()[1] # imag part, I guess
1576 blocks
.append(matrix(field
, 2, [[a
,b
],[-b
,a
]]))
1578 return matrix
.block(field
, n
, blocks
)
1582 def real_unembed(M
):
1584 The inverse of _embed_complex_matrix().
1588 sage: from mjo.eja.eja_algebra import \
1589 ....: ComplexMatrixEuclideanJordanAlgebra
1593 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1594 ....: [-2, 1, -4, 3],
1595 ....: [ 9, 10, 11, 12],
1596 ....: [-10, 9, -12, 11] ])
1597 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(A)
1599 [ 10*I + 9 12*I + 11]
1603 Unembedding is the inverse of embedding::
1605 sage: set_random_seed()
1606 sage: F = QuadraticField(-1, 'I')
1607 sage: M = random_matrix(F, 3)
1608 sage: Me = ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1609 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1615 raise ValueError("the matrix 'M' must be square")
1616 if not n
.mod(2).is_zero():
1617 raise ValueError("the matrix 'M' must be a complex embedding")
1619 # If "M" was normalized, its base ring might have roots
1620 # adjoined and they can stick around after unembedding.
1621 field
= M
.base_ring()
1622 R
= PolynomialRing(field
, 'z')
1625 # Sage doesn't know how to embed AA into QQbar, i.e. how
1626 # to adjoin sqrt(-1) to AA.
1629 F
= field
.extension(z
**2 + 1, 'I', embedding
=CLF(-1).sqrt())
1632 # Go top-left to bottom-right (reading order), converting every
1633 # 2-by-2 block we see to a single complex element.
1635 for k
in range(n
/2):
1636 for j
in range(n
/2):
1637 submat
= M
[2*k
:2*k
+2,2*j
:2*j
+2]
1638 if submat
[0,0] != submat
[1,1]:
1639 raise ValueError('bad on-diagonal submatrix')
1640 if submat
[0,1] != -submat
[1,0]:
1641 raise ValueError('bad off-diagonal submatrix')
1642 z
= submat
[0,0] + submat
[0,1]*i
1645 return matrix(F
, n
/2, elements
)
1649 def matrix_inner_product(cls
,X
,Y
):
1651 Compute a matrix inner product in this algebra directly from
1656 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1660 This gives the same answer as the slow, default method implemented
1661 in :class:`MatrixEuclideanJordanAlgebra`::
1663 sage: set_random_seed()
1664 sage: J = ComplexHermitianEJA.random_instance()
1665 sage: x,y = J.random_elements(2)
1666 sage: Xe = x.to_matrix()
1667 sage: Ye = y.to_matrix()
1668 sage: X = ComplexHermitianEJA.real_unembed(Xe)
1669 sage: Y = ComplexHermitianEJA.real_unembed(Ye)
1670 sage: expected = (X*Y).trace().real()
1671 sage: actual = ComplexHermitianEJA.matrix_inner_product(Xe,Ye)
1672 sage: actual == expected
1676 return RealMatrixEuclideanJordanAlgebra
.matrix_inner_product(X
,Y
)/2
1679 class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra
,
1680 ConcreteEuclideanJordanAlgebra
):
1682 The rank-n simple EJA consisting of complex Hermitian n-by-n
1683 matrices over the real numbers, the usual symmetric Jordan product,
1684 and the real-part-of-trace inner product. It has dimension `n^2` over
1689 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1693 In theory, our "field" can be any subfield of the reals::
1695 sage: ComplexHermitianEJA(2, RDF)
1696 Euclidean Jordan algebra of dimension 4 over Real Double Field
1697 sage: ComplexHermitianEJA(2, RR)
1698 Euclidean Jordan algebra of dimension 4 over Real Field with
1699 53 bits of precision
1703 The dimension of this algebra is `n^2`::
1705 sage: set_random_seed()
1706 sage: n_max = ComplexHermitianEJA._max_random_instance_size()
1707 sage: n = ZZ.random_element(1, n_max)
1708 sage: J = ComplexHermitianEJA(n)
1709 sage: J.dimension() == n^2
1712 The Jordan multiplication is what we think it is::
1714 sage: set_random_seed()
1715 sage: J = ComplexHermitianEJA.random_instance()
1716 sage: x,y = J.random_elements(2)
1717 sage: actual = (x*y).to_matrix()
1718 sage: X = x.to_matrix()
1719 sage: Y = y.to_matrix()
1720 sage: expected = (X*Y + Y*X)/2
1721 sage: actual == expected
1723 sage: J(expected) == x*y
1726 We can change the generator prefix::
1728 sage: ComplexHermitianEJA(2, prefix='z').gens()
1731 We can construct the (trivial) algebra of rank zero::
1733 sage: ComplexHermitianEJA(0)
1734 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1739 def _denormalized_basis(cls
, n
, field
):
1741 Returns a basis for the space of complex Hermitian n-by-n matrices.
1743 Why do we embed these? Basically, because all of numerical linear
1744 algebra assumes that you're working with vectors consisting of `n`
1745 entries from a field and scalars from the same field. There's no way
1746 to tell SageMath that (for example) the vectors contain complex
1747 numbers, while the scalar field is real.
1751 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1755 sage: set_random_seed()
1756 sage: n = ZZ.random_element(1,5)
1757 sage: field = QuadraticField(2, 'sqrt2')
1758 sage: B = ComplexHermitianEJA._denormalized_basis(n, field)
1759 sage: all( M.is_symmetric() for M in B)
1763 R
= PolynomialRing(field
, 'z')
1765 F
= field
.extension(z
**2 + 1, 'I')
1768 # This is like the symmetric case, but we need to be careful:
1770 # * We want conjugate-symmetry, not just symmetry.
1771 # * The diagonal will (as a result) be real.
1775 for j
in range(i
+1):
1776 Eij
= matrix(F
, n
, lambda k
,l
: k
==i
and l
==j
)
1778 Sij
= cls
.real_embed(Eij
)
1781 # The second one has a minus because it's conjugated.
1782 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
1784 Sij_imag
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
1787 # Since we embedded these, we can drop back to the "field" that we
1788 # started with instead of the complex extension "F".
1789 return ( s
.change_ring(field
) for s
in S
)
1792 def __init__(self
, n
, field
=AA
, **kwargs
):
1793 basis
= self
._denormalized
_basis
(n
,field
)
1794 super(ComplexHermitianEJA
,self
).__init
__(field
,
1798 self
.rank
.set_cache(n
)
1801 def _max_random_instance_size():
1802 return 3 # Dimension 9
1805 def random_instance(cls
, field
=AA
, **kwargs
):
1807 Return a random instance of this type of algebra.
1809 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1810 return cls(n
, field
, **kwargs
)
1812 class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1816 Embed the n-by-n quaternion matrix ``M`` into the space of real
1817 matrices of size 4n-by-4n by first sending each quaternion entry `z
1818 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
1819 c+di],[-c + di, a-bi]]`, and then embedding those into a real
1824 sage: from mjo.eja.eja_algebra import \
1825 ....: QuaternionMatrixEuclideanJordanAlgebra
1829 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1830 sage: i,j,k = Q.gens()
1831 sage: x = 1 + 2*i + 3*j + 4*k
1832 sage: M = matrix(Q, 1, [[x]])
1833 sage: QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1839 Embedding is a homomorphism (isomorphism, in fact)::
1841 sage: set_random_seed()
1842 sage: n = ZZ.random_element(2)
1843 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1844 sage: X = random_matrix(Q, n)
1845 sage: Y = random_matrix(Q, n)
1846 sage: Xe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X)
1847 sage: Ye = QuaternionMatrixEuclideanJordanAlgebra.real_embed(Y)
1848 sage: XYe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1853 quaternions
= M
.base_ring()
1856 raise ValueError("the matrix 'M' must be square")
1858 F
= QuadraticField(-1, 'I')
1863 t
= z
.coefficient_tuple()
1868 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
1869 [-c
+ d
*i
, a
- b
*i
]])
1870 realM
= ComplexMatrixEuclideanJordanAlgebra
.real_embed(cplxM
)
1871 blocks
.append(realM
)
1873 # We should have real entries by now, so use the realest field
1874 # we've got for the return value.
1875 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
1880 def real_unembed(M
):
1882 The inverse of _embed_quaternion_matrix().
1886 sage: from mjo.eja.eja_algebra import \
1887 ....: QuaternionMatrixEuclideanJordanAlgebra
1891 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
1892 ....: [-2, 1, -4, 3],
1893 ....: [-3, 4, 1, -2],
1894 ....: [-4, -3, 2, 1]])
1895 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(M)
1896 [1 + 2*i + 3*j + 4*k]
1900 Unembedding is the inverse of embedding::
1902 sage: set_random_seed()
1903 sage: Q = QuaternionAlgebra(QQ, -1, -1)
1904 sage: M = random_matrix(Q, 3)
1905 sage: Me = QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1906 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1912 raise ValueError("the matrix 'M' must be square")
1913 if not n
.mod(4).is_zero():
1914 raise ValueError("the matrix 'M' must be a quaternion embedding")
1916 # Use the base ring of the matrix to ensure that its entries can be
1917 # multiplied by elements of the quaternion algebra.
1918 field
= M
.base_ring()
1919 Q
= QuaternionAlgebra(field
,-1,-1)
1922 # Go top-left to bottom-right (reading order), converting every
1923 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
1926 for l
in range(n
/4):
1927 for m
in range(n
/4):
1928 submat
= ComplexMatrixEuclideanJordanAlgebra
.real_unembed(
1929 M
[4*l
:4*l
+4,4*m
:4*m
+4] )
1930 if submat
[0,0] != submat
[1,1].conjugate():
1931 raise ValueError('bad on-diagonal submatrix')
1932 if submat
[0,1] != -submat
[1,0].conjugate():
1933 raise ValueError('bad off-diagonal submatrix')
1934 z
= submat
[0,0].real()
1935 z
+= submat
[0,0].imag()*i
1936 z
+= submat
[0,1].real()*j
1937 z
+= submat
[0,1].imag()*k
1940 return matrix(Q
, n
/4, elements
)
1944 def matrix_inner_product(cls
,X
,Y
):
1946 Compute a matrix inner product in this algebra directly from
1951 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1955 This gives the same answer as the slow, default method implemented
1956 in :class:`MatrixEuclideanJordanAlgebra`::
1958 sage: set_random_seed()
1959 sage: J = QuaternionHermitianEJA.random_instance()
1960 sage: x,y = J.random_elements(2)
1961 sage: Xe = x.to_matrix()
1962 sage: Ye = y.to_matrix()
1963 sage: X = QuaternionHermitianEJA.real_unembed(Xe)
1964 sage: Y = QuaternionHermitianEJA.real_unembed(Ye)
1965 sage: expected = (X*Y).trace().coefficient_tuple()[0]
1966 sage: actual = QuaternionHermitianEJA.matrix_inner_product(Xe,Ye)
1967 sage: actual == expected
1971 return RealMatrixEuclideanJordanAlgebra
.matrix_inner_product(X
,Y
)/4
1974 class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra
,
1975 ConcreteEuclideanJordanAlgebra
):
1977 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
1978 matrices, the usual symmetric Jordan product, and the
1979 real-part-of-trace inner product. It has dimension `2n^2 - n` over
1984 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1988 In theory, our "field" can be any subfield of the reals::
1990 sage: QuaternionHermitianEJA(2, RDF)
1991 Euclidean Jordan algebra of dimension 6 over Real Double Field
1992 sage: QuaternionHermitianEJA(2, RR)
1993 Euclidean Jordan algebra of dimension 6 over Real Field with
1994 53 bits of precision
1998 The dimension of this algebra is `2*n^2 - n`::
2000 sage: set_random_seed()
2001 sage: n_max = QuaternionHermitianEJA._max_random_instance_size()
2002 sage: n = ZZ.random_element(1, n_max)
2003 sage: J = QuaternionHermitianEJA(n)
2004 sage: J.dimension() == 2*(n^2) - n
2007 The Jordan multiplication is what we think it is::
2009 sage: set_random_seed()
2010 sage: J = QuaternionHermitianEJA.random_instance()
2011 sage: x,y = J.random_elements(2)
2012 sage: actual = (x*y).to_matrix()
2013 sage: X = x.to_matrix()
2014 sage: Y = y.to_matrix()
2015 sage: expected = (X*Y + Y*X)/2
2016 sage: actual == expected
2018 sage: J(expected) == x*y
2021 We can change the generator prefix::
2023 sage: QuaternionHermitianEJA(2, prefix='a').gens()
2024 (a0, a1, a2, a3, a4, a5)
2026 We can construct the (trivial) algebra of rank zero::
2028 sage: QuaternionHermitianEJA(0)
2029 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2033 def _denormalized_basis(cls
, n
, field
):
2035 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
2037 Why do we embed these? Basically, because all of numerical
2038 linear algebra assumes that you're working with vectors consisting
2039 of `n` entries from a field and scalars from the same field. There's
2040 no way to tell SageMath that (for example) the vectors contain
2041 complex numbers, while the scalar field is real.
2045 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2049 sage: set_random_seed()
2050 sage: n = ZZ.random_element(1,5)
2051 sage: B = QuaternionHermitianEJA._denormalized_basis(n,QQ)
2052 sage: all( M.is_symmetric() for M in B )
2056 Q
= QuaternionAlgebra(QQ
,-1,-1)
2059 # This is like the symmetric case, but we need to be careful:
2061 # * We want conjugate-symmetry, not just symmetry.
2062 # * The diagonal will (as a result) be real.
2066 for j
in range(i
+1):
2067 Eij
= matrix(Q
, n
, lambda k
,l
: k
==i
and l
==j
)
2069 Sij
= cls
.real_embed(Eij
)
2072 # The second, third, and fourth ones have a minus
2073 # because they're conjugated.
2074 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
2076 Sij_I
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
2078 Sij_J
= cls
.real_embed(J
*Eij
- J
*Eij
.transpose())
2080 Sij_K
= cls
.real_embed(K
*Eij
- K
*Eij
.transpose())
2083 # Since we embedded these, we can drop back to the "field" that we
2084 # started with instead of the quaternion algebra "Q".
2085 return ( s
.change_ring(field
) for s
in S
)
2088 def __init__(self
, n
, field
=AA
, **kwargs
):
2089 basis
= self
._denormalized
_basis
(n
,field
)
2090 super(QuaternionHermitianEJA
,self
).__init
__(field
,
2094 self
.rank
.set_cache(n
)
2097 def _max_random_instance_size():
2099 The maximum rank of a random QuaternionHermitianEJA.
2101 return 2 # Dimension 6
2104 def random_instance(cls
, field
=AA
, **kwargs
):
2106 Return a random instance of this type of algebra.
2108 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2109 return cls(n
, field
, **kwargs
)
2112 class HadamardEJA(RationalBasisEuclideanJordanAlgebra
,
2113 ConcreteEuclideanJordanAlgebra
):
2115 Return the Euclidean Jordan Algebra corresponding to the set
2116 `R^n` under the Hadamard product.
2118 Note: this is nothing more than the Cartesian product of ``n``
2119 copies of the spin algebra. Once Cartesian product algebras
2120 are implemented, this can go.
2124 sage: from mjo.eja.eja_algebra import HadamardEJA
2128 This multiplication table can be verified by hand::
2130 sage: J = HadamardEJA(3)
2131 sage: e0,e1,e2 = J.gens()
2147 We can change the generator prefix::
2149 sage: HadamardEJA(3, prefix='r').gens()
2153 def __init__(self
, n
, field
=AA
, **kwargs
):
2154 V
= VectorSpace(field
, n
)
2155 mult_table
= [ [ V
.gen(i
)*(i
== j
) for j
in range(n
) ]
2158 # Inner products are real numbers and not algebra
2159 # elements, so once we turn the algebra element
2160 # into a vector in inner_product(), we never go
2161 # back. As a result -- contrary to what we do with
2162 # self._multiplication_table -- we store the inner
2163 # product table as a plain old matrix and not as
2164 # an algebra operator.
2165 ip_table
= matrix
.identity(field
,n
)
2166 self
._inner
_product
_matrix
= ip_table
2168 super(HadamardEJA
, self
).__init
__(field
,
2172 self
.rank
.set_cache(n
)
2175 self
.one
.set_cache( self
.zero() )
2177 self
.one
.set_cache( sum(self
.gens()) )
2180 def _max_random_instance_size():
2182 The maximum dimension of a random HadamardEJA.
2187 def random_instance(cls
, field
=AA
, **kwargs
):
2189 Return a random instance of this type of algebra.
2191 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2192 return cls(n
, field
, **kwargs
)
2195 class BilinearFormEJA(RationalBasisEuclideanJordanAlgebra
,
2196 ConcreteEuclideanJordanAlgebra
):
2198 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2199 with the half-trace inner product and jordan product ``x*y =
2200 (<Bx,y>,y_bar>, x0*y_bar + y0*x_bar)`` where `B = 1 \times B22` is
2201 a symmetric positive-definite "bilinear form" matrix. Its
2202 dimension is the size of `B`, and it has rank two in dimensions
2203 larger than two. It reduces to the ``JordanSpinEJA`` when `B` is
2204 the identity matrix of order ``n``.
2206 We insist that the one-by-one upper-left identity block of `B` be
2207 passed in as well so that we can be passed a matrix of size zero
2208 to construct a trivial algebra.
2212 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
2213 ....: JordanSpinEJA)
2217 When no bilinear form is specified, the identity matrix is used,
2218 and the resulting algebra is the Jordan spin algebra::
2220 sage: B = matrix.identity(AA,3)
2221 sage: J0 = BilinearFormEJA(B)
2222 sage: J1 = JordanSpinEJA(3)
2223 sage: J0.multiplication_table() == J0.multiplication_table()
2226 An error is raised if the matrix `B` does not correspond to a
2227 positive-definite bilinear form::
2229 sage: B = matrix.random(QQ,2,3)
2230 sage: J = BilinearFormEJA(B)
2231 Traceback (most recent call last):
2233 ValueError: bilinear form is not positive-definite
2234 sage: B = matrix.zero(QQ,3)
2235 sage: J = BilinearFormEJA(B)
2236 Traceback (most recent call last):
2238 ValueError: bilinear form is not positive-definite
2242 We can create a zero-dimensional algebra::
2244 sage: B = matrix.identity(AA,0)
2245 sage: J = BilinearFormEJA(B)
2249 We can check the multiplication condition given in the Jordan, von
2250 Neumann, and Wigner paper (and also discussed on my "On the
2251 symmetry..." paper). Note that this relies heavily on the standard
2252 choice of basis, as does anything utilizing the bilinear form matrix::
2254 sage: set_random_seed()
2255 sage: n = ZZ.random_element(5)
2256 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2257 sage: B11 = matrix.identity(QQ,1)
2258 sage: B22 = M.transpose()*M
2259 sage: B = block_matrix(2,2,[ [B11,0 ],
2261 sage: J = BilinearFormEJA(B)
2262 sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
2263 sage: V = J.vector_space()
2264 sage: sis = [ J.from_vector(V([0] + (M.inverse()*ei).list()))
2265 ....: for ei in eis ]
2266 sage: actual = [ sis[i]*sis[j]
2267 ....: for i in range(n-1)
2268 ....: for j in range(n-1) ]
2269 sage: expected = [ J.one() if i == j else J.zero()
2270 ....: for i in range(n-1)
2271 ....: for j in range(n-1) ]
2272 sage: actual == expected
2275 def __init__(self
, B
, field
=AA
, **kwargs
):
2278 if not B
.is_positive_definite():
2279 raise ValueError("bilinear form is not positive-definite")
2281 V
= VectorSpace(field
, n
)
2282 mult_table
= [[V
.zero() for j
in range(n
)] for i
in range(n
)]
2291 z0
= (B
*x
).inner_product(y
)
2292 zbar
= y0
*xbar
+ x0
*ybar
2293 z
= V([z0
] + zbar
.list())
2294 mult_table
[i
][j
] = z
2296 # Inner products are real numbers and not algebra
2297 # elements, so once we turn the algebra element
2298 # into a vector in inner_product(), we never go
2299 # back. As a result -- contrary to what we do with
2300 # self._multiplication_table -- we store the inner
2301 # product table as a plain old matrix and not as
2302 # an algebra operator.
2304 self
._inner
_product
_matrix
= ip_table
2306 super(BilinearFormEJA
, self
).__init
__(field
,
2311 # The rank of this algebra is two, unless we're in a
2312 # one-dimensional ambient space (because the rank is bounded
2313 # by the ambient dimension).
2314 self
.rank
.set_cache(min(n
,2))
2317 self
.one
.set_cache( self
.zero() )
2319 self
.one
.set_cache( self
.monomial(0) )
2322 def _max_random_instance_size():
2324 The maximum dimension of a random BilinearFormEJA.
2329 def random_instance(cls
, field
=AA
, **kwargs
):
2331 Return a random instance of this algebra.
2333 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2335 B
= matrix
.identity(field
, n
)
2336 return cls(B
, field
, **kwargs
)
2338 B11
= matrix
.identity(field
,1)
2339 M
= matrix
.random(field
, n
-1)
2340 I
= matrix
.identity(field
, n
-1)
2341 alpha
= field
.zero()
2342 while alpha
.is_zero():
2343 alpha
= field
.random_element().abs()
2344 B22
= M
.transpose()*M
+ alpha
*I
2346 from sage
.matrix
.special
import block_matrix
2347 B
= block_matrix(2,2, [ [B11
, ZZ(0) ],
2350 return cls(B
, field
, **kwargs
)
2353 class JordanSpinEJA(BilinearFormEJA
):
2355 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2356 with the usual inner product and jordan product ``x*y =
2357 (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2362 sage: from mjo.eja.eja_algebra import JordanSpinEJA
2366 This multiplication table can be verified by hand::
2368 sage: J = JordanSpinEJA(4)
2369 sage: e0,e1,e2,e3 = J.gens()
2385 We can change the generator prefix::
2387 sage: JordanSpinEJA(2, prefix='B').gens()
2392 Ensure that we have the usual inner product on `R^n`::
2394 sage: set_random_seed()
2395 sage: J = JordanSpinEJA.random_instance()
2396 sage: x,y = J.random_elements(2)
2397 sage: actual = x.inner_product(y)
2398 sage: expected = x.to_vector().inner_product(y.to_vector())
2399 sage: actual == expected
2403 def __init__(self
, n
, field
=AA
, **kwargs
):
2404 # This is a special case of the BilinearFormEJA with the identity
2405 # matrix as its bilinear form.
2406 B
= matrix
.identity(field
, n
)
2407 super(JordanSpinEJA
, self
).__init
__(B
, field
, **kwargs
)
2410 def _max_random_instance_size():
2412 The maximum dimension of a random JordanSpinEJA.
2417 def random_instance(cls
, field
=AA
, **kwargs
):
2419 Return a random instance of this type of algebra.
2421 Needed here to override the implementation for ``BilinearFormEJA``.
2423 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2424 return cls(n
, field
, **kwargs
)
2427 class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra
,
2428 ConcreteEuclideanJordanAlgebra
):
2430 The trivial Euclidean Jordan algebra consisting of only a zero element.
2434 sage: from mjo.eja.eja_algebra import TrivialEJA
2438 sage: J = TrivialEJA()
2445 sage: 7*J.one()*12*J.one()
2447 sage: J.one().inner_product(J.one())
2449 sage: J.one().norm()
2451 sage: J.one().subalgebra_generated_by()
2452 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2457 def __init__(self
, field
=AA
, **kwargs
):
2459 self
._inner
_product
_matrix
= matrix(field
,0)
2460 super(TrivialEJA
, self
).__init
__(field
,
2464 # The rank is zero using my definition, namely the dimension of the
2465 # largest subalgebra generated by any element.
2466 self
.rank
.set_cache(0)
2467 self
.one
.set_cache( self
.zero() )
2470 def random_instance(cls
, field
=AA
, **kwargs
):
2471 # We don't take a "size" argument so the superclass method is
2472 # inappropriate for us.
2473 return cls(field
, **kwargs
)
2475 class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra
):
2477 The external (orthogonal) direct sum of two other Euclidean Jordan
2478 algebras. Essentially the Cartesian product of its two factors.
2479 Every Euclidean Jordan algebra decomposes into an orthogonal
2480 direct sum of simple Euclidean Jordan algebras, so no generality
2481 is lost by providing only this construction.
2485 sage: from mjo.eja.eja_algebra import (random_eja,
2487 ....: RealSymmetricEJA,
2492 sage: J1 = HadamardEJA(2)
2493 sage: J2 = RealSymmetricEJA(3)
2494 sage: J = DirectSumEJA(J1,J2)
2502 The external direct sum construction is only valid when the two factors
2503 have the same base ring; an error is raised otherwise::
2505 sage: set_random_seed()
2506 sage: J1 = random_eja(AA)
2507 sage: J2 = random_eja(QQ)
2508 sage: J = DirectSumEJA(J1,J2)
2509 Traceback (most recent call last):
2511 ValueError: algebras must share the same base field
2514 def __init__(self
, J1
, J2
, **kwargs
):
2515 if J1
.base_ring() != J2
.base_ring():
2516 raise ValueError("algebras must share the same base field")
2517 field
= J1
.base_ring()
2519 self
._factors
= (J1
, J2
)
2523 V
= VectorSpace(field
, n
)
2524 mult_table
= [ [ V
.zero() for j
in range(n
) ]
2528 p
= (J1
.monomial(i
)*J1
.monomial(j
)).to_vector()
2529 mult_table
[i
][j
] = V(p
.list() + [field
.zero()]*n2
)
2533 p
= (J2
.monomial(i
)*J2
.monomial(j
)).to_vector()
2534 mult_table
[n1
+i
][n1
+j
] = V([field
.zero()]*n1
+ p
.list())
2536 super(DirectSumEJA
, self
).__init
__(field
,
2540 self
.rank
.set_cache(J1
.rank() + J2
.rank())
2545 Return the pair of this algebra's factors.
2549 sage: from mjo.eja.eja_algebra import (HadamardEJA,
2550 ....: JordanSpinEJA,
2555 sage: J1 = HadamardEJA(2,QQ)
2556 sage: J2 = JordanSpinEJA(3,QQ)
2557 sage: J = DirectSumEJA(J1,J2)
2559 (Euclidean Jordan algebra of dimension 2 over Rational Field,
2560 Euclidean Jordan algebra of dimension 3 over Rational Field)
2563 return self
._factors
2565 def projections(self
):
2567 Return a pair of projections onto this algebra's factors.
2571 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
2572 ....: ComplexHermitianEJA,
2577 sage: J1 = JordanSpinEJA(2)
2578 sage: J2 = ComplexHermitianEJA(2)
2579 sage: J = DirectSumEJA(J1,J2)
2580 sage: (pi_left, pi_right) = J.projections()
2581 sage: J.one().to_vector()
2583 sage: pi_left(J.one()).to_vector()
2585 sage: pi_right(J.one()).to_vector()
2589 (J1
,J2
) = self
.factors()
2592 V_basis
= self
.vector_space().basis()
2593 # Need to specify the dimensions explicitly so that we don't
2594 # wind up with a zero-by-zero matrix when we want e.g. a
2595 # zero-by-two matrix (important for composing things).
2596 P1
= matrix(self
.base_ring(), m
, m
+n
, V_basis
[:m
])
2597 P2
= matrix(self
.base_ring(), n
, m
+n
, V_basis
[m
:])
2598 pi_left
= FiniteDimensionalEuclideanJordanAlgebraOperator(self
,J1
,P1
)
2599 pi_right
= FiniteDimensionalEuclideanJordanAlgebraOperator(self
,J2
,P2
)
2600 return (pi_left
, pi_right
)
2602 def inclusions(self
):
2604 Return the pair of inclusion maps from our factors into us.
2608 sage: from mjo.eja.eja_algebra import (random_eja,
2609 ....: JordanSpinEJA,
2610 ....: RealSymmetricEJA,
2615 sage: J1 = JordanSpinEJA(3)
2616 sage: J2 = RealSymmetricEJA(2)
2617 sage: J = DirectSumEJA(J1,J2)
2618 sage: (iota_left, iota_right) = J.inclusions()
2619 sage: iota_left(J1.zero()) == J.zero()
2621 sage: iota_right(J2.zero()) == J.zero()
2623 sage: J1.one().to_vector()
2625 sage: iota_left(J1.one()).to_vector()
2627 sage: J2.one().to_vector()
2629 sage: iota_right(J2.one()).to_vector()
2631 sage: J.one().to_vector()
2636 Composing a projection with the corresponding inclusion should
2637 produce the identity map, and mismatching them should produce
2640 sage: set_random_seed()
2641 sage: J1 = random_eja()
2642 sage: J2 = random_eja()
2643 sage: J = DirectSumEJA(J1,J2)
2644 sage: (iota_left, iota_right) = J.inclusions()
2645 sage: (pi_left, pi_right) = J.projections()
2646 sage: pi_left*iota_left == J1.one().operator()
2648 sage: pi_right*iota_right == J2.one().operator()
2650 sage: (pi_left*iota_right).is_zero()
2652 sage: (pi_right*iota_left).is_zero()
2656 (J1
,J2
) = self
.factors()
2659 V_basis
= self
.vector_space().basis()
2660 # Need to specify the dimensions explicitly so that we don't
2661 # wind up with a zero-by-zero matrix when we want e.g. a
2662 # two-by-zero matrix (important for composing things).
2663 I1
= matrix
.column(self
.base_ring(), m
, m
+n
, V_basis
[:m
])
2664 I2
= matrix
.column(self
.base_ring(), n
, m
+n
, V_basis
[m
:])
2665 iota_left
= FiniteDimensionalEuclideanJordanAlgebraOperator(J1
,self
,I1
)
2666 iota_right
= FiniteDimensionalEuclideanJordanAlgebraOperator(J2
,self
,I2
)
2667 return (iota_left
, iota_right
)
2669 def inner_product(self
, x
, y
):
2671 The standard Cartesian inner-product.
2673 We project ``x`` and ``y`` onto our factors, and add up the
2674 inner-products from the subalgebras.
2679 sage: from mjo.eja.eja_algebra import (HadamardEJA,
2680 ....: QuaternionHermitianEJA,
2685 sage: J1 = HadamardEJA(3,QQ)
2686 sage: J2 = QuaternionHermitianEJA(2,QQ,normalize_basis=False)
2687 sage: J = DirectSumEJA(J1,J2)
2692 sage: x1.inner_product(x2)
2694 sage: y1.inner_product(y2)
2696 sage: J.one().inner_product(J.one())
2700 (pi_left
, pi_right
) = self
.projections()
2706 return (x1
.inner_product(y1
) + x2
.inner_product(y2
))
2710 random_eja
= ConcreteEuclideanJordanAlgebra
.random_instance