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
.lazy_import
import lazy_import
28 from sage
.misc
.table
import table
29 from sage
.modules
.free_module
import FreeModule
, VectorSpace
30 from sage
.rings
.all
import (ZZ
, QQ
, AA
, QQbar
, RR
, RLF
, CLF
,
33 from mjo
.eja
.eja_element
import FiniteDimensionalEuclideanJordanAlgebraElement
34 lazy_import('mjo.eja.eja_subalgebra',
35 'FiniteDimensionalEuclideanJordanSubalgebra')
36 from mjo
.eja
.eja_operator
import FiniteDimensionalEuclideanJordanAlgebraOperator
37 from mjo
.eja
.eja_utils
import _mat2vec
39 class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule
):
41 def _coerce_map_from_base_ring(self
):
43 Disable the map from the base ring into the algebra.
45 Performing a nonsense conversion like this automatically
46 is counterpedagogical. The fallback is to try the usual
47 element constructor, which should also fail.
51 sage: from mjo.eja.eja_algebra import random_eja
55 sage: set_random_seed()
56 sage: J = random_eja()
58 Traceback (most recent call last):
60 ValueError: not a naturally-represented algebra element
76 sage: from mjo.eja.eja_algebra import (
77 ....: FiniteDimensionalEuclideanJordanAlgebra,
83 By definition, Jordan multiplication commutes::
85 sage: set_random_seed()
86 sage: J = random_eja()
87 sage: x,y = J.random_elements(2)
93 The ``field`` we're given must be real with ``check_field=True``::
95 sage: JordanSpinEJA(2,QQbar)
96 Traceback (most recent call last):
98 ValueError: scalar field is not real
100 The multiplication table must be square with ``check_axioms=True``::
102 sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,((),()))
103 Traceback (most recent call last):
105 ValueError: multiplication table is not square
109 if not field
.is_subring(RR
):
110 # Note: this does return true for the real algebraic
111 # field, the rationals, and any quadratic field where
112 # we've specified a real embedding.
113 raise ValueError("scalar field is not real")
115 # The multiplication table had better be square
118 if not all( len(l
) == n
for l
in mult_table
):
119 raise ValueError("multiplication table is not square")
121 self
._natural
_basis
= natural_basis
124 category
= MagmaticAlgebras(field
).FiniteDimensional()
125 category
= category
.WithBasis().Unital()
127 fda
= super(FiniteDimensionalEuclideanJordanAlgebra
, self
)
132 self
.print_options(bracket
='')
134 # The multiplication table we're given is necessarily in terms
135 # of vectors, because we don't have an algebra yet for
136 # anything to be an element of. However, it's faster in the
137 # long run to have the multiplication table be in terms of
138 # algebra elements. We do this after calling the superclass
139 # constructor so that from_vector() knows what to do.
140 self
._multiplication
_table
= [
141 list(map(lambda x
: self
.from_vector(x
), ls
))
146 if not self
._is
_commutative
():
147 raise ValueError("algebra is not commutative")
148 if not self
._is
_jordanian
():
149 raise ValueError("Jordan identity does not hold")
150 if not self
._inner
_product
_is
_associative
():
151 raise ValueError("inner product is not associative")
153 def _element_constructor_(self
, elt
):
155 Construct an element of this algebra from its natural
158 This gets called only after the parent element _call_ method
159 fails to find a coercion for the argument.
163 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
165 ....: RealSymmetricEJA)
169 The identity in `S^n` is converted to the identity in the EJA::
171 sage: J = RealSymmetricEJA(3)
172 sage: I = matrix.identity(QQ,3)
173 sage: J(I) == J.one()
176 This skew-symmetric matrix can't be represented in the EJA::
178 sage: J = RealSymmetricEJA(3)
179 sage: A = matrix(QQ,3, lambda i,j: i-j)
181 Traceback (most recent call last):
183 ArithmeticError: vector is not in free module
187 Ensure that we can convert any element of the two non-matrix
188 simple algebras (whose natural representations are their usual
189 vector representations) back and forth faithfully::
191 sage: set_random_seed()
192 sage: J = HadamardEJA.random_instance()
193 sage: x = J.random_element()
194 sage: J(x.to_vector().column()) == x
196 sage: J = JordanSpinEJA.random_instance()
197 sage: x = J.random_element()
198 sage: J(x.to_vector().column()) == x
202 msg
= "not a naturally-represented algebra element"
204 # The superclass implementation of random_element()
205 # needs to be able to coerce "0" into the algebra.
207 elif elt
in self
.base_ring():
208 # Ensure that no base ring -> algebra coercion is performed
209 # by this method. There's some stupidity in sage that would
210 # otherwise propagate to this method; for example, sage thinks
211 # that the integer 3 belongs to the space of 2-by-2 matrices.
212 raise ValueError(msg
)
214 natural_basis
= self
.natural_basis()
215 basis_space
= natural_basis
[0].matrix_space()
216 if elt
not in basis_space
:
217 raise ValueError(msg
)
219 # Thanks for nothing! Matrix spaces aren't vector spaces in
220 # Sage, so we have to figure out its natural-basis coordinates
221 # ourselves. We use the basis space's ring instead of the
222 # element's ring because the basis space might be an algebraic
223 # closure whereas the base ring of the 3-by-3 identity matrix
224 # could be QQ instead of QQbar.
225 V
= VectorSpace(basis_space
.base_ring(), elt
.nrows()*elt
.ncols())
226 W
= V
.span_of_basis( _mat2vec(s
) for s
in natural_basis
)
227 coords
= W
.coordinate_vector(_mat2vec(elt
))
228 return self
.from_vector(coords
)
232 Return a string representation of ``self``.
236 sage: from mjo.eja.eja_algebra import JordanSpinEJA
240 Ensure that it says what we think it says::
242 sage: JordanSpinEJA(2, field=AA)
243 Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
244 sage: JordanSpinEJA(3, field=RDF)
245 Euclidean Jordan algebra of dimension 3 over Real Double Field
248 fmt
= "Euclidean Jordan algebra of dimension {} over {}"
249 return fmt
.format(self
.dimension(), self
.base_ring())
251 def product_on_basis(self
, i
, j
):
252 return self
._multiplication
_table
[i
][j
]
254 def _is_commutative(self
):
256 Whether or not this algebra's multiplication table is commutative.
258 This method should of course always return ``True``, unless
259 this algebra was constructed with ``check_axioms=False`` and
260 passed an invalid multiplication table.
262 return all( self
.product_on_basis(i
,j
) == self
.product_on_basis(i
,j
)
263 for i
in range(self
.dimension())
264 for j
in range(self
.dimension()) )
266 def _is_jordanian(self
):
268 Whether or not this algebra's multiplication table respects the
269 Jordan identity `(x^{2})(xy) = x(x^{2}y)`.
271 We only check one arrangement of `x` and `y`, so for a
272 ``True`` result to be truly true, you should also check
273 :meth:`_is_commutative`. This method should of course always
274 return ``True``, unless this algebra was constructed with
275 ``check_axioms=False`` and passed an invalid multiplication table.
277 return all( (self
.monomial(i
)**2)*(self
.monomial(i
)*self
.monomial(j
))
279 (self
.monomial(i
))*((self
.monomial(i
)**2)*self
.monomial(j
))
280 for i
in range(self
.dimension())
281 for j
in range(self
.dimension()) )
283 def _inner_product_is_associative(self
):
285 Return whether or not this algebra's inner product `B` is
286 associative; that is, whether or not `B(xy,z) = B(x,yz)`.
288 This method should of course always return ``True``, unless
289 this algebra was constructed with ``check_axioms=False`` and
290 passed an invalid multiplication table.
293 # Used to check whether or not something is zero in an inexact
294 # ring. This number is sufficient to allow the construction of
295 # QuaternionHermitianEJA(2, RDF) with check_axioms=True.
298 for i
in range(self
.dimension()):
299 for j
in range(self
.dimension()):
300 for k
in range(self
.dimension()):
304 diff
= (x
*y
).inner_product(z
) - x
.inner_product(y
*z
)
306 if self
.base_ring().is_exact():
310 if diff
.abs() > epsilon
:
316 def characteristic_polynomial_of(self
):
318 Return the algebra's "characteristic polynomial of" function,
319 which is itself a multivariate polynomial that, when evaluated
320 at the coordinates of some algebra element, returns that
321 element's characteristic polynomial.
323 The resulting polynomial has `n+1` variables, where `n` is the
324 dimension of this algebra. The first `n` variables correspond to
325 the coordinates of an algebra element: when evaluated at the
326 coordinates of an algebra element with respect to a certain
327 basis, the result is a univariate polynomial (in the one
328 remaining variable ``t``), namely the characteristic polynomial
333 sage: from mjo.eja.eja_algebra import JordanSpinEJA, TrivialEJA
337 The characteristic polynomial in the spin algebra is given in
338 Alizadeh, Example 11.11::
340 sage: J = JordanSpinEJA(3)
341 sage: p = J.characteristic_polynomial_of(); p
342 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
343 sage: xvec = J.one().to_vector()
347 By definition, the characteristic polynomial is a monic
348 degree-zero polynomial in a rank-zero algebra. Note that
349 Cayley-Hamilton is indeed satisfied since the polynomial
350 ``1`` evaluates to the identity element of the algebra on
353 sage: J = TrivialEJA()
354 sage: J.characteristic_polynomial_of()
361 # The list of coefficient polynomials a_0, a_1, a_2, ..., a_(r-1).
362 a
= self
._charpoly
_coefficients
()
364 # We go to a bit of trouble here to reorder the
365 # indeterminates, so that it's easier to evaluate the
366 # characteristic polynomial at x's coordinates and get back
367 # something in terms of t, which is what we want.
368 S
= PolynomialRing(self
.base_ring(),'t')
372 S
= PolynomialRing(S
, R
.variable_names())
375 return (t
**r
+ sum( a
[k
]*(t
**k
) for k
in range(r
) ))
377 def coordinate_polynomial_ring(self
):
379 The multivariate polynomial ring in which this algebra's
380 :meth:`characteristic_polynomial_of` lives.
384 sage: from mjo.eja.eja_algebra import (HadamardEJA,
385 ....: RealSymmetricEJA)
389 sage: J = HadamardEJA(2)
390 sage: J.coordinate_polynomial_ring()
391 Multivariate Polynomial Ring in X1, X2...
392 sage: J = RealSymmetricEJA(3,QQ)
393 sage: J.coordinate_polynomial_ring()
394 Multivariate Polynomial Ring in X1, X2, X3, X4, X5, X6...
397 var_names
= tuple( "X%d" % z
for z
in range(1, self
.dimension()+1) )
398 return PolynomialRing(self
.base_ring(), var_names
)
400 def inner_product(self
, x
, y
):
402 The inner product associated with this Euclidean Jordan algebra.
404 Defaults to the trace inner product, but can be overridden by
405 subclasses if they are sure that the necessary properties are
410 sage: from mjo.eja.eja_algebra import (random_eja,
412 ....: BilinearFormEJA)
416 Our inner product is "associative," which means the following for
417 a symmetric bilinear form::
419 sage: set_random_seed()
420 sage: J = random_eja()
421 sage: x,y,z = J.random_elements(3)
422 sage: (x*y).inner_product(z) == y.inner_product(x*z)
427 Ensure that this is the usual inner product for the algebras
430 sage: set_random_seed()
431 sage: J = HadamardEJA.random_instance()
432 sage: x,y = J.random_elements(2)
433 sage: actual = x.inner_product(y)
434 sage: expected = x.to_vector().inner_product(y.to_vector())
435 sage: actual == expected
438 Ensure that this is one-half of the trace inner-product in a
439 BilinearFormEJA that isn't just the reals (when ``n`` isn't
440 one). This is in Faraut and Koranyi, and also my "On the
443 sage: set_random_seed()
444 sage: J = BilinearFormEJA.random_instance()
445 sage: n = J.dimension()
446 sage: x = J.random_element()
447 sage: y = J.random_element()
448 sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
451 B
= self
._inner
_product
_matrix
452 return (B
*x
.to_vector()).inner_product(y
.to_vector())
455 def is_trivial(self
):
457 Return whether or not this algebra is trivial.
459 A trivial algebra contains only the zero element.
463 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
468 sage: J = ComplexHermitianEJA(3)
474 sage: J = TrivialEJA()
479 return self
.dimension() == 0
482 def multiplication_table(self
):
484 Return a visual representation of this algebra's multiplication
485 table (on basis elements).
489 sage: from mjo.eja.eja_algebra import JordanSpinEJA
493 sage: J = JordanSpinEJA(4)
494 sage: J.multiplication_table()
495 +----++----+----+----+----+
496 | * || e0 | e1 | e2 | e3 |
497 +====++====+====+====+====+
498 | e0 || e0 | e1 | e2 | e3 |
499 +----++----+----+----+----+
500 | e1 || e1 | e0 | 0 | 0 |
501 +----++----+----+----+----+
502 | e2 || e2 | 0 | e0 | 0 |
503 +----++----+----+----+----+
504 | e3 || e3 | 0 | 0 | e0 |
505 +----++----+----+----+----+
508 M
= list(self
._multiplication
_table
) # copy
509 for i
in range(len(M
)):
510 # M had better be "square"
511 M
[i
] = [self
.monomial(i
)] + M
[i
]
512 M
= [["*"] + list(self
.gens())] + M
513 return table(M
, header_row
=True, header_column
=True, frame
=True)
516 def natural_basis(self
):
518 Return a more-natural representation of this algebra's basis.
520 Every finite-dimensional Euclidean Jordan Algebra is a direct
521 sum of five simple algebras, four of which comprise Hermitian
522 matrices. This method returns the original "natural" basis
523 for our underlying vector space. (Typically, the natural basis
524 is used to construct the multiplication table in the first place.)
526 Note that this will always return a matrix. The standard basis
527 in `R^n` will be returned as `n`-by-`1` column matrices.
531 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
532 ....: RealSymmetricEJA)
536 sage: J = RealSymmetricEJA(2)
538 Finite family {0: e0, 1: e1, 2: e2}
539 sage: J.natural_basis()
541 [1 0] [ 0 0.7071067811865475?] [0 0]
542 [0 0], [0.7071067811865475? 0], [0 1]
547 sage: J = JordanSpinEJA(2)
549 Finite family {0: e0, 1: e1}
550 sage: J.natural_basis()
557 if self
._natural
_basis
is None:
558 M
= self
.natural_basis_space()
559 return tuple( M(b
.to_vector()) for b
in self
.basis() )
561 return self
._natural
_basis
564 def natural_basis_space(self
):
566 Return the matrix space in which this algebra's natural basis
569 Generally this will be an `n`-by-`1` column-vector space,
570 except when the algebra is trivial. There it's `n`-by-`n`
571 (where `n` is zero), to ensure that two elements of the
572 natural basis space (empty matrices) can be multiplied.
574 if self
.is_trivial():
575 return MatrixSpace(self
.base_ring(), 0)
576 elif self
._natural
_basis
is None or len(self
._natural
_basis
) == 0:
577 return MatrixSpace(self
.base_ring(), self
.dimension(), 1)
579 return self
._natural
_basis
[0].matrix_space()
585 Return the unit element of this algebra.
589 sage: from mjo.eja.eja_algebra import (HadamardEJA,
594 sage: J = HadamardEJA(5)
596 e0 + e1 + e2 + e3 + e4
600 The identity element acts like the identity::
602 sage: set_random_seed()
603 sage: J = random_eja()
604 sage: x = J.random_element()
605 sage: J.one()*x == x and x*J.one() == x
608 The matrix of the unit element's operator is the identity::
610 sage: set_random_seed()
611 sage: J = random_eja()
612 sage: actual = J.one().operator().matrix()
613 sage: expected = matrix.identity(J.base_ring(), J.dimension())
614 sage: actual == expected
617 Ensure that the cached unit element (often precomputed by
618 hand) agrees with the computed one::
620 sage: set_random_seed()
621 sage: J = random_eja()
622 sage: cached = J.one()
623 sage: J.one.clear_cache()
624 sage: J.one() == cached
628 # We can brute-force compute the matrices of the operators
629 # that correspond to the basis elements of this algebra.
630 # If some linear combination of those basis elements is the
631 # algebra identity, then the same linear combination of
632 # their matrices has to be the identity matrix.
634 # Of course, matrices aren't vectors in sage, so we have to
635 # appeal to the "long vectors" isometry.
636 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
638 # Now we use basic linear algebra to find the coefficients,
639 # of the matrices-as-vectors-linear-combination, which should
640 # work for the original algebra basis too.
641 A
= matrix(self
.base_ring(), oper_vecs
)
643 # We used the isometry on the left-hand side already, but we
644 # still need to do it for the right-hand side. Recall that we
645 # wanted something that summed to the identity matrix.
646 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
648 # Now if there's an identity element in the algebra, this
649 # should work. We solve on the left to avoid having to
650 # transpose the matrix "A".
651 return self
.from_vector(A
.solve_left(b
))
654 def peirce_decomposition(self
, c
):
656 The Peirce decomposition of this algebra relative to the
659 In the future, this can be extended to a complete system of
660 orthogonal idempotents.
664 - ``c`` -- an idempotent of this algebra.
668 A triple (J0, J5, J1) containing two subalgebras and one subspace
671 - ``J0`` -- the algebra on the eigenspace of ``c.operator()``
672 corresponding to the eigenvalue zero.
674 - ``J5`` -- the eigenspace (NOT a subalgebra) of ``c.operator()``
675 corresponding to the eigenvalue one-half.
677 - ``J1`` -- the algebra on the eigenspace of ``c.operator()``
678 corresponding to the eigenvalue one.
680 These are the only possible eigenspaces for that operator, and this
681 algebra is a direct sum of them. The spaces ``J0`` and ``J1`` are
682 orthogonal, and are subalgebras of this algebra with the appropriate
687 sage: from mjo.eja.eja_algebra import random_eja, RealSymmetricEJA
691 The canonical example comes from the symmetric matrices, which
692 decompose into diagonal and off-diagonal parts::
694 sage: J = RealSymmetricEJA(3)
695 sage: C = matrix(QQ, [ [1,0,0],
699 sage: J0,J5,J1 = J.peirce_decomposition(c)
701 Euclidean Jordan algebra of dimension 1...
703 Vector space of degree 6 and dimension 2...
705 Euclidean Jordan algebra of dimension 3...
706 sage: J0.one().natural_representation()
710 sage: orig_df = AA.options.display_format
711 sage: AA.options.display_format = 'radical'
712 sage: J.from_vector(J5.basis()[0]).natural_representation()
716 sage: J.from_vector(J5.basis()[1]).natural_representation()
720 sage: AA.options.display_format = orig_df
721 sage: J1.one().natural_representation()
728 Every algebra decomposes trivially with respect to its identity
731 sage: set_random_seed()
732 sage: J = random_eja()
733 sage: J0,J5,J1 = J.peirce_decomposition(J.one())
734 sage: J0.dimension() == 0 and J5.dimension() == 0
736 sage: J1.superalgebra() == J and J1.dimension() == J.dimension()
739 The decomposition is into eigenspaces, and its components are
740 therefore necessarily orthogonal. Moreover, the identity
741 elements in the two subalgebras are the projections onto their
742 respective subspaces of the superalgebra's identity element::
744 sage: set_random_seed()
745 sage: J = random_eja()
746 sage: x = J.random_element()
747 sage: if not J.is_trivial():
748 ....: while x.is_nilpotent():
749 ....: x = J.random_element()
750 sage: c = x.subalgebra_idempotent()
751 sage: J0,J5,J1 = J.peirce_decomposition(c)
753 sage: for (w,y,z) in zip(J0.basis(), J5.basis(), J1.basis()):
754 ....: w = w.superalgebra_element()
755 ....: y = J.from_vector(y)
756 ....: z = z.superalgebra_element()
757 ....: ipsum += w.inner_product(y).abs()
758 ....: ipsum += w.inner_product(z).abs()
759 ....: ipsum += y.inner_product(z).abs()
762 sage: J1(c) == J1.one()
764 sage: J0(J.one() - c) == J0.one()
768 if not c
.is_idempotent():
769 raise ValueError("element is not idempotent: %s" % c
)
771 # Default these to what they should be if they turn out to be
772 # trivial, because eigenspaces_left() won't return eigenvalues
773 # corresponding to trivial spaces (e.g. it returns only the
774 # eigenspace corresponding to lambda=1 if you take the
775 # decomposition relative to the identity element).
776 trivial
= FiniteDimensionalEuclideanJordanSubalgebra(self
, ())
777 J0
= trivial
# eigenvalue zero
778 J5
= VectorSpace(self
.base_ring(), 0) # eigenvalue one-half
779 J1
= trivial
# eigenvalue one
781 for (eigval
, eigspace
) in c
.operator().matrix().right_eigenspaces():
782 if eigval
== ~
(self
.base_ring()(2)):
785 gens
= tuple( self
.from_vector(b
) for b
in eigspace
.basis() )
786 subalg
= FiniteDimensionalEuclideanJordanSubalgebra(self
,
794 raise ValueError("unexpected eigenvalue: %s" % eigval
)
799 def random_element(self
, thorough
=False):
801 Return a random element of this algebra.
803 Our algebra superclass method only returns a linear
804 combination of at most two basis elements. We instead
805 want the vector space "random element" method that
806 returns a more diverse selection.
810 - ``thorough`` -- (boolean; default False) whether or not we
811 should generate irrational coefficients for the random
812 element when our base ring is irrational; this slows the
813 algebra operations to a crawl, but any truly random method
817 # For a general base ring... maybe we can trust this to do the
818 # right thing? Unlikely, but.
819 V
= self
.vector_space()
820 v
= V
.random_element()
822 if self
.base_ring() is AA
:
823 # The "random element" method of the algebraic reals is
824 # stupid at the moment, and only returns integers between
825 # -2 and 2, inclusive:
827 # https://trac.sagemath.org/ticket/30875
829 # Instead, we implement our own "random vector" method,
830 # and then coerce that into the algebra. We use the vector
831 # space degree here instead of the dimension because a
832 # subalgebra could (for example) be spanned by only two
833 # vectors, each with five coordinates. We need to
834 # generate all five coordinates.
836 v
*= QQbar
.random_element().real()
838 v
*= QQ
.random_element()
840 return self
.from_vector(V
.coordinate_vector(v
))
842 def random_elements(self
, count
, thorough
=False):
844 Return ``count`` random elements as a tuple.
848 - ``thorough`` -- (boolean; default False) whether or not we
849 should generate irrational coefficients for the random
850 elements when our base ring is irrational; this slows the
851 algebra operations to a crawl, but any truly random method
856 sage: from mjo.eja.eja_algebra import JordanSpinEJA
860 sage: J = JordanSpinEJA(3)
861 sage: x,y,z = J.random_elements(3)
862 sage: all( [ x in J, y in J, z in J ])
864 sage: len( J.random_elements(10) ) == 10
868 return tuple( self
.random_element(thorough
)
869 for idx
in range(count
) )
873 def _charpoly_coefficients(self
):
875 The `r` polynomial coefficients of the "characteristic polynomial
879 R
= self
.coordinate_polynomial_ring()
881 F
= R
.fraction_field()
884 # From a result in my book, these are the entries of the
885 # basis representation of L_x.
886 return sum( vars[k
]*self
.monomial(k
).operator().matrix()[i
,j
]
889 L_x
= matrix(F
, n
, n
, L_x_i_j
)
892 if self
.rank
.is_in_cache():
894 # There's no need to pad the system with redundant
895 # columns if we *know* they'll be redundant.
898 # Compute an extra power in case the rank is equal to
899 # the dimension (otherwise, we would stop at x^(r-1)).
900 x_powers
= [ (L_x
**k
)*self
.one().to_vector()
901 for k
in range(n
+1) ]
902 A
= matrix
.column(F
, x_powers
[:n
])
903 AE
= A
.extended_echelon_form()
910 # The theory says that only the first "r" coefficients are
911 # nonzero, and they actually live in the original polynomial
912 # ring and not the fraction field. We negate them because
913 # in the actual characteristic polynomial, they get moved
914 # to the other side where x^r lives.
915 return -A_rref
.solve_right(E
*b
).change_ring(R
)[:r
]
920 Return the rank of this EJA.
922 This is a cached method because we know the rank a priori for
923 all of the algebras we can construct. Thus we can avoid the
924 expensive ``_charpoly_coefficients()`` call unless we truly
925 need to compute the whole characteristic polynomial.
929 sage: from mjo.eja.eja_algebra import (HadamardEJA,
931 ....: RealSymmetricEJA,
932 ....: ComplexHermitianEJA,
933 ....: QuaternionHermitianEJA,
938 The rank of the Jordan spin algebra is always two::
940 sage: JordanSpinEJA(2).rank()
942 sage: JordanSpinEJA(3).rank()
944 sage: JordanSpinEJA(4).rank()
947 The rank of the `n`-by-`n` Hermitian real, complex, or
948 quaternion matrices is `n`::
950 sage: RealSymmetricEJA(4).rank()
952 sage: ComplexHermitianEJA(3).rank()
954 sage: QuaternionHermitianEJA(2).rank()
959 Ensure that every EJA that we know how to construct has a
960 positive integer rank, unless the algebra is trivial in
961 which case its rank will be zero::
963 sage: set_random_seed()
964 sage: J = random_eja()
968 sage: r > 0 or (r == 0 and J.is_trivial())
971 Ensure that computing the rank actually works, since the ranks
972 of all simple algebras are known and will be cached by default::
974 sage: J = HadamardEJA(4)
975 sage: J.rank.clear_cache()
981 sage: J = JordanSpinEJA(4)
982 sage: J.rank.clear_cache()
988 sage: J = RealSymmetricEJA(3)
989 sage: J.rank.clear_cache()
995 sage: J = ComplexHermitianEJA(2)
996 sage: J.rank.clear_cache()
1002 sage: J = QuaternionHermitianEJA(2)
1003 sage: J.rank.clear_cache()
1007 return len(self
._charpoly
_coefficients
())
1010 def vector_space(self
):
1012 Return the vector space that underlies this algebra.
1016 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1020 sage: J = RealSymmetricEJA(2)
1021 sage: J.vector_space()
1022 Vector space of dimension 3 over...
1025 return self
.zero().to_vector().parent().ambient_vector_space()
1028 Element
= FiniteDimensionalEuclideanJordanAlgebraElement
1031 class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra
):
1033 Algebras whose basis consists of vectors with rational
1034 entries. Equivalently, algebras whose multiplication tables
1035 contain only rational coefficients.
1037 When an EJA has a basis that can be made rational, we can speed up
1038 the computation of its characteristic polynomial by doing it over
1039 ``QQ``. All of the named EJA constructors that we provide fall
1043 def _charpoly_coefficients(self
):
1045 Override the parent method with something that tries to compute
1046 over a faster (non-extension) field.
1050 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1054 The base ring of the resulting polynomial coefficients is what
1055 it should be, and not the rationals (unless the algebra was
1056 already over the rationals)::
1058 sage: J = JordanSpinEJA(3)
1059 sage: J._charpoly_coefficients()
1060 (X1^2 - X2^2 - X3^2, -2*X1)
1061 sage: a0 = J._charpoly_coefficients()[0]
1063 Algebraic Real Field
1064 sage: a0.base_ring()
1065 Algebraic Real Field
1068 if self
.base_ring() is QQ
:
1069 # There's no need to construct *another* algebra over the
1070 # rationals if this one is already over the rationals.
1071 superclass
= super(RationalBasisEuclideanJordanAlgebra
, self
)
1072 return superclass
._charpoly
_coefficients
()
1075 map(lambda x
: x
.to_vector(), ls
)
1076 for ls
in self
._multiplication
_table
1079 # Do the computation over the rationals. The answer will be
1080 # the same, because our basis coordinates are (essentially)
1082 J
= FiniteDimensionalEuclideanJordanAlgebra(QQ
,
1086 a
= J
._charpoly
_coefficients
()
1087 return tuple(map(lambda x
: x
.change_ring(self
.base_ring()), a
))
1090 class ConcreteEuclideanJordanAlgebra
:
1092 A class for the Euclidean Jordan algebras that we know by name.
1094 These are the Jordan algebras whose basis, multiplication table,
1095 rank, and so on are known a priori. More to the point, they are
1096 the Euclidean Jordan algebras for which we are able to conjure up
1097 a "random instance."
1101 sage: from mjo.eja.eja_algebra import ConcreteEuclideanJordanAlgebra
1105 Our natural basis is normalized with respect to the natural inner
1106 product unless we specify otherwise::
1108 sage: set_random_seed()
1109 sage: J = ConcreteEuclideanJordanAlgebra.random_instance()
1110 sage: all( b.norm() == 1 for b in J.gens() )
1113 Since our natural basis is normalized with respect to the natural
1114 inner product, and since we know that this algebra is an EJA, any
1115 left-multiplication operator's matrix will be symmetric because
1116 natural->EJA basis representation is an isometry and within the EJA
1117 the operator is self-adjoint by the Jordan axiom::
1119 sage: set_random_seed()
1120 sage: J = ConcreteEuclideanJordanAlgebra.random_instance()
1121 sage: x = J.random_element()
1122 sage: x.operator().matrix().is_symmetric()
1128 def _max_random_instance_size():
1130 Return an integer "size" that is an upper bound on the size of
1131 this algebra when it is used in a random test
1132 case. Unfortunately, the term "size" is ambiguous -- when
1133 dealing with `R^n` under either the Hadamard or Jordan spin
1134 product, the "size" refers to the dimension `n`. When dealing
1135 with a matrix algebra (real symmetric or complex/quaternion
1136 Hermitian), it refers to the size of the matrix, which is far
1137 less than the dimension of the underlying vector space.
1139 This method must be implemented in each subclass.
1141 raise NotImplementedError
1144 def random_instance(cls
, field
=AA
, **kwargs
):
1146 Return a random instance of this type of algebra.
1148 This method should be implemented in each subclass.
1150 from sage
.misc
.prandom
import choice
1151 eja_class
= choice(cls
.__subclasses
__())
1152 return eja_class
.random_instance(field
)
1155 class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra
):
1157 def __init__(self
, field
, basis
, normalize_basis
=True, **kwargs
):
1159 Compared to the superclass constructor, we take a basis instead of
1160 a multiplication table because the latter can be computed in terms
1161 of the former when the product is known (like it is here).
1163 # Used in this class's fast _charpoly_coefficients() override.
1164 self
._basis
_normalizers
= None
1166 # We're going to loop through this a few times, so now's a good
1167 # time to ensure that it isn't a generator expression.
1168 basis
= tuple(basis
)
1170 algebra_dim
= len(basis
)
1171 degree
= 0 # size of the matrices
1173 degree
= basis
[0].nrows()
1175 if algebra_dim
> 1 and normalize_basis
:
1176 # We'll need sqrt(2) to normalize the basis, and this
1177 # winds up in the multiplication table, so the whole
1178 # algebra needs to be over the field extension.
1179 R
= PolynomialRing(field
, 'z')
1182 if p
.is_irreducible():
1183 field
= field
.extension(p
, 'sqrt2', embedding
=RLF(2).sqrt())
1184 basis
= tuple( s
.change_ring(field
) for s
in basis
)
1185 self
._basis
_normalizers
= tuple(
1186 ~
(self
.natural_inner_product(s
,s
).sqrt()) for s
in basis
)
1187 basis
= tuple(s
*c
for (s
,c
) in zip(basis
,self
._basis
_normalizers
))
1189 # Now compute the multiplication and inner product tables.
1190 # We have to do this *after* normalizing the basis, because
1191 # scaling affects the answers.
1192 V
= VectorSpace(field
, degree
**2)
1193 W
= V
.span_of_basis( _mat2vec(s
) for s
in basis
)
1194 mult_table
= [[W
.zero() for j
in range(algebra_dim
)]
1195 for i
in range(algebra_dim
)]
1196 ip_table
= [[W
.zero() for j
in range(algebra_dim
)]
1197 for i
in range(algebra_dim
)]
1198 for i
in range(algebra_dim
):
1199 for j
in range(algebra_dim
):
1200 mat_entry
= (basis
[i
]*basis
[j
] + basis
[j
]*basis
[i
])/2
1201 mult_table
[i
][j
] = W
.coordinate_vector(_mat2vec(mat_entry
))
1204 # HACK: ignore the error here if we don't need the
1205 # inner product (as is the case when we construct
1206 # a dummy QQ-algebra for fast charpoly coefficients.
1207 ip_table
[i
][j
] = self
.natural_inner_product(basis
[i
],
1214 self
._inner
_product
_matrix
= matrix(field
,ip_table
)
1218 super(MatrixEuclideanJordanAlgebra
, self
).__init
__(field
,
1220 natural_basis
=basis
,
1223 if algebra_dim
== 0:
1224 self
.one
.set_cache(self
.zero())
1226 n
= basis
[0].nrows()
1227 # The identity wrt (A,B) -> (AB + BA)/2 is independent of the
1228 # details of this algebra.
1229 self
.one
.set_cache(self(matrix
.identity(field
,n
)))
1233 def _charpoly_coefficients(self
):
1235 Override the parent method with something that tries to compute
1236 over a faster (non-extension) field.
1238 if self
._basis
_normalizers
is None or self
.base_ring() is QQ
:
1239 # We didn't normalize, or the basis we started with had
1240 # entries in a nice field already. Just compute the thing.
1241 return super(MatrixEuclideanJordanAlgebra
, self
)._charpoly
_coefficients
()
1243 basis
= ( (b
/n
) for (b
,n
) in zip(self
.natural_basis(),
1244 self
._basis
_normalizers
) )
1246 # Do this over the rationals and convert back at the end.
1247 # Only works because we know the entries of the basis are
1248 # integers. The argument ``check_axioms=False`` is required
1249 # because the trace inner-product method for this
1250 # class is a stub and can't actually be checked.
1251 J
= MatrixEuclideanJordanAlgebra(QQ
,
1253 normalize_basis
=False,
1256 a
= J
._charpoly
_coefficients
()
1258 # Unfortunately, changing the basis does change the
1259 # coefficients of the characteristic polynomial, but since
1260 # these are really the coefficients of the "characteristic
1261 # polynomial of" function, everything is still nice and
1262 # unevaluated. It's therefore "obvious" how scaling the
1263 # basis affects the coordinate variables X1, X2, et
1264 # cetera. Scaling the first basis vector up by "n" adds a
1265 # factor of 1/n into every "X1" term, for example. So here
1266 # we simply undo the basis_normalizer scaling that we
1267 # performed earlier.
1269 # The a[0] access here is safe because trivial algebras
1270 # won't have any basis normalizers and therefore won't
1271 # make it to this "else" branch.
1272 XS
= a
[0].parent().gens()
1273 subs_dict
= { XS
[i
]: self
._basis
_normalizers
[i
]*XS
[i
]
1274 for i
in range(len(XS
)) }
1275 return tuple( a_i
.subs(subs_dict
) for a_i
in a
)
1281 Embed the matrix ``M`` into a space of real matrices.
1283 The matrix ``M`` can have entries in any field at the moment:
1284 the real numbers, complex numbers, or quaternions. And although
1285 they are not a field, we can probably support octonions at some
1286 point, too. This function returns a real matrix that "acts like"
1287 the original with respect to matrix multiplication; i.e.
1289 real_embed(M*N) = real_embed(M)*real_embed(N)
1292 raise NotImplementedError
1296 def real_unembed(M
):
1298 The inverse of :meth:`real_embed`.
1300 raise NotImplementedError
1303 def natural_inner_product(cls
,X
,Y
):
1304 Xu
= cls
.real_unembed(X
)
1305 Yu
= cls
.real_unembed(Y
)
1306 tr
= (Xu
*Yu
).trace()
1309 # Works in QQ, AA, RDF, et cetera.
1311 except AttributeError:
1312 # A quaternion doesn't have a real() method, but does
1313 # have coefficient_tuple() method that returns the
1314 # coefficients of 1, i, j, and k -- in that order.
1315 return tr
.coefficient_tuple()[0]
1318 class RealMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1322 The identity function, for embedding real matrices into real
1328 def real_unembed(M
):
1330 The identity function, for unembedding real matrices from real
1336 class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra
,
1337 ConcreteEuclideanJordanAlgebra
):
1339 The rank-n simple EJA consisting of real symmetric n-by-n
1340 matrices, the usual symmetric Jordan product, and the trace inner
1341 product. It has dimension `(n^2 + n)/2` over the reals.
1345 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1349 sage: J = RealSymmetricEJA(2)
1350 sage: e0, e1, e2 = J.gens()
1358 In theory, our "field" can be any subfield of the reals::
1360 sage: RealSymmetricEJA(2, RDF)
1361 Euclidean Jordan algebra of dimension 3 over Real Double Field
1362 sage: RealSymmetricEJA(2, RR)
1363 Euclidean Jordan algebra of dimension 3 over Real Field with
1364 53 bits of precision
1368 The dimension of this algebra is `(n^2 + n) / 2`::
1370 sage: set_random_seed()
1371 sage: n_max = RealSymmetricEJA._max_random_instance_size()
1372 sage: n = ZZ.random_element(1, n_max)
1373 sage: J = RealSymmetricEJA(n)
1374 sage: J.dimension() == (n^2 + n)/2
1377 The Jordan multiplication is what we think it is::
1379 sage: set_random_seed()
1380 sage: J = RealSymmetricEJA.random_instance()
1381 sage: x,y = J.random_elements(2)
1382 sage: actual = (x*y).natural_representation()
1383 sage: X = x.natural_representation()
1384 sage: Y = y.natural_representation()
1385 sage: expected = (X*Y + Y*X)/2
1386 sage: actual == expected
1388 sage: J(expected) == x*y
1391 We can change the generator prefix::
1393 sage: RealSymmetricEJA(3, prefix='q').gens()
1394 (q0, q1, q2, q3, q4, q5)
1396 We can construct the (trivial) algebra of rank zero::
1398 sage: RealSymmetricEJA(0)
1399 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1403 def _denormalized_basis(cls
, n
, field
):
1405 Return a basis for the space of real symmetric n-by-n matrices.
1409 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1413 sage: set_random_seed()
1414 sage: n = ZZ.random_element(1,5)
1415 sage: B = RealSymmetricEJA._denormalized_basis(n,QQ)
1416 sage: all( M.is_symmetric() for M in B)
1420 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1424 for j
in range(i
+1):
1425 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
1429 Sij
= Eij
+ Eij
.transpose()
1435 def _max_random_instance_size():
1436 return 4 # Dimension 10
1439 def random_instance(cls
, field
=AA
, **kwargs
):
1441 Return a random instance of this type of algebra.
1443 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1444 return cls(n
, field
, **kwargs
)
1446 def __init__(self
, n
, field
=AA
, **kwargs
):
1447 basis
= self
._denormalized
_basis
(n
, field
)
1448 super(RealSymmetricEJA
, self
).__init
__(field
,
1452 self
.rank
.set_cache(n
)
1455 class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1459 Embed the n-by-n complex matrix ``M`` into the space of real
1460 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1461 bi` to the block matrix ``[[a,b],[-b,a]]``.
1465 sage: from mjo.eja.eja_algebra import \
1466 ....: ComplexMatrixEuclideanJordanAlgebra
1470 sage: F = QuadraticField(-1, 'I')
1471 sage: x1 = F(4 - 2*i)
1472 sage: x2 = F(1 + 2*i)
1475 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1476 sage: ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1485 Embedding is a homomorphism (isomorphism, in fact)::
1487 sage: set_random_seed()
1488 sage: n = ZZ.random_element(3)
1489 sage: F = QuadraticField(-1, 'I')
1490 sage: X = random_matrix(F, n)
1491 sage: Y = random_matrix(F, n)
1492 sage: Xe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X)
1493 sage: Ye = ComplexMatrixEuclideanJordanAlgebra.real_embed(Y)
1494 sage: XYe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1501 raise ValueError("the matrix 'M' must be square")
1503 # We don't need any adjoined elements...
1504 field
= M
.base_ring().base_ring()
1508 a
= z
.list()[0] # real part, I guess
1509 b
= z
.list()[1] # imag part, I guess
1510 blocks
.append(matrix(field
, 2, [[a
,b
],[-b
,a
]]))
1512 return matrix
.block(field
, n
, blocks
)
1516 def real_unembed(M
):
1518 The inverse of _embed_complex_matrix().
1522 sage: from mjo.eja.eja_algebra import \
1523 ....: ComplexMatrixEuclideanJordanAlgebra
1527 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1528 ....: [-2, 1, -4, 3],
1529 ....: [ 9, 10, 11, 12],
1530 ....: [-10, 9, -12, 11] ])
1531 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(A)
1533 [ 10*I + 9 12*I + 11]
1537 Unembedding is the inverse of embedding::
1539 sage: set_random_seed()
1540 sage: F = QuadraticField(-1, 'I')
1541 sage: M = random_matrix(F, 3)
1542 sage: Me = ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1543 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1549 raise ValueError("the matrix 'M' must be square")
1550 if not n
.mod(2).is_zero():
1551 raise ValueError("the matrix 'M' must be a complex embedding")
1553 # If "M" was normalized, its base ring might have roots
1554 # adjoined and they can stick around after unembedding.
1555 field
= M
.base_ring()
1556 R
= PolynomialRing(field
, 'z')
1559 # Sage doesn't know how to embed AA into QQbar, i.e. how
1560 # to adjoin sqrt(-1) to AA.
1563 F
= field
.extension(z
**2 + 1, 'I', embedding
=CLF(-1).sqrt())
1566 # Go top-left to bottom-right (reading order), converting every
1567 # 2-by-2 block we see to a single complex element.
1569 for k
in range(n
/2):
1570 for j
in range(n
/2):
1571 submat
= M
[2*k
:2*k
+2,2*j
:2*j
+2]
1572 if submat
[0,0] != submat
[1,1]:
1573 raise ValueError('bad on-diagonal submatrix')
1574 if submat
[0,1] != -submat
[1,0]:
1575 raise ValueError('bad off-diagonal submatrix')
1576 z
= submat
[0,0] + submat
[0,1]*i
1579 return matrix(F
, n
/2, elements
)
1583 def natural_inner_product(cls
,X
,Y
):
1585 Compute a natural inner product in this algebra directly from
1590 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1594 This gives the same answer as the slow, default method implemented
1595 in :class:`MatrixEuclideanJordanAlgebra`::
1597 sage: set_random_seed()
1598 sage: J = ComplexHermitianEJA.random_instance()
1599 sage: x,y = J.random_elements(2)
1600 sage: Xe = x.natural_representation()
1601 sage: Ye = y.natural_representation()
1602 sage: X = ComplexHermitianEJA.real_unembed(Xe)
1603 sage: Y = ComplexHermitianEJA.real_unembed(Ye)
1604 sage: expected = (X*Y).trace().real()
1605 sage: actual = ComplexHermitianEJA.natural_inner_product(Xe,Ye)
1606 sage: actual == expected
1610 return RealMatrixEuclideanJordanAlgebra
.natural_inner_product(X
,Y
)/2
1613 class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra
,
1614 ConcreteEuclideanJordanAlgebra
):
1616 The rank-n simple EJA consisting of complex Hermitian n-by-n
1617 matrices over the real numbers, the usual symmetric Jordan product,
1618 and the real-part-of-trace inner product. It has dimension `n^2` over
1623 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1627 In theory, our "field" can be any subfield of the reals::
1629 sage: ComplexHermitianEJA(2, RDF)
1630 Euclidean Jordan algebra of dimension 4 over Real Double Field
1631 sage: ComplexHermitianEJA(2, RR)
1632 Euclidean Jordan algebra of dimension 4 over Real Field with
1633 53 bits of precision
1637 The dimension of this algebra is `n^2`::
1639 sage: set_random_seed()
1640 sage: n_max = ComplexHermitianEJA._max_random_instance_size()
1641 sage: n = ZZ.random_element(1, n_max)
1642 sage: J = ComplexHermitianEJA(n)
1643 sage: J.dimension() == n^2
1646 The Jordan multiplication is what we think it is::
1648 sage: set_random_seed()
1649 sage: J = ComplexHermitianEJA.random_instance()
1650 sage: x,y = J.random_elements(2)
1651 sage: actual = (x*y).natural_representation()
1652 sage: X = x.natural_representation()
1653 sage: Y = y.natural_representation()
1654 sage: expected = (X*Y + Y*X)/2
1655 sage: actual == expected
1657 sage: J(expected) == x*y
1660 We can change the generator prefix::
1662 sage: ComplexHermitianEJA(2, prefix='z').gens()
1665 We can construct the (trivial) algebra of rank zero::
1667 sage: ComplexHermitianEJA(0)
1668 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1673 def _denormalized_basis(cls
, n
, field
):
1675 Returns a basis for the space of complex Hermitian n-by-n matrices.
1677 Why do we embed these? Basically, because all of numerical linear
1678 algebra assumes that you're working with vectors consisting of `n`
1679 entries from a field and scalars from the same field. There's no way
1680 to tell SageMath that (for example) the vectors contain complex
1681 numbers, while the scalar field is real.
1685 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1689 sage: set_random_seed()
1690 sage: n = ZZ.random_element(1,5)
1691 sage: field = QuadraticField(2, 'sqrt2')
1692 sage: B = ComplexHermitianEJA._denormalized_basis(n, field)
1693 sage: all( M.is_symmetric() for M in B)
1697 R
= PolynomialRing(field
, 'z')
1699 F
= field
.extension(z
**2 + 1, 'I')
1702 # This is like the symmetric case, but we need to be careful:
1704 # * We want conjugate-symmetry, not just symmetry.
1705 # * The diagonal will (as a result) be real.
1709 for j
in range(i
+1):
1710 Eij
= matrix(F
, n
, lambda k
,l
: k
==i
and l
==j
)
1712 Sij
= cls
.real_embed(Eij
)
1715 # The second one has a minus because it's conjugated.
1716 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
1718 Sij_imag
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
1721 # Since we embedded these, we can drop back to the "field" that we
1722 # started with instead of the complex extension "F".
1723 return ( s
.change_ring(field
) for s
in S
)
1726 def __init__(self
, n
, field
=AA
, **kwargs
):
1727 basis
= self
._denormalized
_basis
(n
,field
)
1728 super(ComplexHermitianEJA
,self
).__init
__(field
,
1732 self
.rank
.set_cache(n
)
1735 def _max_random_instance_size():
1736 return 3 # Dimension 9
1739 def random_instance(cls
, field
=AA
, **kwargs
):
1741 Return a random instance of this type of algebra.
1743 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1744 return cls(n
, field
, **kwargs
)
1746 class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1750 Embed the n-by-n quaternion matrix ``M`` into the space of real
1751 matrices of size 4n-by-4n by first sending each quaternion entry `z
1752 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
1753 c+di],[-c + di, a-bi]]`, and then embedding those into a real
1758 sage: from mjo.eja.eja_algebra import \
1759 ....: QuaternionMatrixEuclideanJordanAlgebra
1763 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1764 sage: i,j,k = Q.gens()
1765 sage: x = 1 + 2*i + 3*j + 4*k
1766 sage: M = matrix(Q, 1, [[x]])
1767 sage: QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1773 Embedding is a homomorphism (isomorphism, in fact)::
1775 sage: set_random_seed()
1776 sage: n = ZZ.random_element(2)
1777 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1778 sage: X = random_matrix(Q, n)
1779 sage: Y = random_matrix(Q, n)
1780 sage: Xe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X)
1781 sage: Ye = QuaternionMatrixEuclideanJordanAlgebra.real_embed(Y)
1782 sage: XYe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1787 quaternions
= M
.base_ring()
1790 raise ValueError("the matrix 'M' must be square")
1792 F
= QuadraticField(-1, 'I')
1797 t
= z
.coefficient_tuple()
1802 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
1803 [-c
+ d
*i
, a
- b
*i
]])
1804 realM
= ComplexMatrixEuclideanJordanAlgebra
.real_embed(cplxM
)
1805 blocks
.append(realM
)
1807 # We should have real entries by now, so use the realest field
1808 # we've got for the return value.
1809 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
1814 def real_unembed(M
):
1816 The inverse of _embed_quaternion_matrix().
1820 sage: from mjo.eja.eja_algebra import \
1821 ....: QuaternionMatrixEuclideanJordanAlgebra
1825 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
1826 ....: [-2, 1, -4, 3],
1827 ....: [-3, 4, 1, -2],
1828 ....: [-4, -3, 2, 1]])
1829 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(M)
1830 [1 + 2*i + 3*j + 4*k]
1834 Unembedding is the inverse of embedding::
1836 sage: set_random_seed()
1837 sage: Q = QuaternionAlgebra(QQ, -1, -1)
1838 sage: M = random_matrix(Q, 3)
1839 sage: Me = QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1840 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1846 raise ValueError("the matrix 'M' must be square")
1847 if not n
.mod(4).is_zero():
1848 raise ValueError("the matrix 'M' must be a quaternion embedding")
1850 # Use the base ring of the matrix to ensure that its entries can be
1851 # multiplied by elements of the quaternion algebra.
1852 field
= M
.base_ring()
1853 Q
= QuaternionAlgebra(field
,-1,-1)
1856 # Go top-left to bottom-right (reading order), converting every
1857 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
1860 for l
in range(n
/4):
1861 for m
in range(n
/4):
1862 submat
= ComplexMatrixEuclideanJordanAlgebra
.real_unembed(
1863 M
[4*l
:4*l
+4,4*m
:4*m
+4] )
1864 if submat
[0,0] != submat
[1,1].conjugate():
1865 raise ValueError('bad on-diagonal submatrix')
1866 if submat
[0,1] != -submat
[1,0].conjugate():
1867 raise ValueError('bad off-diagonal submatrix')
1868 z
= submat
[0,0].real()
1869 z
+= submat
[0,0].imag()*i
1870 z
+= submat
[0,1].real()*j
1871 z
+= submat
[0,1].imag()*k
1874 return matrix(Q
, n
/4, elements
)
1878 def natural_inner_product(cls
,X
,Y
):
1880 Compute a natural inner product in this algebra directly from
1885 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1889 This gives the same answer as the slow, default method implemented
1890 in :class:`MatrixEuclideanJordanAlgebra`::
1892 sage: set_random_seed()
1893 sage: J = QuaternionHermitianEJA.random_instance()
1894 sage: x,y = J.random_elements(2)
1895 sage: Xe = x.natural_representation()
1896 sage: Ye = y.natural_representation()
1897 sage: X = QuaternionHermitianEJA.real_unembed(Xe)
1898 sage: Y = QuaternionHermitianEJA.real_unembed(Ye)
1899 sage: expected = (X*Y).trace().coefficient_tuple()[0]
1900 sage: actual = QuaternionHermitianEJA.natural_inner_product(Xe,Ye)
1901 sage: actual == expected
1905 return RealMatrixEuclideanJordanAlgebra
.natural_inner_product(X
,Y
)/4
1908 class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra
,
1909 ConcreteEuclideanJordanAlgebra
):
1911 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
1912 matrices, the usual symmetric Jordan product, and the
1913 real-part-of-trace inner product. It has dimension `2n^2 - n` over
1918 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1922 In theory, our "field" can be any subfield of the reals::
1924 sage: QuaternionHermitianEJA(2, RDF)
1925 Euclidean Jordan algebra of dimension 6 over Real Double Field
1926 sage: QuaternionHermitianEJA(2, RR)
1927 Euclidean Jordan algebra of dimension 6 over Real Field with
1928 53 bits of precision
1932 The dimension of this algebra is `2*n^2 - n`::
1934 sage: set_random_seed()
1935 sage: n_max = QuaternionHermitianEJA._max_random_instance_size()
1936 sage: n = ZZ.random_element(1, n_max)
1937 sage: J = QuaternionHermitianEJA(n)
1938 sage: J.dimension() == 2*(n^2) - n
1941 The Jordan multiplication is what we think it is::
1943 sage: set_random_seed()
1944 sage: J = QuaternionHermitianEJA.random_instance()
1945 sage: x,y = J.random_elements(2)
1946 sage: actual = (x*y).natural_representation()
1947 sage: X = x.natural_representation()
1948 sage: Y = y.natural_representation()
1949 sage: expected = (X*Y + Y*X)/2
1950 sage: actual == expected
1952 sage: J(expected) == x*y
1955 We can change the generator prefix::
1957 sage: QuaternionHermitianEJA(2, prefix='a').gens()
1958 (a0, a1, a2, a3, a4, a5)
1960 We can construct the (trivial) algebra of rank zero::
1962 sage: QuaternionHermitianEJA(0)
1963 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1967 def _denormalized_basis(cls
, n
, field
):
1969 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
1971 Why do we embed these? Basically, because all of numerical
1972 linear algebra assumes that you're working with vectors consisting
1973 of `n` entries from a field and scalars from the same field. There's
1974 no way to tell SageMath that (for example) the vectors contain
1975 complex numbers, while the scalar field is real.
1979 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1983 sage: set_random_seed()
1984 sage: n = ZZ.random_element(1,5)
1985 sage: B = QuaternionHermitianEJA._denormalized_basis(n,QQ)
1986 sage: all( M.is_symmetric() for M in B )
1990 Q
= QuaternionAlgebra(QQ
,-1,-1)
1993 # This is like the symmetric case, but we need to be careful:
1995 # * We want conjugate-symmetry, not just symmetry.
1996 # * The diagonal will (as a result) be real.
2000 for j
in range(i
+1):
2001 Eij
= matrix(Q
, n
, lambda k
,l
: k
==i
and l
==j
)
2003 Sij
= cls
.real_embed(Eij
)
2006 # The second, third, and fourth ones have a minus
2007 # because they're conjugated.
2008 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
2010 Sij_I
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
2012 Sij_J
= cls
.real_embed(J
*Eij
- J
*Eij
.transpose())
2014 Sij_K
= cls
.real_embed(K
*Eij
- K
*Eij
.transpose())
2017 # Since we embedded these, we can drop back to the "field" that we
2018 # started with instead of the quaternion algebra "Q".
2019 return ( s
.change_ring(field
) for s
in S
)
2022 def __init__(self
, n
, field
=AA
, **kwargs
):
2023 basis
= self
._denormalized
_basis
(n
,field
)
2024 super(QuaternionHermitianEJA
,self
).__init
__(field
,
2028 self
.rank
.set_cache(n
)
2031 def _max_random_instance_size():
2033 The maximum rank of a random QuaternionHermitianEJA.
2035 return 2 # Dimension 6
2038 def random_instance(cls
, field
=AA
, **kwargs
):
2040 Return a random instance of this type of algebra.
2042 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2043 return cls(n
, field
, **kwargs
)
2046 class HadamardEJA(RationalBasisEuclideanJordanAlgebra
,
2047 ConcreteEuclideanJordanAlgebra
):
2049 Return the Euclidean Jordan Algebra corresponding to the set
2050 `R^n` under the Hadamard product.
2052 Note: this is nothing more than the Cartesian product of ``n``
2053 copies of the spin algebra. Once Cartesian product algebras
2054 are implemented, this can go.
2058 sage: from mjo.eja.eja_algebra import HadamardEJA
2062 This multiplication table can be verified by hand::
2064 sage: J = HadamardEJA(3)
2065 sage: e0,e1,e2 = J.gens()
2081 We can change the generator prefix::
2083 sage: HadamardEJA(3, prefix='r').gens()
2087 def __init__(self
, n
, field
=AA
, **kwargs
):
2088 V
= VectorSpace(field
, n
)
2089 mult_table
= [ [ V
.gen(i
)*(i
== j
) for j
in range(n
) ]
2092 # Inner products are real numbers and not algebra
2093 # elements, so once we turn the algebra element
2094 # into a vector in inner_product(), we never go
2095 # back. As a result -- contrary to what we do with
2096 # self._multiplication_table -- we store the inner
2097 # product table as a plain old matrix and not as
2098 # an algebra operator.
2099 ip_table
= matrix
.identity(field
,n
)
2100 self
._inner
_product
_matrix
= ip_table
2102 super(HadamardEJA
, self
).__init
__(field
,
2106 self
.rank
.set_cache(n
)
2109 self
.one
.set_cache( self
.zero() )
2111 self
.one
.set_cache( sum(self
.gens()) )
2114 def _max_random_instance_size():
2116 The maximum dimension of a random HadamardEJA.
2121 def random_instance(cls
, field
=AA
, **kwargs
):
2123 Return a random instance of this type of algebra.
2125 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2126 return cls(n
, field
, **kwargs
)
2129 class BilinearFormEJA(RationalBasisEuclideanJordanAlgebra
,
2130 ConcreteEuclideanJordanAlgebra
):
2132 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2133 with the half-trace inner product and jordan product ``x*y =
2134 (<Bx,y>,y_bar>, x0*y_bar + y0*x_bar)`` where `B = 1 \times B22` is
2135 a symmetric positive-definite "bilinear form" matrix. Its
2136 dimension is the size of `B`, and it has rank two in dimensions
2137 larger than two. It reduces to the ``JordanSpinEJA`` when `B` is
2138 the identity matrix of order ``n``.
2140 We insist that the one-by-one upper-left identity block of `B` be
2141 passed in as well so that we can be passed a matrix of size zero
2142 to construct a trivial algebra.
2146 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
2147 ....: JordanSpinEJA)
2151 When no bilinear form is specified, the identity matrix is used,
2152 and the resulting algebra is the Jordan spin algebra::
2154 sage: B = matrix.identity(AA,3)
2155 sage: J0 = BilinearFormEJA(B)
2156 sage: J1 = JordanSpinEJA(3)
2157 sage: J0.multiplication_table() == J0.multiplication_table()
2160 An error is raised if the matrix `B` does not correspond to a
2161 positive-definite bilinear form::
2163 sage: B = matrix.random(QQ,2,3)
2164 sage: J = BilinearFormEJA(B)
2165 Traceback (most recent call last):
2167 ValueError: bilinear form is not positive-definite
2168 sage: B = matrix.zero(QQ,3)
2169 sage: J = BilinearFormEJA(B)
2170 Traceback (most recent call last):
2172 ValueError: bilinear form is not positive-definite
2176 We can create a zero-dimensional algebra::
2178 sage: B = matrix.identity(AA,0)
2179 sage: J = BilinearFormEJA(B)
2183 We can check the multiplication condition given in the Jordan, von
2184 Neumann, and Wigner paper (and also discussed on my "On the
2185 symmetry..." paper). Note that this relies heavily on the standard
2186 choice of basis, as does anything utilizing the bilinear form matrix::
2188 sage: set_random_seed()
2189 sage: n = ZZ.random_element(5)
2190 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2191 sage: B11 = matrix.identity(QQ,1)
2192 sage: B22 = M.transpose()*M
2193 sage: B = block_matrix(2,2,[ [B11,0 ],
2195 sage: J = BilinearFormEJA(B)
2196 sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
2197 sage: V = J.vector_space()
2198 sage: sis = [ J.from_vector(V([0] + (M.inverse()*ei).list()))
2199 ....: for ei in eis ]
2200 sage: actual = [ sis[i]*sis[j]
2201 ....: for i in range(n-1)
2202 ....: for j in range(n-1) ]
2203 sage: expected = [ J.one() if i == j else J.zero()
2204 ....: for i in range(n-1)
2205 ....: for j in range(n-1) ]
2206 sage: actual == expected
2209 def __init__(self
, B
, field
=AA
, **kwargs
):
2212 if not B
.is_positive_definite():
2213 raise ValueError("bilinear form is not positive-definite")
2215 V
= VectorSpace(field
, n
)
2216 mult_table
= [[V
.zero() for j
in range(n
)] for i
in range(n
)]
2225 z0
= (B
*x
).inner_product(y
)
2226 zbar
= y0
*xbar
+ x0
*ybar
2227 z
= V([z0
] + zbar
.list())
2228 mult_table
[i
][j
] = z
2230 # Inner products are real numbers and not algebra
2231 # elements, so once we turn the algebra element
2232 # into a vector in inner_product(), we never go
2233 # back. As a result -- contrary to what we do with
2234 # self._multiplication_table -- we store the inner
2235 # product table as a plain old matrix and not as
2236 # an algebra operator.
2238 self
._inner
_product
_matrix
= ip_table
2240 super(BilinearFormEJA
, self
).__init
__(field
,
2245 # The rank of this algebra is two, unless we're in a
2246 # one-dimensional ambient space (because the rank is bounded
2247 # by the ambient dimension).
2248 self
.rank
.set_cache(min(n
,2))
2251 self
.one
.set_cache( self
.zero() )
2253 self
.one
.set_cache( self
.monomial(0) )
2256 def _max_random_instance_size():
2258 The maximum dimension of a random BilinearFormEJA.
2263 def random_instance(cls
, field
=AA
, **kwargs
):
2265 Return a random instance of this algebra.
2267 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2269 B
= matrix
.identity(field
, n
)
2270 return cls(B
, field
, **kwargs
)
2272 B11
= matrix
.identity(field
,1)
2273 M
= matrix
.random(field
, n
-1)
2274 I
= matrix
.identity(field
, n
-1)
2275 alpha
= field
.zero()
2276 while alpha
.is_zero():
2277 alpha
= field
.random_element().abs()
2278 B22
= M
.transpose()*M
+ alpha
*I
2280 from sage
.matrix
.special
import block_matrix
2281 B
= block_matrix(2,2, [ [B11
, ZZ(0) ],
2284 return cls(B
, field
, **kwargs
)
2287 class JordanSpinEJA(BilinearFormEJA
):
2289 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2290 with the usual inner product and jordan product ``x*y =
2291 (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2296 sage: from mjo.eja.eja_algebra import JordanSpinEJA
2300 This multiplication table can be verified by hand::
2302 sage: J = JordanSpinEJA(4)
2303 sage: e0,e1,e2,e3 = J.gens()
2319 We can change the generator prefix::
2321 sage: JordanSpinEJA(2, prefix='B').gens()
2326 Ensure that we have the usual inner product on `R^n`::
2328 sage: set_random_seed()
2329 sage: J = JordanSpinEJA.random_instance()
2330 sage: x,y = J.random_elements(2)
2331 sage: actual = x.inner_product(y)
2332 sage: expected = x.to_vector().inner_product(y.to_vector())
2333 sage: actual == expected
2337 def __init__(self
, n
, field
=AA
, **kwargs
):
2338 # This is a special case of the BilinearFormEJA with the identity
2339 # matrix as its bilinear form.
2340 B
= matrix
.identity(field
, n
)
2341 super(JordanSpinEJA
, self
).__init
__(B
, field
, **kwargs
)
2344 def _max_random_instance_size():
2346 The maximum dimension of a random JordanSpinEJA.
2351 def random_instance(cls
, field
=AA
, **kwargs
):
2353 Return a random instance of this type of algebra.
2355 Needed here to override the implementation for ``BilinearFormEJA``.
2357 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2358 return cls(n
, field
, **kwargs
)
2361 class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra
,
2362 ConcreteEuclideanJordanAlgebra
):
2364 The trivial Euclidean Jordan algebra consisting of only a zero element.
2368 sage: from mjo.eja.eja_algebra import TrivialEJA
2372 sage: J = TrivialEJA()
2379 sage: 7*J.one()*12*J.one()
2381 sage: J.one().inner_product(J.one())
2383 sage: J.one().norm()
2385 sage: J.one().subalgebra_generated_by()
2386 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2391 def __init__(self
, field
=AA
, **kwargs
):
2393 self
._inner
_product
_matrix
= matrix(field
,0)
2394 super(TrivialEJA
, self
).__init
__(field
,
2398 # The rank is zero using my definition, namely the dimension of the
2399 # largest subalgebra generated by any element.
2400 self
.rank
.set_cache(0)
2401 self
.one
.set_cache( self
.zero() )
2404 def random_instance(cls
, field
=AA
, **kwargs
):
2405 # We don't take a "size" argument so the superclass method is
2406 # inappropriate for us.
2407 return cls(field
, **kwargs
)
2409 class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra
):
2411 The external (orthogonal) direct sum of two other Euclidean Jordan
2412 algebras. Essentially the Cartesian product of its two factors.
2413 Every Euclidean Jordan algebra decomposes into an orthogonal
2414 direct sum of simple Euclidean Jordan algebras, so no generality
2415 is lost by providing only this construction.
2419 sage: from mjo.eja.eja_algebra import (random_eja,
2421 ....: RealSymmetricEJA,
2426 sage: J1 = HadamardEJA(2)
2427 sage: J2 = RealSymmetricEJA(3)
2428 sage: J = DirectSumEJA(J1,J2)
2436 The external direct sum construction is only valid when the two factors
2437 have the same base ring; an error is raised otherwise::
2439 sage: set_random_seed()
2440 sage: J1 = random_eja(AA)
2441 sage: J2 = random_eja(QQ)
2442 sage: J = DirectSumEJA(J1,J2)
2443 Traceback (most recent call last):
2445 ValueError: algebras must share the same base field
2448 def __init__(self
, J1
, J2
, **kwargs
):
2449 if J1
.base_ring() != J2
.base_ring():
2450 raise ValueError("algebras must share the same base field")
2451 field
= J1
.base_ring()
2453 self
._factors
= (J1
, J2
)
2457 V
= VectorSpace(field
, n
)
2458 mult_table
= [ [ V
.zero() for j
in range(n
) ]
2462 p
= (J1
.monomial(i
)*J1
.monomial(j
)).to_vector()
2463 mult_table
[i
][j
] = V(p
.list() + [field
.zero()]*n2
)
2467 p
= (J2
.monomial(i
)*J2
.monomial(j
)).to_vector()
2468 mult_table
[n1
+i
][n1
+j
] = V([field
.zero()]*n1
+ p
.list())
2470 super(DirectSumEJA
, self
).__init
__(field
,
2474 self
.rank
.set_cache(J1
.rank() + J2
.rank())
2479 Return the pair of this algebra's factors.
2483 sage: from mjo.eja.eja_algebra import (HadamardEJA,
2484 ....: JordanSpinEJA,
2489 sage: J1 = HadamardEJA(2,QQ)
2490 sage: J2 = JordanSpinEJA(3,QQ)
2491 sage: J = DirectSumEJA(J1,J2)
2493 (Euclidean Jordan algebra of dimension 2 over Rational Field,
2494 Euclidean Jordan algebra of dimension 3 over Rational Field)
2497 return self
._factors
2499 def projections(self
):
2501 Return a pair of projections onto this algebra's factors.
2505 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
2506 ....: ComplexHermitianEJA,
2511 sage: J1 = JordanSpinEJA(2)
2512 sage: J2 = ComplexHermitianEJA(2)
2513 sage: J = DirectSumEJA(J1,J2)
2514 sage: (pi_left, pi_right) = J.projections()
2515 sage: J.one().to_vector()
2517 sage: pi_left(J.one()).to_vector()
2519 sage: pi_right(J.one()).to_vector()
2523 (J1
,J2
) = self
.factors()
2526 V_basis
= self
.vector_space().basis()
2527 # Need to specify the dimensions explicitly so that we don't
2528 # wind up with a zero-by-zero matrix when we want e.g. a
2529 # zero-by-two matrix (important for composing things).
2530 P1
= matrix(self
.base_ring(), m
, m
+n
, V_basis
[:m
])
2531 P2
= matrix(self
.base_ring(), n
, m
+n
, V_basis
[m
:])
2532 pi_left
= FiniteDimensionalEuclideanJordanAlgebraOperator(self
,J1
,P1
)
2533 pi_right
= FiniteDimensionalEuclideanJordanAlgebraOperator(self
,J2
,P2
)
2534 return (pi_left
, pi_right
)
2536 def inclusions(self
):
2538 Return the pair of inclusion maps from our factors into us.
2542 sage: from mjo.eja.eja_algebra import (random_eja,
2543 ....: JordanSpinEJA,
2544 ....: RealSymmetricEJA,
2549 sage: J1 = JordanSpinEJA(3)
2550 sage: J2 = RealSymmetricEJA(2)
2551 sage: J = DirectSumEJA(J1,J2)
2552 sage: (iota_left, iota_right) = J.inclusions()
2553 sage: iota_left(J1.zero()) == J.zero()
2555 sage: iota_right(J2.zero()) == J.zero()
2557 sage: J1.one().to_vector()
2559 sage: iota_left(J1.one()).to_vector()
2561 sage: J2.one().to_vector()
2563 sage: iota_right(J2.one()).to_vector()
2565 sage: J.one().to_vector()
2570 Composing a projection with the corresponding inclusion should
2571 produce the identity map, and mismatching them should produce
2574 sage: set_random_seed()
2575 sage: J1 = random_eja()
2576 sage: J2 = random_eja()
2577 sage: J = DirectSumEJA(J1,J2)
2578 sage: (iota_left, iota_right) = J.inclusions()
2579 sage: (pi_left, pi_right) = J.projections()
2580 sage: pi_left*iota_left == J1.one().operator()
2582 sage: pi_right*iota_right == J2.one().operator()
2584 sage: (pi_left*iota_right).is_zero()
2586 sage: (pi_right*iota_left).is_zero()
2590 (J1
,J2
) = self
.factors()
2593 V_basis
= self
.vector_space().basis()
2594 # Need to specify the dimensions explicitly so that we don't
2595 # wind up with a zero-by-zero matrix when we want e.g. a
2596 # two-by-zero matrix (important for composing things).
2597 I1
= matrix
.column(self
.base_ring(), m
, m
+n
, V_basis
[:m
])
2598 I2
= matrix
.column(self
.base_ring(), n
, m
+n
, V_basis
[m
:])
2599 iota_left
= FiniteDimensionalEuclideanJordanAlgebraOperator(J1
,self
,I1
)
2600 iota_right
= FiniteDimensionalEuclideanJordanAlgebraOperator(J2
,self
,I2
)
2601 return (iota_left
, iota_right
)
2603 def inner_product(self
, x
, y
):
2605 The standard Cartesian inner-product.
2607 We project ``x`` and ``y`` onto our factors, and add up the
2608 inner-products from the subalgebras.
2613 sage: from mjo.eja.eja_algebra import (HadamardEJA,
2614 ....: QuaternionHermitianEJA,
2619 sage: J1 = HadamardEJA(3,QQ)
2620 sage: J2 = QuaternionHermitianEJA(2,QQ,normalize_basis=False)
2621 sage: J = DirectSumEJA(J1,J2)
2626 sage: x1.inner_product(x2)
2628 sage: y1.inner_product(y2)
2630 sage: J.one().inner_product(J.one())
2634 (pi_left
, pi_right
) = self
.projections()
2640 return (x1
.inner_product(y1
) + x2
.inner_product(y2
))
2644 random_eja
= ConcreteEuclideanJordanAlgebra
.random_instance