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 A
= matrix(field
, vector_basis
)
1066 vector_basis
= gram_schmidt(vector_basis
, inner_product
)
1067 W
= V
.span_of_basis( vector_basis
)
1068 Q
= matrix(field
, vector_basis
)
1069 # A = QR <==> A.T == R.T*Q.T
1070 # So, Q.solve_right() is equivalent to the Q.T.solve_left()
1072 self
._deorthonormalization
_matrix
= Q
.solve_right(A
)
1074 if basis_is_matrices
:
1075 from mjo
.eja
.eja_utils
import _vec2mat
1076 basis
= tuple( map(_vec2mat
,vector_basis
) )
1078 W
= V
.span_of_basis( vector_basis
)
1080 mult_table
= [ [0 for i
in range(n
)] for j
in range(n
) ]
1081 ip_table
= [ [0 for i
in range(n
)] for j
in range(n
) ]
1083 # Note: the Jordan and inner- products are defined in terms
1084 # of the ambient basis. It's important that their arguments
1085 # are in ambient coordinates as well.
1087 for j
in range(i
+1):
1088 # ortho basis w.r.t. ambient coords
1089 q_i
= vector_basis
[i
]
1090 q_j
= vector_basis
[j
]
1092 if basis_is_matrices
:
1096 elt
= jordan_product(q_i
, q_j
)
1097 ip
= inner_product(q_i
, q_j
)
1099 if basis_is_matrices
:
1100 # do another mat2vec because the multiplication
1101 # table is in terms of vectors
1104 elt
= W
.coordinate_vector(elt
)
1105 mult_table
[i
][j
] = elt
1106 mult_table
[j
][i
] = elt
1110 self
._inner
_product
_matrix
= matrix(field
,ip_table
)
1112 if basis_is_matrices
:
1116 basis
= tuple( x
.column() for x
in basis
)
1118 super().__init
__(field
,
1122 basis
, # matrix basis
1126 class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra
):
1128 Algebras whose basis consists of vectors with rational
1129 entries. Equivalently, algebras whose multiplication tables
1130 contain only rational coefficients.
1132 When an EJA has a basis that can be made rational, we can speed up
1133 the computation of its characteristic polynomial by doing it over
1134 ``QQ``. All of the named EJA constructors that we provide fall
1138 def _charpoly_coefficients(self
):
1140 Override the parent method with something that tries to compute
1141 over a faster (non-extension) field.
1145 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1149 The base ring of the resulting polynomial coefficients is what
1150 it should be, and not the rationals (unless the algebra was
1151 already over the rationals)::
1153 sage: J = JordanSpinEJA(3)
1154 sage: J._charpoly_coefficients()
1155 (X1^2 - X2^2 - X3^2, -2*X1)
1156 sage: a0 = J._charpoly_coefficients()[0]
1158 Algebraic Real Field
1159 sage: a0.base_ring()
1160 Algebraic Real Field
1163 if self
.base_ring() is QQ
:
1164 # There's no need to construct *another* algebra over the
1165 # rationals if this one is already over the rationals.
1166 superclass
= super(RationalBasisEuclideanJordanAlgebra
, self
)
1167 return superclass
._charpoly
_coefficients
()
1170 tuple(map(lambda x
: x
.to_vector(), ls
))
1171 for ls
in self
._multiplication
_table
1174 # Do the computation over the rationals. The answer will be
1175 # the same, because our basis coordinates are (essentially)
1177 J
= FiniteDimensionalEuclideanJordanAlgebra(QQ
,
1181 a
= J
._charpoly
_coefficients
()
1182 return tuple(map(lambda x
: x
.change_ring(self
.base_ring()), a
))
1185 class ConcreteEuclideanJordanAlgebra
:
1187 A class for the Euclidean Jordan algebras that we know by name.
1189 These are the Jordan algebras whose basis, multiplication table,
1190 rank, and so on are known a priori. More to the point, they are
1191 the Euclidean Jordan algebras for which we are able to conjure up
1192 a "random instance."
1196 sage: from mjo.eja.eja_algebra import ConcreteEuclideanJordanAlgebra
1200 Our basis is normalized with respect to the algebra's inner
1201 product, unless we specify otherwise::
1203 sage: set_random_seed()
1204 sage: J = ConcreteEuclideanJordanAlgebra.random_instance()
1205 sage: all( b.norm() == 1 for b in J.gens() )
1208 Since our basis is orthonormal with respect to the algebra's inner
1209 product, and since we know that this algebra is an EJA, any
1210 left-multiplication operator's matrix will be symmetric because
1211 natural->EJA basis representation is an isometry and within the
1212 EJA the operator is self-adjoint by the Jordan axiom::
1214 sage: set_random_seed()
1215 sage: J = ConcreteEuclideanJordanAlgebra.random_instance()
1216 sage: x = J.random_element()
1217 sage: x.operator().is_self_adjoint()
1222 def _max_random_instance_size():
1224 Return an integer "size" that is an upper bound on the size of
1225 this algebra when it is used in a random test
1226 case. Unfortunately, the term "size" is ambiguous -- when
1227 dealing with `R^n` under either the Hadamard or Jordan spin
1228 product, the "size" refers to the dimension `n`. When dealing
1229 with a matrix algebra (real symmetric or complex/quaternion
1230 Hermitian), it refers to the size of the matrix, which is far
1231 less than the dimension of the underlying vector space.
1233 This method must be implemented in each subclass.
1235 raise NotImplementedError
1238 def random_instance(cls
, field
=AA
, **kwargs
):
1240 Return a random instance of this type of algebra.
1242 This method should be implemented in each subclass.
1244 from sage
.misc
.prandom
import choice
1245 eja_class
= choice(cls
.__subclasses
__())
1246 return eja_class
.random_instance(field
)
1249 class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra
):
1251 def __init__(self
, field
, basis
, normalize_basis
=True, **kwargs
):
1253 Compared to the superclass constructor, we take a basis instead of
1254 a multiplication table because the latter can be computed in terms
1255 of the former when the product is known (like it is here).
1257 # Used in this class's fast _charpoly_coefficients() override.
1258 self
._basis
_normalizers
= None
1260 # We're going to loop through this a few times, so now's a good
1261 # time to ensure that it isn't a generator expression.
1262 basis
= tuple(basis
)
1264 algebra_dim
= len(basis
)
1265 degree
= 0 # size of the matrices
1267 degree
= basis
[0].nrows()
1269 if algebra_dim
> 1 and normalize_basis
:
1270 # We'll need sqrt(2) to normalize the basis, and this
1271 # winds up in the multiplication table, so the whole
1272 # algebra needs to be over the field extension.
1273 R
= PolynomialRing(field
, 'z')
1276 if p
.is_irreducible():
1277 field
= field
.extension(p
, 'sqrt2', embedding
=RLF(2).sqrt())
1278 basis
= tuple( s
.change_ring(field
) for s
in basis
)
1279 self
._basis
_normalizers
= tuple(
1280 ~
(self
.matrix_inner_product(s
,s
).sqrt()) for s
in basis
)
1281 basis
= tuple(s
*c
for (s
,c
) in zip(basis
,self
._basis
_normalizers
))
1283 # Now compute the multiplication and inner product tables.
1284 # We have to do this *after* normalizing the basis, because
1285 # scaling affects the answers.
1286 V
= VectorSpace(field
, degree
**2)
1287 W
= V
.span_of_basis( _mat2vec(s
) for s
in basis
)
1288 mult_table
= [[W
.zero() for j
in range(algebra_dim
)]
1289 for i
in range(algebra_dim
)]
1290 ip_table
= [[W
.zero() for j
in range(algebra_dim
)]
1291 for i
in range(algebra_dim
)]
1292 for i
in range(algebra_dim
):
1293 for j
in range(algebra_dim
):
1294 mat_entry
= (basis
[i
]*basis
[j
] + basis
[j
]*basis
[i
])/2
1295 mult_table
[i
][j
] = W
.coordinate_vector(_mat2vec(mat_entry
))
1298 # HACK: ignore the error here if we don't need the
1299 # inner product (as is the case when we construct
1300 # a dummy QQ-algebra for fast charpoly coefficients.
1301 ip_table
[i
][j
] = self
.matrix_inner_product(basis
[i
],
1308 self
._inner
_product
_matrix
= matrix(field
,ip_table
)
1312 super(MatrixEuclideanJordanAlgebra
, self
).__init
__(field
,
1317 if algebra_dim
== 0:
1318 self
.one
.set_cache(self
.zero())
1320 n
= basis
[0].nrows()
1321 # The identity wrt (A,B) -> (AB + BA)/2 is independent of the
1322 # details of this algebra.
1323 self
.one
.set_cache(self(matrix
.identity(field
,n
)))
1327 def _charpoly_coefficients(self
):
1329 Override the parent method with something that tries to compute
1330 over a faster (non-extension) field.
1332 if self
._basis
_normalizers
is None or self
.base_ring() is QQ
:
1333 # We didn't normalize, or the basis we started with had
1334 # entries in a nice field already. Just compute the thing.
1335 return super(MatrixEuclideanJordanAlgebra
, self
)._charpoly
_coefficients
()
1337 basis
= ( (b
/n
) for (b
,n
) in zip(self
.matrix_basis(),
1338 self
._basis
_normalizers
) )
1340 # Do this over the rationals and convert back at the end.
1341 # Only works because we know the entries of the basis are
1342 # integers. The argument ``check_axioms=False`` is required
1343 # because the trace inner-product method for this
1344 # class is a stub and can't actually be checked.
1345 J
= MatrixEuclideanJordanAlgebra(QQ
,
1347 normalize_basis
=False,
1350 a
= J
._charpoly
_coefficients
()
1352 # Unfortunately, changing the basis does change the
1353 # coefficients of the characteristic polynomial, but since
1354 # these are really the coefficients of the "characteristic
1355 # polynomial of" function, everything is still nice and
1356 # unevaluated. It's therefore "obvious" how scaling the
1357 # basis affects the coordinate variables X1, X2, et
1358 # cetera. Scaling the first basis vector up by "n" adds a
1359 # factor of 1/n into every "X1" term, for example. So here
1360 # we simply undo the basis_normalizer scaling that we
1361 # performed earlier.
1363 # The a[0] access here is safe because trivial algebras
1364 # won't have any basis normalizers and therefore won't
1365 # make it to this "else" branch.
1366 XS
= a
[0].parent().gens()
1367 subs_dict
= { XS
[i
]: self
._basis
_normalizers
[i
]*XS
[i
]
1368 for i
in range(len(XS
)) }
1369 return tuple( a_i
.subs(subs_dict
) for a_i
in a
)
1375 Embed the matrix ``M`` into a space of real matrices.
1377 The matrix ``M`` can have entries in any field at the moment:
1378 the real numbers, complex numbers, or quaternions. And although
1379 they are not a field, we can probably support octonions at some
1380 point, too. This function returns a real matrix that "acts like"
1381 the original with respect to matrix multiplication; i.e.
1383 real_embed(M*N) = real_embed(M)*real_embed(N)
1386 raise NotImplementedError
1390 def real_unembed(M
):
1392 The inverse of :meth:`real_embed`.
1394 raise NotImplementedError
1397 def matrix_inner_product(cls
,X
,Y
):
1398 Xu
= cls
.real_unembed(X
)
1399 Yu
= cls
.real_unembed(Y
)
1400 tr
= (Xu
*Yu
).trace()
1403 # Works in QQ, AA, RDF, et cetera.
1405 except AttributeError:
1406 # A quaternion doesn't have a real() method, but does
1407 # have coefficient_tuple() method that returns the
1408 # coefficients of 1, i, j, and k -- in that order.
1409 return tr
.coefficient_tuple()[0]
1412 class RealMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1416 The identity function, for embedding real matrices into real
1422 def real_unembed(M
):
1424 The identity function, for unembedding real matrices from real
1430 class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra
,
1431 ConcreteEuclideanJordanAlgebra
):
1433 The rank-n simple EJA consisting of real symmetric n-by-n
1434 matrices, the usual symmetric Jordan product, and the trace inner
1435 product. It has dimension `(n^2 + n)/2` over the reals.
1439 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1443 sage: J = RealSymmetricEJA(2)
1444 sage: e0, e1, e2 = J.gens()
1452 In theory, our "field" can be any subfield of the reals::
1454 sage: RealSymmetricEJA(2, RDF)
1455 Euclidean Jordan algebra of dimension 3 over Real Double Field
1456 sage: RealSymmetricEJA(2, RR)
1457 Euclidean Jordan algebra of dimension 3 over Real Field with
1458 53 bits of precision
1462 The dimension of this algebra is `(n^2 + n) / 2`::
1464 sage: set_random_seed()
1465 sage: n_max = RealSymmetricEJA._max_random_instance_size()
1466 sage: n = ZZ.random_element(1, n_max)
1467 sage: J = RealSymmetricEJA(n)
1468 sage: J.dimension() == (n^2 + n)/2
1471 The Jordan multiplication is what we think it is::
1473 sage: set_random_seed()
1474 sage: J = RealSymmetricEJA.random_instance()
1475 sage: x,y = J.random_elements(2)
1476 sage: actual = (x*y).to_matrix()
1477 sage: X = x.to_matrix()
1478 sage: Y = y.to_matrix()
1479 sage: expected = (X*Y + Y*X)/2
1480 sage: actual == expected
1482 sage: J(expected) == x*y
1485 We can change the generator prefix::
1487 sage: RealSymmetricEJA(3, prefix='q').gens()
1488 (q0, q1, q2, q3, q4, q5)
1490 We can construct the (trivial) algebra of rank zero::
1492 sage: RealSymmetricEJA(0)
1493 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1497 def _denormalized_basis(cls
, n
, field
):
1499 Return a basis for the space of real symmetric n-by-n matrices.
1503 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1507 sage: set_random_seed()
1508 sage: n = ZZ.random_element(1,5)
1509 sage: B = RealSymmetricEJA._denormalized_basis(n,QQ)
1510 sage: all( M.is_symmetric() for M in B)
1514 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1518 for j
in range(i
+1):
1519 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
1523 Sij
= Eij
+ Eij
.transpose()
1529 def _max_random_instance_size():
1530 return 4 # Dimension 10
1533 def random_instance(cls
, field
=AA
, **kwargs
):
1535 Return a random instance of this type of algebra.
1537 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1538 return cls(n
, field
, **kwargs
)
1540 def __init__(self
, n
, field
=AA
, **kwargs
):
1541 basis
= self
._denormalized
_basis
(n
, field
)
1542 super(RealSymmetricEJA
, self
).__init
__(field
,
1546 self
.rank
.set_cache(n
)
1549 class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1553 Embed the n-by-n complex matrix ``M`` into the space of real
1554 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1555 bi` to the block matrix ``[[a,b],[-b,a]]``.
1559 sage: from mjo.eja.eja_algebra import \
1560 ....: ComplexMatrixEuclideanJordanAlgebra
1564 sage: F = QuadraticField(-1, 'I')
1565 sage: x1 = F(4 - 2*i)
1566 sage: x2 = F(1 + 2*i)
1569 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1570 sage: ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1579 Embedding is a homomorphism (isomorphism, in fact)::
1581 sage: set_random_seed()
1582 sage: n = ZZ.random_element(3)
1583 sage: F = QuadraticField(-1, 'I')
1584 sage: X = random_matrix(F, n)
1585 sage: Y = random_matrix(F, n)
1586 sage: Xe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X)
1587 sage: Ye = ComplexMatrixEuclideanJordanAlgebra.real_embed(Y)
1588 sage: XYe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1595 raise ValueError("the matrix 'M' must be square")
1597 # We don't need any adjoined elements...
1598 field
= M
.base_ring().base_ring()
1602 a
= z
.list()[0] # real part, I guess
1603 b
= z
.list()[1] # imag part, I guess
1604 blocks
.append(matrix(field
, 2, [[a
,b
],[-b
,a
]]))
1606 return matrix
.block(field
, n
, blocks
)
1610 def real_unembed(M
):
1612 The inverse of _embed_complex_matrix().
1616 sage: from mjo.eja.eja_algebra import \
1617 ....: ComplexMatrixEuclideanJordanAlgebra
1621 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1622 ....: [-2, 1, -4, 3],
1623 ....: [ 9, 10, 11, 12],
1624 ....: [-10, 9, -12, 11] ])
1625 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(A)
1627 [ 10*I + 9 12*I + 11]
1631 Unembedding is the inverse of embedding::
1633 sage: set_random_seed()
1634 sage: F = QuadraticField(-1, 'I')
1635 sage: M = random_matrix(F, 3)
1636 sage: Me = ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1637 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1643 raise ValueError("the matrix 'M' must be square")
1644 if not n
.mod(2).is_zero():
1645 raise ValueError("the matrix 'M' must be a complex embedding")
1647 # If "M" was normalized, its base ring might have roots
1648 # adjoined and they can stick around after unembedding.
1649 field
= M
.base_ring()
1650 R
= PolynomialRing(field
, 'z')
1653 # Sage doesn't know how to embed AA into QQbar, i.e. how
1654 # to adjoin sqrt(-1) to AA.
1657 F
= field
.extension(z
**2 + 1, 'I', embedding
=CLF(-1).sqrt())
1660 # Go top-left to bottom-right (reading order), converting every
1661 # 2-by-2 block we see to a single complex element.
1663 for k
in range(n
/2):
1664 for j
in range(n
/2):
1665 submat
= M
[2*k
:2*k
+2,2*j
:2*j
+2]
1666 if submat
[0,0] != submat
[1,1]:
1667 raise ValueError('bad on-diagonal submatrix')
1668 if submat
[0,1] != -submat
[1,0]:
1669 raise ValueError('bad off-diagonal submatrix')
1670 z
= submat
[0,0] + submat
[0,1]*i
1673 return matrix(F
, n
/2, elements
)
1677 def matrix_inner_product(cls
,X
,Y
):
1679 Compute a matrix inner product in this algebra directly from
1684 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1688 This gives the same answer as the slow, default method implemented
1689 in :class:`MatrixEuclideanJordanAlgebra`::
1691 sage: set_random_seed()
1692 sage: J = ComplexHermitianEJA.random_instance()
1693 sage: x,y = J.random_elements(2)
1694 sage: Xe = x.to_matrix()
1695 sage: Ye = y.to_matrix()
1696 sage: X = ComplexHermitianEJA.real_unembed(Xe)
1697 sage: Y = ComplexHermitianEJA.real_unembed(Ye)
1698 sage: expected = (X*Y).trace().real()
1699 sage: actual = ComplexHermitianEJA.matrix_inner_product(Xe,Ye)
1700 sage: actual == expected
1704 return RealMatrixEuclideanJordanAlgebra
.matrix_inner_product(X
,Y
)/2
1707 class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra
,
1708 ConcreteEuclideanJordanAlgebra
):
1710 The rank-n simple EJA consisting of complex Hermitian n-by-n
1711 matrices over the real numbers, the usual symmetric Jordan product,
1712 and the real-part-of-trace inner product. It has dimension `n^2` over
1717 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1721 In theory, our "field" can be any subfield of the reals::
1723 sage: ComplexHermitianEJA(2, RDF)
1724 Euclidean Jordan algebra of dimension 4 over Real Double Field
1725 sage: ComplexHermitianEJA(2, RR)
1726 Euclidean Jordan algebra of dimension 4 over Real Field with
1727 53 bits of precision
1731 The dimension of this algebra is `n^2`::
1733 sage: set_random_seed()
1734 sage: n_max = ComplexHermitianEJA._max_random_instance_size()
1735 sage: n = ZZ.random_element(1, n_max)
1736 sage: J = ComplexHermitianEJA(n)
1737 sage: J.dimension() == n^2
1740 The Jordan multiplication is what we think it is::
1742 sage: set_random_seed()
1743 sage: J = ComplexHermitianEJA.random_instance()
1744 sage: x,y = J.random_elements(2)
1745 sage: actual = (x*y).to_matrix()
1746 sage: X = x.to_matrix()
1747 sage: Y = y.to_matrix()
1748 sage: expected = (X*Y + Y*X)/2
1749 sage: actual == expected
1751 sage: J(expected) == x*y
1754 We can change the generator prefix::
1756 sage: ComplexHermitianEJA(2, prefix='z').gens()
1759 We can construct the (trivial) algebra of rank zero::
1761 sage: ComplexHermitianEJA(0)
1762 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1767 def _denormalized_basis(cls
, n
, field
):
1769 Returns a basis for the space of complex Hermitian n-by-n matrices.
1771 Why do we embed these? Basically, because all of numerical linear
1772 algebra assumes that you're working with vectors consisting of `n`
1773 entries from a field and scalars from the same field. There's no way
1774 to tell SageMath that (for example) the vectors contain complex
1775 numbers, while the scalar field is real.
1779 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1783 sage: set_random_seed()
1784 sage: n = ZZ.random_element(1,5)
1785 sage: field = QuadraticField(2, 'sqrt2')
1786 sage: B = ComplexHermitianEJA._denormalized_basis(n, field)
1787 sage: all( M.is_symmetric() for M in B)
1791 R
= PolynomialRing(field
, 'z')
1793 F
= field
.extension(z
**2 + 1, 'I')
1796 # This is like the symmetric case, but we need to be careful:
1798 # * We want conjugate-symmetry, not just symmetry.
1799 # * The diagonal will (as a result) be real.
1803 for j
in range(i
+1):
1804 Eij
= matrix(F
, n
, lambda k
,l
: k
==i
and l
==j
)
1806 Sij
= cls
.real_embed(Eij
)
1809 # The second one has a minus because it's conjugated.
1810 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
1812 Sij_imag
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
1815 # Since we embedded these, we can drop back to the "field" that we
1816 # started with instead of the complex extension "F".
1817 return ( s
.change_ring(field
) for s
in S
)
1820 def __init__(self
, n
, field
=AA
, **kwargs
):
1821 basis
= self
._denormalized
_basis
(n
,field
)
1822 super(ComplexHermitianEJA
,self
).__init
__(field
,
1826 self
.rank
.set_cache(n
)
1829 def _max_random_instance_size():
1830 return 3 # Dimension 9
1833 def random_instance(cls
, field
=AA
, **kwargs
):
1835 Return a random instance of this type of algebra.
1837 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1838 return cls(n
, field
, **kwargs
)
1840 class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1844 Embed the n-by-n quaternion matrix ``M`` into the space of real
1845 matrices of size 4n-by-4n by first sending each quaternion entry `z
1846 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
1847 c+di],[-c + di, a-bi]]`, and then embedding those into a real
1852 sage: from mjo.eja.eja_algebra import \
1853 ....: QuaternionMatrixEuclideanJordanAlgebra
1857 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1858 sage: i,j,k = Q.gens()
1859 sage: x = 1 + 2*i + 3*j + 4*k
1860 sage: M = matrix(Q, 1, [[x]])
1861 sage: QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1867 Embedding is a homomorphism (isomorphism, in fact)::
1869 sage: set_random_seed()
1870 sage: n = ZZ.random_element(2)
1871 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1872 sage: X = random_matrix(Q, n)
1873 sage: Y = random_matrix(Q, n)
1874 sage: Xe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X)
1875 sage: Ye = QuaternionMatrixEuclideanJordanAlgebra.real_embed(Y)
1876 sage: XYe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1881 quaternions
= M
.base_ring()
1884 raise ValueError("the matrix 'M' must be square")
1886 F
= QuadraticField(-1, 'I')
1891 t
= z
.coefficient_tuple()
1896 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
1897 [-c
+ d
*i
, a
- b
*i
]])
1898 realM
= ComplexMatrixEuclideanJordanAlgebra
.real_embed(cplxM
)
1899 blocks
.append(realM
)
1901 # We should have real entries by now, so use the realest field
1902 # we've got for the return value.
1903 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
1908 def real_unembed(M
):
1910 The inverse of _embed_quaternion_matrix().
1914 sage: from mjo.eja.eja_algebra import \
1915 ....: QuaternionMatrixEuclideanJordanAlgebra
1919 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
1920 ....: [-2, 1, -4, 3],
1921 ....: [-3, 4, 1, -2],
1922 ....: [-4, -3, 2, 1]])
1923 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(M)
1924 [1 + 2*i + 3*j + 4*k]
1928 Unembedding is the inverse of embedding::
1930 sage: set_random_seed()
1931 sage: Q = QuaternionAlgebra(QQ, -1, -1)
1932 sage: M = random_matrix(Q, 3)
1933 sage: Me = QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1934 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1940 raise ValueError("the matrix 'M' must be square")
1941 if not n
.mod(4).is_zero():
1942 raise ValueError("the matrix 'M' must be a quaternion embedding")
1944 # Use the base ring of the matrix to ensure that its entries can be
1945 # multiplied by elements of the quaternion algebra.
1946 field
= M
.base_ring()
1947 Q
= QuaternionAlgebra(field
,-1,-1)
1950 # Go top-left to bottom-right (reading order), converting every
1951 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
1954 for l
in range(n
/4):
1955 for m
in range(n
/4):
1956 submat
= ComplexMatrixEuclideanJordanAlgebra
.real_unembed(
1957 M
[4*l
:4*l
+4,4*m
:4*m
+4] )
1958 if submat
[0,0] != submat
[1,1].conjugate():
1959 raise ValueError('bad on-diagonal submatrix')
1960 if submat
[0,1] != -submat
[1,0].conjugate():
1961 raise ValueError('bad off-diagonal submatrix')
1962 z
= submat
[0,0].real()
1963 z
+= submat
[0,0].imag()*i
1964 z
+= submat
[0,1].real()*j
1965 z
+= submat
[0,1].imag()*k
1968 return matrix(Q
, n
/4, elements
)
1972 def matrix_inner_product(cls
,X
,Y
):
1974 Compute a matrix inner product in this algebra directly from
1979 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1983 This gives the same answer as the slow, default method implemented
1984 in :class:`MatrixEuclideanJordanAlgebra`::
1986 sage: set_random_seed()
1987 sage: J = QuaternionHermitianEJA.random_instance()
1988 sage: x,y = J.random_elements(2)
1989 sage: Xe = x.to_matrix()
1990 sage: Ye = y.to_matrix()
1991 sage: X = QuaternionHermitianEJA.real_unembed(Xe)
1992 sage: Y = QuaternionHermitianEJA.real_unembed(Ye)
1993 sage: expected = (X*Y).trace().coefficient_tuple()[0]
1994 sage: actual = QuaternionHermitianEJA.matrix_inner_product(Xe,Ye)
1995 sage: actual == expected
1999 return RealMatrixEuclideanJordanAlgebra
.matrix_inner_product(X
,Y
)/4
2002 class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra
,
2003 ConcreteEuclideanJordanAlgebra
):
2005 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
2006 matrices, the usual symmetric Jordan product, and the
2007 real-part-of-trace inner product. It has dimension `2n^2 - n` over
2012 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2016 In theory, our "field" can be any subfield of the reals::
2018 sage: QuaternionHermitianEJA(2, RDF)
2019 Euclidean Jordan algebra of dimension 6 over Real Double Field
2020 sage: QuaternionHermitianEJA(2, RR)
2021 Euclidean Jordan algebra of dimension 6 over Real Field with
2022 53 bits of precision
2026 The dimension of this algebra is `2*n^2 - n`::
2028 sage: set_random_seed()
2029 sage: n_max = QuaternionHermitianEJA._max_random_instance_size()
2030 sage: n = ZZ.random_element(1, n_max)
2031 sage: J = QuaternionHermitianEJA(n)
2032 sage: J.dimension() == 2*(n^2) - n
2035 The Jordan multiplication is what we think it is::
2037 sage: set_random_seed()
2038 sage: J = QuaternionHermitianEJA.random_instance()
2039 sage: x,y = J.random_elements(2)
2040 sage: actual = (x*y).to_matrix()
2041 sage: X = x.to_matrix()
2042 sage: Y = y.to_matrix()
2043 sage: expected = (X*Y + Y*X)/2
2044 sage: actual == expected
2046 sage: J(expected) == x*y
2049 We can change the generator prefix::
2051 sage: QuaternionHermitianEJA(2, prefix='a').gens()
2052 (a0, a1, a2, a3, a4, a5)
2054 We can construct the (trivial) algebra of rank zero::
2056 sage: QuaternionHermitianEJA(0)
2057 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2061 def _denormalized_basis(cls
, n
, field
):
2063 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
2065 Why do we embed these? Basically, because all of numerical
2066 linear algebra assumes that you're working with vectors consisting
2067 of `n` entries from a field and scalars from the same field. There's
2068 no way to tell SageMath that (for example) the vectors contain
2069 complex numbers, while the scalar field is real.
2073 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2077 sage: set_random_seed()
2078 sage: n = ZZ.random_element(1,5)
2079 sage: B = QuaternionHermitianEJA._denormalized_basis(n,QQ)
2080 sage: all( M.is_symmetric() for M in B )
2084 Q
= QuaternionAlgebra(QQ
,-1,-1)
2087 # This is like the symmetric case, but we need to be careful:
2089 # * We want conjugate-symmetry, not just symmetry.
2090 # * The diagonal will (as a result) be real.
2094 for j
in range(i
+1):
2095 Eij
= matrix(Q
, n
, lambda k
,l
: k
==i
and l
==j
)
2097 Sij
= cls
.real_embed(Eij
)
2100 # The second, third, and fourth ones have a minus
2101 # because they're conjugated.
2102 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
2104 Sij_I
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
2106 Sij_J
= cls
.real_embed(J
*Eij
- J
*Eij
.transpose())
2108 Sij_K
= cls
.real_embed(K
*Eij
- K
*Eij
.transpose())
2111 # Since we embedded these, we can drop back to the "field" that we
2112 # started with instead of the quaternion algebra "Q".
2113 return ( s
.change_ring(field
) for s
in S
)
2116 def __init__(self
, n
, field
=AA
, **kwargs
):
2117 basis
= self
._denormalized
_basis
(n
,field
)
2118 super(QuaternionHermitianEJA
,self
).__init
__(field
,
2122 self
.rank
.set_cache(n
)
2125 def _max_random_instance_size():
2127 The maximum rank of a random QuaternionHermitianEJA.
2129 return 2 # Dimension 6
2132 def random_instance(cls
, field
=AA
, **kwargs
):
2134 Return a random instance of this type of algebra.
2136 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2137 return cls(n
, field
, **kwargs
)
2140 class HadamardEJA(RationalBasisEuclideanJordanAlgebraNg
,
2141 ConcreteEuclideanJordanAlgebra
):
2143 Return the Euclidean Jordan Algebra corresponding to the set
2144 `R^n` under the Hadamard product.
2146 Note: this is nothing more than the Cartesian product of ``n``
2147 copies of the spin algebra. Once Cartesian product algebras
2148 are implemented, this can go.
2152 sage: from mjo.eja.eja_algebra import HadamardEJA
2156 This multiplication table can be verified by hand::
2158 sage: J = HadamardEJA(3)
2159 sage: e0,e1,e2 = J.gens()
2175 We can change the generator prefix::
2177 sage: HadamardEJA(3, prefix='r').gens()
2181 def __init__(self
, n
, field
=AA
, **kwargs
):
2182 V
= VectorSpace(field
, n
)
2185 def jordan_product(x
,y
):
2186 return V([ xi
*yi
for (xi
,yi
) in zip(x
,y
) ])
2187 def inner_product(x
,y
):
2188 return x
.inner_product(y
)
2190 super(HadamardEJA
, self
).__init
__(field
,
2195 self
.rank
.set_cache(n
)
2198 self
.one
.set_cache( self
.zero() )
2200 self
.one
.set_cache( sum(self
.gens()) )
2203 def _max_random_instance_size():
2205 The maximum dimension of a random HadamardEJA.
2210 def random_instance(cls
, field
=AA
, **kwargs
):
2212 Return a random instance of this type of algebra.
2214 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2215 return cls(n
, field
, **kwargs
)
2218 class BilinearFormEJA(RationalBasisEuclideanJordanAlgebraNg
,
2219 ConcreteEuclideanJordanAlgebra
):
2221 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2222 with the half-trace inner product and jordan product ``x*y =
2223 (<Bx,y>,y_bar>, x0*y_bar + y0*x_bar)`` where `B = 1 \times B22` is
2224 a symmetric positive-definite "bilinear form" matrix. Its
2225 dimension is the size of `B`, and it has rank two in dimensions
2226 larger than two. It reduces to the ``JordanSpinEJA`` when `B` is
2227 the identity matrix of order ``n``.
2229 We insist that the one-by-one upper-left identity block of `B` be
2230 passed in as well so that we can be passed a matrix of size zero
2231 to construct a trivial algebra.
2235 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
2236 ....: JordanSpinEJA)
2240 When no bilinear form is specified, the identity matrix is used,
2241 and the resulting algebra is the Jordan spin algebra::
2243 sage: B = matrix.identity(AA,3)
2244 sage: J0 = BilinearFormEJA(B)
2245 sage: J1 = JordanSpinEJA(3)
2246 sage: J0.multiplication_table() == J0.multiplication_table()
2249 An error is raised if the matrix `B` does not correspond to a
2250 positive-definite bilinear form::
2252 sage: B = matrix.random(QQ,2,3)
2253 sage: J = BilinearFormEJA(B)
2254 Traceback (most recent call last):
2256 ValueError: bilinear form is not positive-definite
2257 sage: B = matrix.zero(QQ,3)
2258 sage: J = BilinearFormEJA(B)
2259 Traceback (most recent call last):
2261 ValueError: bilinear form is not positive-definite
2265 We can create a zero-dimensional algebra::
2267 sage: B = matrix.identity(AA,0)
2268 sage: J = BilinearFormEJA(B)
2272 We can check the multiplication condition given in the Jordan, von
2273 Neumann, and Wigner paper (and also discussed on my "On the
2274 symmetry..." paper). Note that this relies heavily on the standard
2275 choice of basis, as does anything utilizing the bilinear form matrix::
2277 sage: set_random_seed()
2278 sage: n = ZZ.random_element(5)
2279 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2280 sage: B11 = matrix.identity(QQ,1)
2281 sage: B22 = M.transpose()*M
2282 sage: B = block_matrix(2,2,[ [B11,0 ],
2284 sage: J = BilinearFormEJA(B)
2285 sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
2286 sage: V = J.vector_space()
2287 sage: sis = [ J.from_vector(V([0] + (M.inverse()*ei).list()))
2288 ....: for ei in eis ]
2289 sage: actual = [ sis[i]*sis[j]
2290 ....: for i in range(n-1)
2291 ....: for j in range(n-1) ]
2292 sage: expected = [ J.one() if i == j else J.zero()
2293 ....: for i in range(n-1)
2294 ....: for j in range(n-1) ]
2295 sage: actual == expected
2298 def __init__(self
, B
, field
=AA
, **kwargs
):
2299 if not B
.is_positive_definite():
2300 raise ValueError("bilinear form is not positive-definite")
2303 V
= VectorSpace(field
, n
)
2305 def inner_product(x
,y
):
2306 return (B
*x
).inner_product(y
)
2308 def jordan_product(x
,y
):
2313 z0
= inner_product(x
,y
)
2314 zbar
= y0
*xbar
+ x0
*ybar
2315 return V([z0
] + zbar
.list())
2317 super(BilinearFormEJA
, self
).__init
__(field
,
2323 # The rank of this algebra is two, unless we're in a
2324 # one-dimensional ambient space (because the rank is bounded
2325 # by the ambient dimension).
2326 self
.rank
.set_cache(min(n
,2))
2329 self
.one
.set_cache( self
.zero() )
2331 self
.one
.set_cache( self
.monomial(0) )
2334 def _max_random_instance_size():
2336 The maximum dimension of a random BilinearFormEJA.
2341 def random_instance(cls
, field
=AA
, **kwargs
):
2343 Return a random instance of this algebra.
2345 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2347 B
= matrix
.identity(field
, n
)
2348 return cls(B
, field
, **kwargs
)
2350 B11
= matrix
.identity(field
,1)
2351 M
= matrix
.random(field
, n
-1)
2352 I
= matrix
.identity(field
, n
-1)
2353 alpha
= field
.zero()
2354 while alpha
.is_zero():
2355 alpha
= field
.random_element().abs()
2356 B22
= M
.transpose()*M
+ alpha
*I
2358 from sage
.matrix
.special
import block_matrix
2359 B
= block_matrix(2,2, [ [B11
, ZZ(0) ],
2362 return cls(B
, field
, **kwargs
)
2365 class JordanSpinEJA(BilinearFormEJA
):
2367 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2368 with the usual inner product and jordan product ``x*y =
2369 (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2374 sage: from mjo.eja.eja_algebra import JordanSpinEJA
2378 This multiplication table can be verified by hand::
2380 sage: J = JordanSpinEJA(4)
2381 sage: e0,e1,e2,e3 = J.gens()
2397 We can change the generator prefix::
2399 sage: JordanSpinEJA(2, prefix='B').gens()
2404 Ensure that we have the usual inner product on `R^n`::
2406 sage: set_random_seed()
2407 sage: J = JordanSpinEJA.random_instance()
2408 sage: x,y = J.random_elements(2)
2409 sage: actual = x.inner_product(y)
2410 sage: expected = x.to_vector().inner_product(y.to_vector())
2411 sage: actual == expected
2415 def __init__(self
, n
, field
=AA
, **kwargs
):
2416 # This is a special case of the BilinearFormEJA with the identity
2417 # matrix as its bilinear form.
2418 B
= matrix
.identity(field
, n
)
2419 super(JordanSpinEJA
, self
).__init
__(B
, field
, **kwargs
)
2422 def _max_random_instance_size():
2424 The maximum dimension of a random JordanSpinEJA.
2429 def random_instance(cls
, field
=AA
, **kwargs
):
2431 Return a random instance of this type of algebra.
2433 Needed here to override the implementation for ``BilinearFormEJA``.
2435 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2436 return cls(n
, field
, **kwargs
)
2439 class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra
,
2440 ConcreteEuclideanJordanAlgebra
):
2442 The trivial Euclidean Jordan algebra consisting of only a zero element.
2446 sage: from mjo.eja.eja_algebra import TrivialEJA
2450 sage: J = TrivialEJA()
2457 sage: 7*J.one()*12*J.one()
2459 sage: J.one().inner_product(J.one())
2461 sage: J.one().norm()
2463 sage: J.one().subalgebra_generated_by()
2464 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2469 def __init__(self
, field
=AA
, **kwargs
):
2471 self
._inner
_product
_matrix
= matrix(field
,0)
2472 super(TrivialEJA
, self
).__init
__(field
,
2476 # The rank is zero using my definition, namely the dimension of the
2477 # largest subalgebra generated by any element.
2478 self
.rank
.set_cache(0)
2479 self
.one
.set_cache( self
.zero() )
2482 def random_instance(cls
, field
=AA
, **kwargs
):
2483 # We don't take a "size" argument so the superclass method is
2484 # inappropriate for us.
2485 return cls(field
, **kwargs
)
2487 class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra
):
2489 The external (orthogonal) direct sum of two other Euclidean Jordan
2490 algebras. Essentially the Cartesian product of its two factors.
2491 Every Euclidean Jordan algebra decomposes into an orthogonal
2492 direct sum of simple Euclidean Jordan algebras, so no generality
2493 is lost by providing only this construction.
2497 sage: from mjo.eja.eja_algebra import (random_eja,
2499 ....: RealSymmetricEJA,
2504 sage: J1 = HadamardEJA(2)
2505 sage: J2 = RealSymmetricEJA(3)
2506 sage: J = DirectSumEJA(J1,J2)
2514 The external direct sum construction is only valid when the two factors
2515 have the same base ring; an error is raised otherwise::
2517 sage: set_random_seed()
2518 sage: J1 = random_eja(AA)
2519 sage: J2 = random_eja(QQ)
2520 sage: J = DirectSumEJA(J1,J2)
2521 Traceback (most recent call last):
2523 ValueError: algebras must share the same base field
2526 def __init__(self
, J1
, J2
, **kwargs
):
2527 if J1
.base_ring() != J2
.base_ring():
2528 raise ValueError("algebras must share the same base field")
2529 field
= J1
.base_ring()
2531 self
._factors
= (J1
, J2
)
2535 V
= VectorSpace(field
, n
)
2536 mult_table
= [ [ V
.zero() for j
in range(n
) ]
2540 p
= (J1
.monomial(i
)*J1
.monomial(j
)).to_vector()
2541 mult_table
[i
][j
] = V(p
.list() + [field
.zero()]*n2
)
2545 p
= (J2
.monomial(i
)*J2
.monomial(j
)).to_vector()
2546 mult_table
[n1
+i
][n1
+j
] = V([field
.zero()]*n1
+ p
.list())
2548 super(DirectSumEJA
, self
).__init
__(field
,
2552 self
.rank
.set_cache(J1
.rank() + J2
.rank())
2557 Return the pair of this algebra's factors.
2561 sage: from mjo.eja.eja_algebra import (HadamardEJA,
2562 ....: JordanSpinEJA,
2567 sage: J1 = HadamardEJA(2,QQ)
2568 sage: J2 = JordanSpinEJA(3,QQ)
2569 sage: J = DirectSumEJA(J1,J2)
2571 (Euclidean Jordan algebra of dimension 2 over Rational Field,
2572 Euclidean Jordan algebra of dimension 3 over Rational Field)
2575 return self
._factors
2577 def projections(self
):
2579 Return a pair of projections onto this algebra's factors.
2583 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
2584 ....: ComplexHermitianEJA,
2589 sage: J1 = JordanSpinEJA(2)
2590 sage: J2 = ComplexHermitianEJA(2)
2591 sage: J = DirectSumEJA(J1,J2)
2592 sage: (pi_left, pi_right) = J.projections()
2593 sage: J.one().to_vector()
2595 sage: pi_left(J.one()).to_vector()
2597 sage: pi_right(J.one()).to_vector()
2601 (J1
,J2
) = self
.factors()
2604 V_basis
= self
.vector_space().basis()
2605 # Need to specify the dimensions explicitly so that we don't
2606 # wind up with a zero-by-zero matrix when we want e.g. a
2607 # zero-by-two matrix (important for composing things).
2608 P1
= matrix(self
.base_ring(), m
, m
+n
, V_basis
[:m
])
2609 P2
= matrix(self
.base_ring(), n
, m
+n
, V_basis
[m
:])
2610 pi_left
= FiniteDimensionalEuclideanJordanAlgebraOperator(self
,J1
,P1
)
2611 pi_right
= FiniteDimensionalEuclideanJordanAlgebraOperator(self
,J2
,P2
)
2612 return (pi_left
, pi_right
)
2614 def inclusions(self
):
2616 Return the pair of inclusion maps from our factors into us.
2620 sage: from mjo.eja.eja_algebra import (random_eja,
2621 ....: JordanSpinEJA,
2622 ....: RealSymmetricEJA,
2627 sage: J1 = JordanSpinEJA(3)
2628 sage: J2 = RealSymmetricEJA(2)
2629 sage: J = DirectSumEJA(J1,J2)
2630 sage: (iota_left, iota_right) = J.inclusions()
2631 sage: iota_left(J1.zero()) == J.zero()
2633 sage: iota_right(J2.zero()) == J.zero()
2635 sage: J1.one().to_vector()
2637 sage: iota_left(J1.one()).to_vector()
2639 sage: J2.one().to_vector()
2641 sage: iota_right(J2.one()).to_vector()
2643 sage: J.one().to_vector()
2648 Composing a projection with the corresponding inclusion should
2649 produce the identity map, and mismatching them should produce
2652 sage: set_random_seed()
2653 sage: J1 = random_eja()
2654 sage: J2 = random_eja()
2655 sage: J = DirectSumEJA(J1,J2)
2656 sage: (iota_left, iota_right) = J.inclusions()
2657 sage: (pi_left, pi_right) = J.projections()
2658 sage: pi_left*iota_left == J1.one().operator()
2660 sage: pi_right*iota_right == J2.one().operator()
2662 sage: (pi_left*iota_right).is_zero()
2664 sage: (pi_right*iota_left).is_zero()
2668 (J1
,J2
) = self
.factors()
2671 V_basis
= self
.vector_space().basis()
2672 # Need to specify the dimensions explicitly so that we don't
2673 # wind up with a zero-by-zero matrix when we want e.g. a
2674 # two-by-zero matrix (important for composing things).
2675 I1
= matrix
.column(self
.base_ring(), m
, m
+n
, V_basis
[:m
])
2676 I2
= matrix
.column(self
.base_ring(), n
, m
+n
, V_basis
[m
:])
2677 iota_left
= FiniteDimensionalEuclideanJordanAlgebraOperator(J1
,self
,I1
)
2678 iota_right
= FiniteDimensionalEuclideanJordanAlgebraOperator(J2
,self
,I2
)
2679 return (iota_left
, iota_right
)
2681 def inner_product(self
, x
, y
):
2683 The standard Cartesian inner-product.
2685 We project ``x`` and ``y`` onto our factors, and add up the
2686 inner-products from the subalgebras.
2691 sage: from mjo.eja.eja_algebra import (HadamardEJA,
2692 ....: QuaternionHermitianEJA,
2697 sage: J1 = HadamardEJA(3,QQ)
2698 sage: J2 = QuaternionHermitianEJA(2,QQ,normalize_basis=False)
2699 sage: J = DirectSumEJA(J1,J2)
2704 sage: x1.inner_product(x2)
2706 sage: y1.inner_product(y2)
2708 sage: J.one().inner_product(J.one())
2712 (pi_left
, pi_right
) = self
.projections()
2718 return (x1
.inner_product(y1
) + x2
.inner_product(y2
))
2722 random_eja
= ConcreteEuclideanJordanAlgebra
.random_instance