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
._natural
_basis
= natural_basis
121 category
= MagmaticAlgebras(field
).FiniteDimensional()
122 category
= category
.WithBasis().Unital()
124 fda
= super(FiniteDimensionalEuclideanJordanAlgebra
, self
)
129 self
.print_options(bracket
='')
131 # The multiplication table we're given is necessarily in terms
132 # of vectors, because we don't have an algebra yet for
133 # anything to be an element of. However, it's faster in the
134 # long run to have the multiplication table be in terms of
135 # algebra elements. We do this after calling the superclass
136 # constructor so that from_vector() knows what to do.
137 self
._multiplication
_table
= [
138 list(map(lambda x
: self
.from_vector(x
), ls
))
143 if not self
._is
_commutative
():
144 raise ValueError("algebra is not commutative")
145 if not self
._is
_jordanian
():
146 raise ValueError("Jordan identity does not hold")
147 if not self
._inner
_product
_is
_associative
():
148 raise ValueError("inner product is not associative")
150 def _element_constructor_(self
, elt
):
152 Construct an element of this algebra from its natural
155 This gets called only after the parent element _call_ method
156 fails to find a coercion for the argument.
160 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
162 ....: RealSymmetricEJA)
166 The identity in `S^n` is converted to the identity in the EJA::
168 sage: J = RealSymmetricEJA(3)
169 sage: I = matrix.identity(QQ,3)
170 sage: J(I) == J.one()
173 This skew-symmetric matrix can't be represented in the EJA::
175 sage: J = RealSymmetricEJA(3)
176 sage: A = matrix(QQ,3, lambda i,j: i-j)
178 Traceback (most recent call last):
180 ValueError: not an element of this algebra
184 Ensure that we can convert any element of the two non-matrix
185 simple algebras (whose natural representations are their usual
186 vector representations) back and forth faithfully::
188 sage: set_random_seed()
189 sage: J = HadamardEJA.random_instance()
190 sage: x = J.random_element()
191 sage: J(x.to_vector().column()) == x
193 sage: J = JordanSpinEJA.random_instance()
194 sage: x = J.random_element()
195 sage: J(x.to_vector().column()) == x
199 msg
= "not an element of this algebra"
201 # The superclass implementation of random_element()
202 # needs to be able to coerce "0" into the algebra.
204 elif elt
in self
.base_ring():
205 # Ensure that no base ring -> algebra coercion is performed
206 # by this method. There's some stupidity in sage that would
207 # otherwise propagate to this method; for example, sage thinks
208 # that the integer 3 belongs to the space of 2-by-2 matrices.
209 raise ValueError(msg
)
211 natural_basis
= self
.natural_basis()
212 basis_space
= natural_basis
[0].matrix_space()
213 if elt
not in basis_space
:
214 raise ValueError(msg
)
216 # Thanks for nothing! Matrix spaces aren't vector spaces in
217 # Sage, so we have to figure out its natural-basis coordinates
218 # ourselves. We use the basis space's ring instead of the
219 # element's ring because the basis space might be an algebraic
220 # closure whereas the base ring of the 3-by-3 identity matrix
221 # could be QQ instead of QQbar.
222 V
= VectorSpace(basis_space
.base_ring(), elt
.nrows()*elt
.ncols())
223 W
= V
.span_of_basis( _mat2vec(s
) for s
in natural_basis
)
226 coords
= W
.coordinate_vector(_mat2vec(elt
))
227 except ArithmeticError: # vector is not in free module
228 raise ValueError(msg
)
230 return self
.from_vector(coords
)
234 Return a string representation of ``self``.
238 sage: from mjo.eja.eja_algebra import JordanSpinEJA
242 Ensure that it says what we think it says::
244 sage: JordanSpinEJA(2, field=AA)
245 Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
246 sage: JordanSpinEJA(3, field=RDF)
247 Euclidean Jordan algebra of dimension 3 over Real Double Field
250 fmt
= "Euclidean Jordan algebra of dimension {} over {}"
251 return fmt
.format(self
.dimension(), self
.base_ring())
253 def product_on_basis(self
, i
, j
):
254 return self
._multiplication
_table
[i
][j
]
256 def _is_commutative(self
):
258 Whether or not this algebra's multiplication table is commutative.
260 This method should of course always return ``True``, unless
261 this algebra was constructed with ``check_axioms=False`` and
262 passed an invalid multiplication table.
264 return all( self
.product_on_basis(i
,j
) == self
.product_on_basis(i
,j
)
265 for i
in range(self
.dimension())
266 for j
in range(self
.dimension()) )
268 def _is_jordanian(self
):
270 Whether or not this algebra's multiplication table respects the
271 Jordan identity `(x^{2})(xy) = x(x^{2}y)`.
273 We only check one arrangement of `x` and `y`, so for a
274 ``True`` result to be truly true, you should also check
275 :meth:`_is_commutative`. This method should of course always
276 return ``True``, unless this algebra was constructed with
277 ``check_axioms=False`` and passed an invalid multiplication table.
279 return all( (self
.monomial(i
)**2)*(self
.monomial(i
)*self
.monomial(j
))
281 (self
.monomial(i
))*((self
.monomial(i
)**2)*self
.monomial(j
))
282 for i
in range(self
.dimension())
283 for j
in range(self
.dimension()) )
285 def _inner_product_is_associative(self
):
287 Return whether or not this algebra's inner product `B` is
288 associative; that is, whether or not `B(xy,z) = B(x,yz)`.
290 This method should of course always return ``True``, unless
291 this algebra was constructed with ``check_axioms=False`` and
292 passed an invalid multiplication table.
295 # Used to check whether or not something is zero in an inexact
296 # ring. This number is sufficient to allow the construction of
297 # QuaternionHermitianEJA(2, RDF) with check_axioms=True.
300 for i
in range(self
.dimension()):
301 for j
in range(self
.dimension()):
302 for k
in range(self
.dimension()):
306 diff
= (x
*y
).inner_product(z
) - x
.inner_product(y
*z
)
308 if self
.base_ring().is_exact():
312 if diff
.abs() > epsilon
:
318 def characteristic_polynomial_of(self
):
320 Return the algebra's "characteristic polynomial of" function,
321 which is itself a multivariate polynomial that, when evaluated
322 at the coordinates of some algebra element, returns that
323 element's characteristic polynomial.
325 The resulting polynomial has `n+1` variables, where `n` is the
326 dimension of this algebra. The first `n` variables correspond to
327 the coordinates of an algebra element: when evaluated at the
328 coordinates of an algebra element with respect to a certain
329 basis, the result is a univariate polynomial (in the one
330 remaining variable ``t``), namely the characteristic polynomial
335 sage: from mjo.eja.eja_algebra import JordanSpinEJA, TrivialEJA
339 The characteristic polynomial in the spin algebra is given in
340 Alizadeh, Example 11.11::
342 sage: J = JordanSpinEJA(3)
343 sage: p = J.characteristic_polynomial_of(); p
344 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
345 sage: xvec = J.one().to_vector()
349 By definition, the characteristic polynomial is a monic
350 degree-zero polynomial in a rank-zero algebra. Note that
351 Cayley-Hamilton is indeed satisfied since the polynomial
352 ``1`` evaluates to the identity element of the algebra on
355 sage: J = TrivialEJA()
356 sage: J.characteristic_polynomial_of()
363 # The list of coefficient polynomials a_0, a_1, a_2, ..., a_(r-1).
364 a
= self
._charpoly
_coefficients
()
366 # We go to a bit of trouble here to reorder the
367 # indeterminates, so that it's easier to evaluate the
368 # characteristic polynomial at x's coordinates and get back
369 # something in terms of t, which is what we want.
370 S
= PolynomialRing(self
.base_ring(),'t')
374 S
= PolynomialRing(S
, R
.variable_names())
377 return (t
**r
+ sum( a
[k
]*(t
**k
) for k
in range(r
) ))
379 def coordinate_polynomial_ring(self
):
381 The multivariate polynomial ring in which this algebra's
382 :meth:`characteristic_polynomial_of` lives.
386 sage: from mjo.eja.eja_algebra import (HadamardEJA,
387 ....: RealSymmetricEJA)
391 sage: J = HadamardEJA(2)
392 sage: J.coordinate_polynomial_ring()
393 Multivariate Polynomial Ring in X1, X2...
394 sage: J = RealSymmetricEJA(3,QQ)
395 sage: J.coordinate_polynomial_ring()
396 Multivariate Polynomial Ring in X1, X2, X3, X4, X5, X6...
399 var_names
= tuple( "X%d" % z
for z
in range(1, self
.dimension()+1) )
400 return PolynomialRing(self
.base_ring(), var_names
)
402 def inner_product(self
, x
, y
):
404 The inner product associated with this Euclidean Jordan algebra.
406 Defaults to the trace inner product, but can be overridden by
407 subclasses if they are sure that the necessary properties are
412 sage: from mjo.eja.eja_algebra import (random_eja,
414 ....: BilinearFormEJA)
418 Our inner product is "associative," which means the following for
419 a symmetric bilinear form::
421 sage: set_random_seed()
422 sage: J = random_eja()
423 sage: x,y,z = J.random_elements(3)
424 sage: (x*y).inner_product(z) == y.inner_product(x*z)
429 Ensure that this is the usual inner product for the algebras
432 sage: set_random_seed()
433 sage: J = HadamardEJA.random_instance()
434 sage: x,y = J.random_elements(2)
435 sage: actual = x.inner_product(y)
436 sage: expected = x.to_vector().inner_product(y.to_vector())
437 sage: actual == expected
440 Ensure that this is one-half of the trace inner-product in a
441 BilinearFormEJA that isn't just the reals (when ``n`` isn't
442 one). This is in Faraut and Koranyi, and also my "On the
445 sage: set_random_seed()
446 sage: J = BilinearFormEJA.random_instance()
447 sage: n = J.dimension()
448 sage: x = J.random_element()
449 sage: y = J.random_element()
450 sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
453 B
= self
._inner
_product
_matrix
454 return (B
*x
.to_vector()).inner_product(y
.to_vector())
457 def is_trivial(self
):
459 Return whether or not this algebra is trivial.
461 A trivial algebra contains only the zero element.
465 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
470 sage: J = ComplexHermitianEJA(3)
476 sage: J = TrivialEJA()
481 return self
.dimension() == 0
484 def multiplication_table(self
):
486 Return a visual representation of this algebra's multiplication
487 table (on basis elements).
491 sage: from mjo.eja.eja_algebra import JordanSpinEJA
495 sage: J = JordanSpinEJA(4)
496 sage: J.multiplication_table()
497 +----++----+----+----+----+
498 | * || e0 | e1 | e2 | e3 |
499 +====++====+====+====+====+
500 | e0 || e0 | e1 | e2 | e3 |
501 +----++----+----+----+----+
502 | e1 || e1 | e0 | 0 | 0 |
503 +----++----+----+----+----+
504 | e2 || e2 | 0 | e0 | 0 |
505 +----++----+----+----+----+
506 | e3 || e3 | 0 | 0 | e0 |
507 +----++----+----+----+----+
510 M
= list(self
._multiplication
_table
) # copy
511 for i
in range(len(M
)):
512 # M had better be "square"
513 M
[i
] = [self
.monomial(i
)] + M
[i
]
514 M
= [["*"] + list(self
.gens())] + M
515 return table(M
, header_row
=True, header_column
=True, frame
=True)
518 def natural_basis(self
):
520 Return a more-natural representation of this algebra's basis.
522 Every finite-dimensional Euclidean Jordan Algebra is a direct
523 sum of five simple algebras, four of which comprise Hermitian
524 matrices. This method returns the original "natural" basis
525 for our underlying vector space. (Typically, the natural basis
526 is used to construct the multiplication table in the first place.)
528 Note that this will always return a matrix. The standard basis
529 in `R^n` will be returned as `n`-by-`1` column matrices.
533 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
534 ....: RealSymmetricEJA)
538 sage: J = RealSymmetricEJA(2)
540 Finite family {0: e0, 1: e1, 2: e2}
541 sage: J.natural_basis()
543 [1 0] [ 0 0.7071067811865475?] [0 0]
544 [0 0], [0.7071067811865475? 0], [0 1]
549 sage: J = JordanSpinEJA(2)
551 Finite family {0: e0, 1: e1}
552 sage: J.natural_basis()
559 if self
._natural
_basis
is None:
560 M
= self
.natural_basis_space()
561 return tuple( M(b
.to_vector()) for b
in self
.basis() )
563 return self
._natural
_basis
566 def natural_basis_space(self
):
568 Return the matrix space in which this algebra's natural basis
571 Generally this will be an `n`-by-`1` column-vector space,
572 except when the algebra is trivial. There it's `n`-by-`n`
573 (where `n` is zero), to ensure that two elements of the
574 natural basis space (empty matrices) can be multiplied.
576 if self
.is_trivial():
577 return MatrixSpace(self
.base_ring(), 0)
578 elif self
._natural
_basis
is None or len(self
._natural
_basis
) == 0:
579 return MatrixSpace(self
.base_ring(), self
.dimension(), 1)
581 return self
._natural
_basis
[0].matrix_space()
587 Return the unit element of this algebra.
591 sage: from mjo.eja.eja_algebra import (HadamardEJA,
596 sage: J = HadamardEJA(5)
598 e0 + e1 + e2 + e3 + e4
602 The identity element acts like the identity::
604 sage: set_random_seed()
605 sage: J = random_eja()
606 sage: x = J.random_element()
607 sage: J.one()*x == x and x*J.one() == x
610 The matrix of the unit element's operator is the identity::
612 sage: set_random_seed()
613 sage: J = random_eja()
614 sage: actual = J.one().operator().matrix()
615 sage: expected = matrix.identity(J.base_ring(), J.dimension())
616 sage: actual == expected
619 Ensure that the cached unit element (often precomputed by
620 hand) agrees with the computed one::
622 sage: set_random_seed()
623 sage: J = random_eja()
624 sage: cached = J.one()
625 sage: J.one.clear_cache()
626 sage: J.one() == cached
630 # We can brute-force compute the matrices of the operators
631 # that correspond to the basis elements of this algebra.
632 # If some linear combination of those basis elements is the
633 # algebra identity, then the same linear combination of
634 # their matrices has to be the identity matrix.
636 # Of course, matrices aren't vectors in sage, so we have to
637 # appeal to the "long vectors" isometry.
638 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
640 # Now we use basic linear algebra to find the coefficients,
641 # of the matrices-as-vectors-linear-combination, which should
642 # work for the original algebra basis too.
643 A
= matrix(self
.base_ring(), oper_vecs
)
645 # We used the isometry on the left-hand side already, but we
646 # still need to do it for the right-hand side. Recall that we
647 # wanted something that summed to the identity matrix.
648 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
650 # Now if there's an identity element in the algebra, this
651 # should work. We solve on the left to avoid having to
652 # transpose the matrix "A".
653 return self
.from_vector(A
.solve_left(b
))
656 def peirce_decomposition(self
, c
):
658 The Peirce decomposition of this algebra relative to the
661 In the future, this can be extended to a complete system of
662 orthogonal idempotents.
666 - ``c`` -- an idempotent of this algebra.
670 A triple (J0, J5, J1) containing two subalgebras and one subspace
673 - ``J0`` -- the algebra on the eigenspace of ``c.operator()``
674 corresponding to the eigenvalue zero.
676 - ``J5`` -- the eigenspace (NOT a subalgebra) of ``c.operator()``
677 corresponding to the eigenvalue one-half.
679 - ``J1`` -- the algebra on the eigenspace of ``c.operator()``
680 corresponding to the eigenvalue one.
682 These are the only possible eigenspaces for that operator, and this
683 algebra is a direct sum of them. The spaces ``J0`` and ``J1`` are
684 orthogonal, and are subalgebras of this algebra with the appropriate
689 sage: from mjo.eja.eja_algebra import random_eja, RealSymmetricEJA
693 The canonical example comes from the symmetric matrices, which
694 decompose into diagonal and off-diagonal parts::
696 sage: J = RealSymmetricEJA(3)
697 sage: C = matrix(QQ, [ [1,0,0],
701 sage: J0,J5,J1 = J.peirce_decomposition(c)
703 Euclidean Jordan algebra of dimension 1...
705 Vector space of degree 6 and dimension 2...
707 Euclidean Jordan algebra of dimension 3...
708 sage: J0.one().natural_representation()
712 sage: orig_df = AA.options.display_format
713 sage: AA.options.display_format = 'radical'
714 sage: J.from_vector(J5.basis()[0]).natural_representation()
718 sage: J.from_vector(J5.basis()[1]).natural_representation()
722 sage: AA.options.display_format = orig_df
723 sage: J1.one().natural_representation()
730 Every algebra decomposes trivially with respect to its identity
733 sage: set_random_seed()
734 sage: J = random_eja()
735 sage: J0,J5,J1 = J.peirce_decomposition(J.one())
736 sage: J0.dimension() == 0 and J5.dimension() == 0
738 sage: J1.superalgebra() == J and J1.dimension() == J.dimension()
741 The decomposition is into eigenspaces, and its components are
742 therefore necessarily orthogonal. Moreover, the identity
743 elements in the two subalgebras are the projections onto their
744 respective subspaces of the superalgebra's identity element::
746 sage: set_random_seed()
747 sage: J = random_eja()
748 sage: x = J.random_element()
749 sage: if not J.is_trivial():
750 ....: while x.is_nilpotent():
751 ....: x = J.random_element()
752 sage: c = x.subalgebra_idempotent()
753 sage: J0,J5,J1 = J.peirce_decomposition(c)
755 sage: for (w,y,z) in zip(J0.basis(), J5.basis(), J1.basis()):
756 ....: w = w.superalgebra_element()
757 ....: y = J.from_vector(y)
758 ....: z = z.superalgebra_element()
759 ....: ipsum += w.inner_product(y).abs()
760 ....: ipsum += w.inner_product(z).abs()
761 ....: ipsum += y.inner_product(z).abs()
764 sage: J1(c) == J1.one()
766 sage: J0(J.one() - c) == J0.one()
770 if not c
.is_idempotent():
771 raise ValueError("element is not idempotent: %s" % c
)
773 from mjo
.eja
.eja_subalgebra
import FiniteDimensionalEuclideanJordanSubalgebra
775 # Default these to what they should be if they turn out to be
776 # trivial, because eigenspaces_left() won't return eigenvalues
777 # corresponding to trivial spaces (e.g. it returns only the
778 # eigenspace corresponding to lambda=1 if you take the
779 # decomposition relative to the identity element).
780 trivial
= FiniteDimensionalEuclideanJordanSubalgebra(self
, ())
781 J0
= trivial
# eigenvalue zero
782 J5
= VectorSpace(self
.base_ring(), 0) # eigenvalue one-half
783 J1
= trivial
# eigenvalue one
785 for (eigval
, eigspace
) in c
.operator().matrix().right_eigenspaces():
786 if eigval
== ~
(self
.base_ring()(2)):
789 gens
= tuple( self
.from_vector(b
) for b
in eigspace
.basis() )
790 subalg
= FiniteDimensionalEuclideanJordanSubalgebra(self
,
798 raise ValueError("unexpected eigenvalue: %s" % eigval
)
803 def random_element(self
, thorough
=False):
805 Return a random element of this algebra.
807 Our algebra superclass method only returns a linear
808 combination of at most two basis elements. We instead
809 want the vector space "random element" method that
810 returns a more diverse selection.
814 - ``thorough`` -- (boolean; default False) whether or not we
815 should generate irrational coefficients for the random
816 element when our base ring is irrational; this slows the
817 algebra operations to a crawl, but any truly random method
821 # For a general base ring... maybe we can trust this to do the
822 # right thing? Unlikely, but.
823 V
= self
.vector_space()
824 v
= V
.random_element()
826 if self
.base_ring() is AA
:
827 # The "random element" method of the algebraic reals is
828 # stupid at the moment, and only returns integers between
829 # -2 and 2, inclusive:
831 # https://trac.sagemath.org/ticket/30875
833 # Instead, we implement our own "random vector" method,
834 # and then coerce that into the algebra. We use the vector
835 # space degree here instead of the dimension because a
836 # subalgebra could (for example) be spanned by only two
837 # vectors, each with five coordinates. We need to
838 # generate all five coordinates.
840 v
*= QQbar
.random_element().real()
842 v
*= QQ
.random_element()
844 return self
.from_vector(V
.coordinate_vector(v
))
846 def random_elements(self
, count
, thorough
=False):
848 Return ``count`` random elements as a tuple.
852 - ``thorough`` -- (boolean; default False) whether or not we
853 should generate irrational coefficients for the random
854 elements when our base ring is irrational; this slows the
855 algebra operations to a crawl, but any truly random method
860 sage: from mjo.eja.eja_algebra import JordanSpinEJA
864 sage: J = JordanSpinEJA(3)
865 sage: x,y,z = J.random_elements(3)
866 sage: all( [ x in J, y in J, z in J ])
868 sage: len( J.random_elements(10) ) == 10
872 return tuple( self
.random_element(thorough
)
873 for idx
in range(count
) )
877 def _charpoly_coefficients(self
):
879 The `r` polynomial coefficients of the "characteristic polynomial
883 R
= self
.coordinate_polynomial_ring()
885 F
= R
.fraction_field()
888 # From a result in my book, these are the entries of the
889 # basis representation of L_x.
890 return sum( vars[k
]*self
.monomial(k
).operator().matrix()[i
,j
]
893 L_x
= matrix(F
, n
, n
, L_x_i_j
)
896 if self
.rank
.is_in_cache():
898 # There's no need to pad the system with redundant
899 # columns if we *know* they'll be redundant.
902 # Compute an extra power in case the rank is equal to
903 # the dimension (otherwise, we would stop at x^(r-1)).
904 x_powers
= [ (L_x
**k
)*self
.one().to_vector()
905 for k
in range(n
+1) ]
906 A
= matrix
.column(F
, x_powers
[:n
])
907 AE
= A
.extended_echelon_form()
914 # The theory says that only the first "r" coefficients are
915 # nonzero, and they actually live in the original polynomial
916 # ring and not the fraction field. We negate them because
917 # in the actual characteristic polynomial, they get moved
918 # to the other side where x^r lives.
919 return -A_rref
.solve_right(E
*b
).change_ring(R
)[:r
]
924 Return the rank of this EJA.
926 This is a cached method because we know the rank a priori for
927 all of the algebras we can construct. Thus we can avoid the
928 expensive ``_charpoly_coefficients()`` call unless we truly
929 need to compute the whole characteristic polynomial.
933 sage: from mjo.eja.eja_algebra import (HadamardEJA,
935 ....: RealSymmetricEJA,
936 ....: ComplexHermitianEJA,
937 ....: QuaternionHermitianEJA,
942 The rank of the Jordan spin algebra is always two::
944 sage: JordanSpinEJA(2).rank()
946 sage: JordanSpinEJA(3).rank()
948 sage: JordanSpinEJA(4).rank()
951 The rank of the `n`-by-`n` Hermitian real, complex, or
952 quaternion matrices is `n`::
954 sage: RealSymmetricEJA(4).rank()
956 sage: ComplexHermitianEJA(3).rank()
958 sage: QuaternionHermitianEJA(2).rank()
963 Ensure that every EJA that we know how to construct has a
964 positive integer rank, unless the algebra is trivial in
965 which case its rank will be zero::
967 sage: set_random_seed()
968 sage: J = random_eja()
972 sage: r > 0 or (r == 0 and J.is_trivial())
975 Ensure that computing the rank actually works, since the ranks
976 of all simple algebras are known and will be cached by default::
978 sage: J = HadamardEJA(4)
979 sage: J.rank.clear_cache()
985 sage: J = JordanSpinEJA(4)
986 sage: J.rank.clear_cache()
992 sage: J = RealSymmetricEJA(3)
993 sage: J.rank.clear_cache()
999 sage: J = ComplexHermitianEJA(2)
1000 sage: J.rank.clear_cache()
1006 sage: J = QuaternionHermitianEJA(2)
1007 sage: J.rank.clear_cache()
1011 return len(self
._charpoly
_coefficients
())
1014 def vector_space(self
):
1016 Return the vector space that underlies this algebra.
1020 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1024 sage: J = RealSymmetricEJA(2)
1025 sage: J.vector_space()
1026 Vector space of dimension 3 over...
1029 return self
.zero().to_vector().parent().ambient_vector_space()
1032 Element
= FiniteDimensionalEuclideanJordanAlgebraElement
1035 class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra
):
1037 Algebras whose basis consists of vectors with rational
1038 entries. Equivalently, algebras whose multiplication tables
1039 contain only rational coefficients.
1041 When an EJA has a basis that can be made rational, we can speed up
1042 the computation of its characteristic polynomial by doing it over
1043 ``QQ``. All of the named EJA constructors that we provide fall
1047 def _charpoly_coefficients(self
):
1049 Override the parent method with something that tries to compute
1050 over a faster (non-extension) field.
1054 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1058 The base ring of the resulting polynomial coefficients is what
1059 it should be, and not the rationals (unless the algebra was
1060 already over the rationals)::
1062 sage: J = JordanSpinEJA(3)
1063 sage: J._charpoly_coefficients()
1064 (X1^2 - X2^2 - X3^2, -2*X1)
1065 sage: a0 = J._charpoly_coefficients()[0]
1067 Algebraic Real Field
1068 sage: a0.base_ring()
1069 Algebraic Real Field
1072 if self
.base_ring() is QQ
:
1073 # There's no need to construct *another* algebra over the
1074 # rationals if this one is already over the rationals.
1075 superclass
= super(RationalBasisEuclideanJordanAlgebra
, self
)
1076 return superclass
._charpoly
_coefficients
()
1079 map(lambda x
: x
.to_vector(), ls
)
1080 for ls
in self
._multiplication
_table
1083 # Do the computation over the rationals. The answer will be
1084 # the same, because our basis coordinates are (essentially)
1086 J
= FiniteDimensionalEuclideanJordanAlgebra(QQ
,
1090 a
= J
._charpoly
_coefficients
()
1091 return tuple(map(lambda x
: x
.change_ring(self
.base_ring()), a
))
1094 class ConcreteEuclideanJordanAlgebra
:
1096 A class for the Euclidean Jordan algebras that we know by name.
1098 These are the Jordan algebras whose basis, multiplication table,
1099 rank, and so on are known a priori. More to the point, they are
1100 the Euclidean Jordan algebras for which we are able to conjure up
1101 a "random instance."
1105 sage: from mjo.eja.eja_algebra import ConcreteEuclideanJordanAlgebra
1109 Our natural basis is normalized with respect to the natural inner
1110 product unless we specify otherwise::
1112 sage: set_random_seed()
1113 sage: J = ConcreteEuclideanJordanAlgebra.random_instance()
1114 sage: all( b.norm() == 1 for b in J.gens() )
1117 Since our natural basis is normalized with respect to the natural
1118 inner product, and since we know that this algebra is an EJA, any
1119 left-multiplication operator's matrix will be symmetric because
1120 natural->EJA basis representation is an isometry and within the EJA
1121 the operator is self-adjoint by the Jordan axiom::
1123 sage: set_random_seed()
1124 sage: J = ConcreteEuclideanJordanAlgebra.random_instance()
1125 sage: x = J.random_element()
1126 sage: x.operator().matrix().is_symmetric()
1132 def _max_random_instance_size():
1134 Return an integer "size" that is an upper bound on the size of
1135 this algebra when it is used in a random test
1136 case. Unfortunately, the term "size" is ambiguous -- when
1137 dealing with `R^n` under either the Hadamard or Jordan spin
1138 product, the "size" refers to the dimension `n`. When dealing
1139 with a matrix algebra (real symmetric or complex/quaternion
1140 Hermitian), it refers to the size of the matrix, which is far
1141 less than the dimension of the underlying vector space.
1143 This method must be implemented in each subclass.
1145 raise NotImplementedError
1148 def random_instance(cls
, field
=AA
, **kwargs
):
1150 Return a random instance of this type of algebra.
1152 This method should be implemented in each subclass.
1154 from sage
.misc
.prandom
import choice
1155 eja_class
= choice(cls
.__subclasses
__())
1156 return eja_class
.random_instance(field
)
1159 class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra
):
1161 def __init__(self
, field
, basis
, normalize_basis
=True, **kwargs
):
1163 Compared to the superclass constructor, we take a basis instead of
1164 a multiplication table because the latter can be computed in terms
1165 of the former when the product is known (like it is here).
1167 # Used in this class's fast _charpoly_coefficients() override.
1168 self
._basis
_normalizers
= None
1170 # We're going to loop through this a few times, so now's a good
1171 # time to ensure that it isn't a generator expression.
1172 basis
= tuple(basis
)
1174 algebra_dim
= len(basis
)
1175 degree
= 0 # size of the matrices
1177 degree
= basis
[0].nrows()
1179 if algebra_dim
> 1 and normalize_basis
:
1180 # We'll need sqrt(2) to normalize the basis, and this
1181 # winds up in the multiplication table, so the whole
1182 # algebra needs to be over the field extension.
1183 R
= PolynomialRing(field
, 'z')
1186 if p
.is_irreducible():
1187 field
= field
.extension(p
, 'sqrt2', embedding
=RLF(2).sqrt())
1188 basis
= tuple( s
.change_ring(field
) for s
in basis
)
1189 self
._basis
_normalizers
= tuple(
1190 ~
(self
.natural_inner_product(s
,s
).sqrt()) for s
in basis
)
1191 basis
= tuple(s
*c
for (s
,c
) in zip(basis
,self
._basis
_normalizers
))
1193 # Now compute the multiplication and inner product tables.
1194 # We have to do this *after* normalizing the basis, because
1195 # scaling affects the answers.
1196 V
= VectorSpace(field
, degree
**2)
1197 W
= V
.span_of_basis( _mat2vec(s
) for s
in basis
)
1198 mult_table
= [[W
.zero() for j
in range(algebra_dim
)]
1199 for i
in range(algebra_dim
)]
1200 ip_table
= [[W
.zero() for j
in range(algebra_dim
)]
1201 for i
in range(algebra_dim
)]
1202 for i
in range(algebra_dim
):
1203 for j
in range(algebra_dim
):
1204 mat_entry
= (basis
[i
]*basis
[j
] + basis
[j
]*basis
[i
])/2
1205 mult_table
[i
][j
] = W
.coordinate_vector(_mat2vec(mat_entry
))
1208 # HACK: ignore the error here if we don't need the
1209 # inner product (as is the case when we construct
1210 # a dummy QQ-algebra for fast charpoly coefficients.
1211 ip_table
[i
][j
] = self
.natural_inner_product(basis
[i
],
1218 self
._inner
_product
_matrix
= matrix(field
,ip_table
)
1222 super(MatrixEuclideanJordanAlgebra
, self
).__init
__(field
,
1224 natural_basis
=basis
,
1227 if algebra_dim
== 0:
1228 self
.one
.set_cache(self
.zero())
1230 n
= basis
[0].nrows()
1231 # The identity wrt (A,B) -> (AB + BA)/2 is independent of the
1232 # details of this algebra.
1233 self
.one
.set_cache(self(matrix
.identity(field
,n
)))
1237 def _charpoly_coefficients(self
):
1239 Override the parent method with something that tries to compute
1240 over a faster (non-extension) field.
1242 if self
._basis
_normalizers
is None or self
.base_ring() is QQ
:
1243 # We didn't normalize, or the basis we started with had
1244 # entries in a nice field already. Just compute the thing.
1245 return super(MatrixEuclideanJordanAlgebra
, self
)._charpoly
_coefficients
()
1247 basis
= ( (b
/n
) for (b
,n
) in zip(self
.natural_basis(),
1248 self
._basis
_normalizers
) )
1250 # Do this over the rationals and convert back at the end.
1251 # Only works because we know the entries of the basis are
1252 # integers. The argument ``check_axioms=False`` is required
1253 # because the trace inner-product method for this
1254 # class is a stub and can't actually be checked.
1255 J
= MatrixEuclideanJordanAlgebra(QQ
,
1257 normalize_basis
=False,
1260 a
= J
._charpoly
_coefficients
()
1262 # Unfortunately, changing the basis does change the
1263 # coefficients of the characteristic polynomial, but since
1264 # these are really the coefficients of the "characteristic
1265 # polynomial of" function, everything is still nice and
1266 # unevaluated. It's therefore "obvious" how scaling the
1267 # basis affects the coordinate variables X1, X2, et
1268 # cetera. Scaling the first basis vector up by "n" adds a
1269 # factor of 1/n into every "X1" term, for example. So here
1270 # we simply undo the basis_normalizer scaling that we
1271 # performed earlier.
1273 # The a[0] access here is safe because trivial algebras
1274 # won't have any basis normalizers and therefore won't
1275 # make it to this "else" branch.
1276 XS
= a
[0].parent().gens()
1277 subs_dict
= { XS
[i
]: self
._basis
_normalizers
[i
]*XS
[i
]
1278 for i
in range(len(XS
)) }
1279 return tuple( a_i
.subs(subs_dict
) for a_i
in a
)
1285 Embed the matrix ``M`` into a space of real matrices.
1287 The matrix ``M`` can have entries in any field at the moment:
1288 the real numbers, complex numbers, or quaternions. And although
1289 they are not a field, we can probably support octonions at some
1290 point, too. This function returns a real matrix that "acts like"
1291 the original with respect to matrix multiplication; i.e.
1293 real_embed(M*N) = real_embed(M)*real_embed(N)
1296 raise NotImplementedError
1300 def real_unembed(M
):
1302 The inverse of :meth:`real_embed`.
1304 raise NotImplementedError
1307 def natural_inner_product(cls
,X
,Y
):
1308 Xu
= cls
.real_unembed(X
)
1309 Yu
= cls
.real_unembed(Y
)
1310 tr
= (Xu
*Yu
).trace()
1313 # Works in QQ, AA, RDF, et cetera.
1315 except AttributeError:
1316 # A quaternion doesn't have a real() method, but does
1317 # have coefficient_tuple() method that returns the
1318 # coefficients of 1, i, j, and k -- in that order.
1319 return tr
.coefficient_tuple()[0]
1322 class RealMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1326 The identity function, for embedding real matrices into real
1332 def real_unembed(M
):
1334 The identity function, for unembedding real matrices from real
1340 class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra
,
1341 ConcreteEuclideanJordanAlgebra
):
1343 The rank-n simple EJA consisting of real symmetric n-by-n
1344 matrices, the usual symmetric Jordan product, and the trace inner
1345 product. It has dimension `(n^2 + n)/2` over the reals.
1349 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1353 sage: J = RealSymmetricEJA(2)
1354 sage: e0, e1, e2 = J.gens()
1362 In theory, our "field" can be any subfield of the reals::
1364 sage: RealSymmetricEJA(2, RDF)
1365 Euclidean Jordan algebra of dimension 3 over Real Double Field
1366 sage: RealSymmetricEJA(2, RR)
1367 Euclidean Jordan algebra of dimension 3 over Real Field with
1368 53 bits of precision
1372 The dimension of this algebra is `(n^2 + n) / 2`::
1374 sage: set_random_seed()
1375 sage: n_max = RealSymmetricEJA._max_random_instance_size()
1376 sage: n = ZZ.random_element(1, n_max)
1377 sage: J = RealSymmetricEJA(n)
1378 sage: J.dimension() == (n^2 + n)/2
1381 The Jordan multiplication is what we think it is::
1383 sage: set_random_seed()
1384 sage: J = RealSymmetricEJA.random_instance()
1385 sage: x,y = J.random_elements(2)
1386 sage: actual = (x*y).natural_representation()
1387 sage: X = x.natural_representation()
1388 sage: Y = y.natural_representation()
1389 sage: expected = (X*Y + Y*X)/2
1390 sage: actual == expected
1392 sage: J(expected) == x*y
1395 We can change the generator prefix::
1397 sage: RealSymmetricEJA(3, prefix='q').gens()
1398 (q0, q1, q2, q3, q4, q5)
1400 We can construct the (trivial) algebra of rank zero::
1402 sage: RealSymmetricEJA(0)
1403 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1407 def _denormalized_basis(cls
, n
, field
):
1409 Return a basis for the space of real symmetric n-by-n matrices.
1413 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1417 sage: set_random_seed()
1418 sage: n = ZZ.random_element(1,5)
1419 sage: B = RealSymmetricEJA._denormalized_basis(n,QQ)
1420 sage: all( M.is_symmetric() for M in B)
1424 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1428 for j
in range(i
+1):
1429 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
1433 Sij
= Eij
+ Eij
.transpose()
1439 def _max_random_instance_size():
1440 return 4 # Dimension 10
1443 def random_instance(cls
, field
=AA
, **kwargs
):
1445 Return a random instance of this type of algebra.
1447 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1448 return cls(n
, field
, **kwargs
)
1450 def __init__(self
, n
, field
=AA
, **kwargs
):
1451 basis
= self
._denormalized
_basis
(n
, field
)
1452 super(RealSymmetricEJA
, self
).__init
__(field
,
1456 self
.rank
.set_cache(n
)
1459 class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1463 Embed the n-by-n complex matrix ``M`` into the space of real
1464 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1465 bi` to the block matrix ``[[a,b],[-b,a]]``.
1469 sage: from mjo.eja.eja_algebra import \
1470 ....: ComplexMatrixEuclideanJordanAlgebra
1474 sage: F = QuadraticField(-1, 'I')
1475 sage: x1 = F(4 - 2*i)
1476 sage: x2 = F(1 + 2*i)
1479 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1480 sage: ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1489 Embedding is a homomorphism (isomorphism, in fact)::
1491 sage: set_random_seed()
1492 sage: n = ZZ.random_element(3)
1493 sage: F = QuadraticField(-1, 'I')
1494 sage: X = random_matrix(F, n)
1495 sage: Y = random_matrix(F, n)
1496 sage: Xe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X)
1497 sage: Ye = ComplexMatrixEuclideanJordanAlgebra.real_embed(Y)
1498 sage: XYe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1505 raise ValueError("the matrix 'M' must be square")
1507 # We don't need any adjoined elements...
1508 field
= M
.base_ring().base_ring()
1512 a
= z
.list()[0] # real part, I guess
1513 b
= z
.list()[1] # imag part, I guess
1514 blocks
.append(matrix(field
, 2, [[a
,b
],[-b
,a
]]))
1516 return matrix
.block(field
, n
, blocks
)
1520 def real_unembed(M
):
1522 The inverse of _embed_complex_matrix().
1526 sage: from mjo.eja.eja_algebra import \
1527 ....: ComplexMatrixEuclideanJordanAlgebra
1531 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1532 ....: [-2, 1, -4, 3],
1533 ....: [ 9, 10, 11, 12],
1534 ....: [-10, 9, -12, 11] ])
1535 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(A)
1537 [ 10*I + 9 12*I + 11]
1541 Unembedding is the inverse of embedding::
1543 sage: set_random_seed()
1544 sage: F = QuadraticField(-1, 'I')
1545 sage: M = random_matrix(F, 3)
1546 sage: Me = ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1547 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1553 raise ValueError("the matrix 'M' must be square")
1554 if not n
.mod(2).is_zero():
1555 raise ValueError("the matrix 'M' must be a complex embedding")
1557 # If "M" was normalized, its base ring might have roots
1558 # adjoined and they can stick around after unembedding.
1559 field
= M
.base_ring()
1560 R
= PolynomialRing(field
, 'z')
1563 # Sage doesn't know how to embed AA into QQbar, i.e. how
1564 # to adjoin sqrt(-1) to AA.
1567 F
= field
.extension(z
**2 + 1, 'I', embedding
=CLF(-1).sqrt())
1570 # Go top-left to bottom-right (reading order), converting every
1571 # 2-by-2 block we see to a single complex element.
1573 for k
in range(n
/2):
1574 for j
in range(n
/2):
1575 submat
= M
[2*k
:2*k
+2,2*j
:2*j
+2]
1576 if submat
[0,0] != submat
[1,1]:
1577 raise ValueError('bad on-diagonal submatrix')
1578 if submat
[0,1] != -submat
[1,0]:
1579 raise ValueError('bad off-diagonal submatrix')
1580 z
= submat
[0,0] + submat
[0,1]*i
1583 return matrix(F
, n
/2, elements
)
1587 def natural_inner_product(cls
,X
,Y
):
1589 Compute a natural inner product in this algebra directly from
1594 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1598 This gives the same answer as the slow, default method implemented
1599 in :class:`MatrixEuclideanJordanAlgebra`::
1601 sage: set_random_seed()
1602 sage: J = ComplexHermitianEJA.random_instance()
1603 sage: x,y = J.random_elements(2)
1604 sage: Xe = x.natural_representation()
1605 sage: Ye = y.natural_representation()
1606 sage: X = ComplexHermitianEJA.real_unembed(Xe)
1607 sage: Y = ComplexHermitianEJA.real_unembed(Ye)
1608 sage: expected = (X*Y).trace().real()
1609 sage: actual = ComplexHermitianEJA.natural_inner_product(Xe,Ye)
1610 sage: actual == expected
1614 return RealMatrixEuclideanJordanAlgebra
.natural_inner_product(X
,Y
)/2
1617 class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra
,
1618 ConcreteEuclideanJordanAlgebra
):
1620 The rank-n simple EJA consisting of complex Hermitian n-by-n
1621 matrices over the real numbers, the usual symmetric Jordan product,
1622 and the real-part-of-trace inner product. It has dimension `n^2` over
1627 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1631 In theory, our "field" can be any subfield of the reals::
1633 sage: ComplexHermitianEJA(2, RDF)
1634 Euclidean Jordan algebra of dimension 4 over Real Double Field
1635 sage: ComplexHermitianEJA(2, RR)
1636 Euclidean Jordan algebra of dimension 4 over Real Field with
1637 53 bits of precision
1641 The dimension of this algebra is `n^2`::
1643 sage: set_random_seed()
1644 sage: n_max = ComplexHermitianEJA._max_random_instance_size()
1645 sage: n = ZZ.random_element(1, n_max)
1646 sage: J = ComplexHermitianEJA(n)
1647 sage: J.dimension() == n^2
1650 The Jordan multiplication is what we think it is::
1652 sage: set_random_seed()
1653 sage: J = ComplexHermitianEJA.random_instance()
1654 sage: x,y = J.random_elements(2)
1655 sage: actual = (x*y).natural_representation()
1656 sage: X = x.natural_representation()
1657 sage: Y = y.natural_representation()
1658 sage: expected = (X*Y + Y*X)/2
1659 sage: actual == expected
1661 sage: J(expected) == x*y
1664 We can change the generator prefix::
1666 sage: ComplexHermitianEJA(2, prefix='z').gens()
1669 We can construct the (trivial) algebra of rank zero::
1671 sage: ComplexHermitianEJA(0)
1672 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1677 def _denormalized_basis(cls
, n
, field
):
1679 Returns a basis for the space of complex Hermitian n-by-n matrices.
1681 Why do we embed these? Basically, because all of numerical linear
1682 algebra assumes that you're working with vectors consisting of `n`
1683 entries from a field and scalars from the same field. There's no way
1684 to tell SageMath that (for example) the vectors contain complex
1685 numbers, while the scalar field is real.
1689 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1693 sage: set_random_seed()
1694 sage: n = ZZ.random_element(1,5)
1695 sage: field = QuadraticField(2, 'sqrt2')
1696 sage: B = ComplexHermitianEJA._denormalized_basis(n, field)
1697 sage: all( M.is_symmetric() for M in B)
1701 R
= PolynomialRing(field
, 'z')
1703 F
= field
.extension(z
**2 + 1, 'I')
1706 # This is like the symmetric case, but we need to be careful:
1708 # * We want conjugate-symmetry, not just symmetry.
1709 # * The diagonal will (as a result) be real.
1713 for j
in range(i
+1):
1714 Eij
= matrix(F
, n
, lambda k
,l
: k
==i
and l
==j
)
1716 Sij
= cls
.real_embed(Eij
)
1719 # The second one has a minus because it's conjugated.
1720 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
1722 Sij_imag
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
1725 # Since we embedded these, we can drop back to the "field" that we
1726 # started with instead of the complex extension "F".
1727 return ( s
.change_ring(field
) for s
in S
)
1730 def __init__(self
, n
, field
=AA
, **kwargs
):
1731 basis
= self
._denormalized
_basis
(n
,field
)
1732 super(ComplexHermitianEJA
,self
).__init
__(field
,
1736 self
.rank
.set_cache(n
)
1739 def _max_random_instance_size():
1740 return 3 # Dimension 9
1743 def random_instance(cls
, field
=AA
, **kwargs
):
1745 Return a random instance of this type of algebra.
1747 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1748 return cls(n
, field
, **kwargs
)
1750 class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1754 Embed the n-by-n quaternion matrix ``M`` into the space of real
1755 matrices of size 4n-by-4n by first sending each quaternion entry `z
1756 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
1757 c+di],[-c + di, a-bi]]`, and then embedding those into a real
1762 sage: from mjo.eja.eja_algebra import \
1763 ....: QuaternionMatrixEuclideanJordanAlgebra
1767 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1768 sage: i,j,k = Q.gens()
1769 sage: x = 1 + 2*i + 3*j + 4*k
1770 sage: M = matrix(Q, 1, [[x]])
1771 sage: QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1777 Embedding is a homomorphism (isomorphism, in fact)::
1779 sage: set_random_seed()
1780 sage: n = ZZ.random_element(2)
1781 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1782 sage: X = random_matrix(Q, n)
1783 sage: Y = random_matrix(Q, n)
1784 sage: Xe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X)
1785 sage: Ye = QuaternionMatrixEuclideanJordanAlgebra.real_embed(Y)
1786 sage: XYe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1791 quaternions
= M
.base_ring()
1794 raise ValueError("the matrix 'M' must be square")
1796 F
= QuadraticField(-1, 'I')
1801 t
= z
.coefficient_tuple()
1806 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
1807 [-c
+ d
*i
, a
- b
*i
]])
1808 realM
= ComplexMatrixEuclideanJordanAlgebra
.real_embed(cplxM
)
1809 blocks
.append(realM
)
1811 # We should have real entries by now, so use the realest field
1812 # we've got for the return value.
1813 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
1818 def real_unembed(M
):
1820 The inverse of _embed_quaternion_matrix().
1824 sage: from mjo.eja.eja_algebra import \
1825 ....: QuaternionMatrixEuclideanJordanAlgebra
1829 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
1830 ....: [-2, 1, -4, 3],
1831 ....: [-3, 4, 1, -2],
1832 ....: [-4, -3, 2, 1]])
1833 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(M)
1834 [1 + 2*i + 3*j + 4*k]
1838 Unembedding is the inverse of embedding::
1840 sage: set_random_seed()
1841 sage: Q = QuaternionAlgebra(QQ, -1, -1)
1842 sage: M = random_matrix(Q, 3)
1843 sage: Me = QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1844 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1850 raise ValueError("the matrix 'M' must be square")
1851 if not n
.mod(4).is_zero():
1852 raise ValueError("the matrix 'M' must be a quaternion embedding")
1854 # Use the base ring of the matrix to ensure that its entries can be
1855 # multiplied by elements of the quaternion algebra.
1856 field
= M
.base_ring()
1857 Q
= QuaternionAlgebra(field
,-1,-1)
1860 # Go top-left to bottom-right (reading order), converting every
1861 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
1864 for l
in range(n
/4):
1865 for m
in range(n
/4):
1866 submat
= ComplexMatrixEuclideanJordanAlgebra
.real_unembed(
1867 M
[4*l
:4*l
+4,4*m
:4*m
+4] )
1868 if submat
[0,0] != submat
[1,1].conjugate():
1869 raise ValueError('bad on-diagonal submatrix')
1870 if submat
[0,1] != -submat
[1,0].conjugate():
1871 raise ValueError('bad off-diagonal submatrix')
1872 z
= submat
[0,0].real()
1873 z
+= submat
[0,0].imag()*i
1874 z
+= submat
[0,1].real()*j
1875 z
+= submat
[0,1].imag()*k
1878 return matrix(Q
, n
/4, elements
)
1882 def natural_inner_product(cls
,X
,Y
):
1884 Compute a natural inner product in this algebra directly from
1889 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1893 This gives the same answer as the slow, default method implemented
1894 in :class:`MatrixEuclideanJordanAlgebra`::
1896 sage: set_random_seed()
1897 sage: J = QuaternionHermitianEJA.random_instance()
1898 sage: x,y = J.random_elements(2)
1899 sage: Xe = x.natural_representation()
1900 sage: Ye = y.natural_representation()
1901 sage: X = QuaternionHermitianEJA.real_unembed(Xe)
1902 sage: Y = QuaternionHermitianEJA.real_unembed(Ye)
1903 sage: expected = (X*Y).trace().coefficient_tuple()[0]
1904 sage: actual = QuaternionHermitianEJA.natural_inner_product(Xe,Ye)
1905 sage: actual == expected
1909 return RealMatrixEuclideanJordanAlgebra
.natural_inner_product(X
,Y
)/4
1912 class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra
,
1913 ConcreteEuclideanJordanAlgebra
):
1915 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
1916 matrices, the usual symmetric Jordan product, and the
1917 real-part-of-trace inner product. It has dimension `2n^2 - n` over
1922 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1926 In theory, our "field" can be any subfield of the reals::
1928 sage: QuaternionHermitianEJA(2, RDF)
1929 Euclidean Jordan algebra of dimension 6 over Real Double Field
1930 sage: QuaternionHermitianEJA(2, RR)
1931 Euclidean Jordan algebra of dimension 6 over Real Field with
1932 53 bits of precision
1936 The dimension of this algebra is `2*n^2 - n`::
1938 sage: set_random_seed()
1939 sage: n_max = QuaternionHermitianEJA._max_random_instance_size()
1940 sage: n = ZZ.random_element(1, n_max)
1941 sage: J = QuaternionHermitianEJA(n)
1942 sage: J.dimension() == 2*(n^2) - n
1945 The Jordan multiplication is what we think it is::
1947 sage: set_random_seed()
1948 sage: J = QuaternionHermitianEJA.random_instance()
1949 sage: x,y = J.random_elements(2)
1950 sage: actual = (x*y).natural_representation()
1951 sage: X = x.natural_representation()
1952 sage: Y = y.natural_representation()
1953 sage: expected = (X*Y + Y*X)/2
1954 sage: actual == expected
1956 sage: J(expected) == x*y
1959 We can change the generator prefix::
1961 sage: QuaternionHermitianEJA(2, prefix='a').gens()
1962 (a0, a1, a2, a3, a4, a5)
1964 We can construct the (trivial) algebra of rank zero::
1966 sage: QuaternionHermitianEJA(0)
1967 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1971 def _denormalized_basis(cls
, n
, field
):
1973 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
1975 Why do we embed these? Basically, because all of numerical
1976 linear algebra assumes that you're working with vectors consisting
1977 of `n` entries from a field and scalars from the same field. There's
1978 no way to tell SageMath that (for example) the vectors contain
1979 complex numbers, while the scalar field is real.
1983 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1987 sage: set_random_seed()
1988 sage: n = ZZ.random_element(1,5)
1989 sage: B = QuaternionHermitianEJA._denormalized_basis(n,QQ)
1990 sage: all( M.is_symmetric() for M in B )
1994 Q
= QuaternionAlgebra(QQ
,-1,-1)
1997 # This is like the symmetric case, but we need to be careful:
1999 # * We want conjugate-symmetry, not just symmetry.
2000 # * The diagonal will (as a result) be real.
2004 for j
in range(i
+1):
2005 Eij
= matrix(Q
, n
, lambda k
,l
: k
==i
and l
==j
)
2007 Sij
= cls
.real_embed(Eij
)
2010 # The second, third, and fourth ones have a minus
2011 # because they're conjugated.
2012 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
2014 Sij_I
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
2016 Sij_J
= cls
.real_embed(J
*Eij
- J
*Eij
.transpose())
2018 Sij_K
= cls
.real_embed(K
*Eij
- K
*Eij
.transpose())
2021 # Since we embedded these, we can drop back to the "field" that we
2022 # started with instead of the quaternion algebra "Q".
2023 return ( s
.change_ring(field
) for s
in S
)
2026 def __init__(self
, n
, field
=AA
, **kwargs
):
2027 basis
= self
._denormalized
_basis
(n
,field
)
2028 super(QuaternionHermitianEJA
,self
).__init
__(field
,
2032 self
.rank
.set_cache(n
)
2035 def _max_random_instance_size():
2037 The maximum rank of a random QuaternionHermitianEJA.
2039 return 2 # Dimension 6
2042 def random_instance(cls
, field
=AA
, **kwargs
):
2044 Return a random instance of this type of algebra.
2046 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2047 return cls(n
, field
, **kwargs
)
2050 class HadamardEJA(RationalBasisEuclideanJordanAlgebra
,
2051 ConcreteEuclideanJordanAlgebra
):
2053 Return the Euclidean Jordan Algebra corresponding to the set
2054 `R^n` under the Hadamard product.
2056 Note: this is nothing more than the Cartesian product of ``n``
2057 copies of the spin algebra. Once Cartesian product algebras
2058 are implemented, this can go.
2062 sage: from mjo.eja.eja_algebra import HadamardEJA
2066 This multiplication table can be verified by hand::
2068 sage: J = HadamardEJA(3)
2069 sage: e0,e1,e2 = J.gens()
2085 We can change the generator prefix::
2087 sage: HadamardEJA(3, prefix='r').gens()
2091 def __init__(self
, n
, field
=AA
, **kwargs
):
2092 V
= VectorSpace(field
, n
)
2093 mult_table
= [ [ V
.gen(i
)*(i
== j
) for j
in range(n
) ]
2096 # Inner products are real numbers and not algebra
2097 # elements, so once we turn the algebra element
2098 # into a vector in inner_product(), we never go
2099 # back. As a result -- contrary to what we do with
2100 # self._multiplication_table -- we store the inner
2101 # product table as a plain old matrix and not as
2102 # an algebra operator.
2103 ip_table
= matrix
.identity(field
,n
)
2104 self
._inner
_product
_matrix
= ip_table
2106 super(HadamardEJA
, self
).__init
__(field
,
2110 self
.rank
.set_cache(n
)
2113 self
.one
.set_cache( self
.zero() )
2115 self
.one
.set_cache( sum(self
.gens()) )
2118 def _max_random_instance_size():
2120 The maximum dimension of a random HadamardEJA.
2125 def random_instance(cls
, field
=AA
, **kwargs
):
2127 Return a random instance of this type of algebra.
2129 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2130 return cls(n
, field
, **kwargs
)
2133 class BilinearFormEJA(RationalBasisEuclideanJordanAlgebra
,
2134 ConcreteEuclideanJordanAlgebra
):
2136 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2137 with the half-trace inner product and jordan product ``x*y =
2138 (<Bx,y>,y_bar>, x0*y_bar + y0*x_bar)`` where `B = 1 \times B22` is
2139 a symmetric positive-definite "bilinear form" matrix. Its
2140 dimension is the size of `B`, and it has rank two in dimensions
2141 larger than two. It reduces to the ``JordanSpinEJA`` when `B` is
2142 the identity matrix of order ``n``.
2144 We insist that the one-by-one upper-left identity block of `B` be
2145 passed in as well so that we can be passed a matrix of size zero
2146 to construct a trivial algebra.
2150 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
2151 ....: JordanSpinEJA)
2155 When no bilinear form is specified, the identity matrix is used,
2156 and the resulting algebra is the Jordan spin algebra::
2158 sage: B = matrix.identity(AA,3)
2159 sage: J0 = BilinearFormEJA(B)
2160 sage: J1 = JordanSpinEJA(3)
2161 sage: J0.multiplication_table() == J0.multiplication_table()
2164 An error is raised if the matrix `B` does not correspond to a
2165 positive-definite bilinear form::
2167 sage: B = matrix.random(QQ,2,3)
2168 sage: J = BilinearFormEJA(B)
2169 Traceback (most recent call last):
2171 ValueError: bilinear form is not positive-definite
2172 sage: B = matrix.zero(QQ,3)
2173 sage: J = BilinearFormEJA(B)
2174 Traceback (most recent call last):
2176 ValueError: bilinear form is not positive-definite
2180 We can create a zero-dimensional algebra::
2182 sage: B = matrix.identity(AA,0)
2183 sage: J = BilinearFormEJA(B)
2187 We can check the multiplication condition given in the Jordan, von
2188 Neumann, and Wigner paper (and also discussed on my "On the
2189 symmetry..." paper). Note that this relies heavily on the standard
2190 choice of basis, as does anything utilizing the bilinear form matrix::
2192 sage: set_random_seed()
2193 sage: n = ZZ.random_element(5)
2194 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2195 sage: B11 = matrix.identity(QQ,1)
2196 sage: B22 = M.transpose()*M
2197 sage: B = block_matrix(2,2,[ [B11,0 ],
2199 sage: J = BilinearFormEJA(B)
2200 sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
2201 sage: V = J.vector_space()
2202 sage: sis = [ J.from_vector(V([0] + (M.inverse()*ei).list()))
2203 ....: for ei in eis ]
2204 sage: actual = [ sis[i]*sis[j]
2205 ....: for i in range(n-1)
2206 ....: for j in range(n-1) ]
2207 sage: expected = [ J.one() if i == j else J.zero()
2208 ....: for i in range(n-1)
2209 ....: for j in range(n-1) ]
2210 sage: actual == expected
2213 def __init__(self
, B
, field
=AA
, **kwargs
):
2216 if not B
.is_positive_definite():
2217 raise ValueError("bilinear form is not positive-definite")
2219 V
= VectorSpace(field
, n
)
2220 mult_table
= [[V
.zero() for j
in range(n
)] for i
in range(n
)]
2229 z0
= (B
*x
).inner_product(y
)
2230 zbar
= y0
*xbar
+ x0
*ybar
2231 z
= V([z0
] + zbar
.list())
2232 mult_table
[i
][j
] = z
2234 # Inner products are real numbers and not algebra
2235 # elements, so once we turn the algebra element
2236 # into a vector in inner_product(), we never go
2237 # back. As a result -- contrary to what we do with
2238 # self._multiplication_table -- we store the inner
2239 # product table as a plain old matrix and not as
2240 # an algebra operator.
2242 self
._inner
_product
_matrix
= ip_table
2244 super(BilinearFormEJA
, self
).__init
__(field
,
2249 # The rank of this algebra is two, unless we're in a
2250 # one-dimensional ambient space (because the rank is bounded
2251 # by the ambient dimension).
2252 self
.rank
.set_cache(min(n
,2))
2255 self
.one
.set_cache( self
.zero() )
2257 self
.one
.set_cache( self
.monomial(0) )
2260 def _max_random_instance_size():
2262 The maximum dimension of a random BilinearFormEJA.
2267 def random_instance(cls
, field
=AA
, **kwargs
):
2269 Return a random instance of this algebra.
2271 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2273 B
= matrix
.identity(field
, n
)
2274 return cls(B
, field
, **kwargs
)
2276 B11
= matrix
.identity(field
,1)
2277 M
= matrix
.random(field
, n
-1)
2278 I
= matrix
.identity(field
, n
-1)
2279 alpha
= field
.zero()
2280 while alpha
.is_zero():
2281 alpha
= field
.random_element().abs()
2282 B22
= M
.transpose()*M
+ alpha
*I
2284 from sage
.matrix
.special
import block_matrix
2285 B
= block_matrix(2,2, [ [B11
, ZZ(0) ],
2288 return cls(B
, field
, **kwargs
)
2291 class JordanSpinEJA(BilinearFormEJA
):
2293 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2294 with the usual inner product and jordan product ``x*y =
2295 (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2300 sage: from mjo.eja.eja_algebra import JordanSpinEJA
2304 This multiplication table can be verified by hand::
2306 sage: J = JordanSpinEJA(4)
2307 sage: e0,e1,e2,e3 = J.gens()
2323 We can change the generator prefix::
2325 sage: JordanSpinEJA(2, prefix='B').gens()
2330 Ensure that we have the usual inner product on `R^n`::
2332 sage: set_random_seed()
2333 sage: J = JordanSpinEJA.random_instance()
2334 sage: x,y = J.random_elements(2)
2335 sage: actual = x.inner_product(y)
2336 sage: expected = x.to_vector().inner_product(y.to_vector())
2337 sage: actual == expected
2341 def __init__(self
, n
, field
=AA
, **kwargs
):
2342 # This is a special case of the BilinearFormEJA with the identity
2343 # matrix as its bilinear form.
2344 B
= matrix
.identity(field
, n
)
2345 super(JordanSpinEJA
, self
).__init
__(B
, field
, **kwargs
)
2348 def _max_random_instance_size():
2350 The maximum dimension of a random JordanSpinEJA.
2355 def random_instance(cls
, field
=AA
, **kwargs
):
2357 Return a random instance of this type of algebra.
2359 Needed here to override the implementation for ``BilinearFormEJA``.
2361 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2362 return cls(n
, field
, **kwargs
)
2365 class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra
,
2366 ConcreteEuclideanJordanAlgebra
):
2368 The trivial Euclidean Jordan algebra consisting of only a zero element.
2372 sage: from mjo.eja.eja_algebra import TrivialEJA
2376 sage: J = TrivialEJA()
2383 sage: 7*J.one()*12*J.one()
2385 sage: J.one().inner_product(J.one())
2387 sage: J.one().norm()
2389 sage: J.one().subalgebra_generated_by()
2390 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2395 def __init__(self
, field
=AA
, **kwargs
):
2397 self
._inner
_product
_matrix
= matrix(field
,0)
2398 super(TrivialEJA
, self
).__init
__(field
,
2402 # The rank is zero using my definition, namely the dimension of the
2403 # largest subalgebra generated by any element.
2404 self
.rank
.set_cache(0)
2405 self
.one
.set_cache( self
.zero() )
2408 def random_instance(cls
, field
=AA
, **kwargs
):
2409 # We don't take a "size" argument so the superclass method is
2410 # inappropriate for us.
2411 return cls(field
, **kwargs
)
2413 class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra
):
2415 The external (orthogonal) direct sum of two other Euclidean Jordan
2416 algebras. Essentially the Cartesian product of its two factors.
2417 Every Euclidean Jordan algebra decomposes into an orthogonal
2418 direct sum of simple Euclidean Jordan algebras, so no generality
2419 is lost by providing only this construction.
2423 sage: from mjo.eja.eja_algebra import (random_eja,
2425 ....: RealSymmetricEJA,
2430 sage: J1 = HadamardEJA(2)
2431 sage: J2 = RealSymmetricEJA(3)
2432 sage: J = DirectSumEJA(J1,J2)
2440 The external direct sum construction is only valid when the two factors
2441 have the same base ring; an error is raised otherwise::
2443 sage: set_random_seed()
2444 sage: J1 = random_eja(AA)
2445 sage: J2 = random_eja(QQ)
2446 sage: J = DirectSumEJA(J1,J2)
2447 Traceback (most recent call last):
2449 ValueError: algebras must share the same base field
2452 def __init__(self
, J1
, J2
, **kwargs
):
2453 if J1
.base_ring() != J2
.base_ring():
2454 raise ValueError("algebras must share the same base field")
2455 field
= J1
.base_ring()
2457 self
._factors
= (J1
, J2
)
2461 V
= VectorSpace(field
, n
)
2462 mult_table
= [ [ V
.zero() for j
in range(n
) ]
2466 p
= (J1
.monomial(i
)*J1
.monomial(j
)).to_vector()
2467 mult_table
[i
][j
] = V(p
.list() + [field
.zero()]*n2
)
2471 p
= (J2
.monomial(i
)*J2
.monomial(j
)).to_vector()
2472 mult_table
[n1
+i
][n1
+j
] = V([field
.zero()]*n1
+ p
.list())
2474 super(DirectSumEJA
, self
).__init
__(field
,
2478 self
.rank
.set_cache(J1
.rank() + J2
.rank())
2483 Return the pair of this algebra's factors.
2487 sage: from mjo.eja.eja_algebra import (HadamardEJA,
2488 ....: JordanSpinEJA,
2493 sage: J1 = HadamardEJA(2,QQ)
2494 sage: J2 = JordanSpinEJA(3,QQ)
2495 sage: J = DirectSumEJA(J1,J2)
2497 (Euclidean Jordan algebra of dimension 2 over Rational Field,
2498 Euclidean Jordan algebra of dimension 3 over Rational Field)
2501 return self
._factors
2503 def projections(self
):
2505 Return a pair of projections onto this algebra's factors.
2509 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
2510 ....: ComplexHermitianEJA,
2515 sage: J1 = JordanSpinEJA(2)
2516 sage: J2 = ComplexHermitianEJA(2)
2517 sage: J = DirectSumEJA(J1,J2)
2518 sage: (pi_left, pi_right) = J.projections()
2519 sage: J.one().to_vector()
2521 sage: pi_left(J.one()).to_vector()
2523 sage: pi_right(J.one()).to_vector()
2527 (J1
,J2
) = self
.factors()
2530 V_basis
= self
.vector_space().basis()
2531 # Need to specify the dimensions explicitly so that we don't
2532 # wind up with a zero-by-zero matrix when we want e.g. a
2533 # zero-by-two matrix (important for composing things).
2534 P1
= matrix(self
.base_ring(), m
, m
+n
, V_basis
[:m
])
2535 P2
= matrix(self
.base_ring(), n
, m
+n
, V_basis
[m
:])
2536 pi_left
= FiniteDimensionalEuclideanJordanAlgebraOperator(self
,J1
,P1
)
2537 pi_right
= FiniteDimensionalEuclideanJordanAlgebraOperator(self
,J2
,P2
)
2538 return (pi_left
, pi_right
)
2540 def inclusions(self
):
2542 Return the pair of inclusion maps from our factors into us.
2546 sage: from mjo.eja.eja_algebra import (random_eja,
2547 ....: JordanSpinEJA,
2548 ....: RealSymmetricEJA,
2553 sage: J1 = JordanSpinEJA(3)
2554 sage: J2 = RealSymmetricEJA(2)
2555 sage: J = DirectSumEJA(J1,J2)
2556 sage: (iota_left, iota_right) = J.inclusions()
2557 sage: iota_left(J1.zero()) == J.zero()
2559 sage: iota_right(J2.zero()) == J.zero()
2561 sage: J1.one().to_vector()
2563 sage: iota_left(J1.one()).to_vector()
2565 sage: J2.one().to_vector()
2567 sage: iota_right(J2.one()).to_vector()
2569 sage: J.one().to_vector()
2574 Composing a projection with the corresponding inclusion should
2575 produce the identity map, and mismatching them should produce
2578 sage: set_random_seed()
2579 sage: J1 = random_eja()
2580 sage: J2 = random_eja()
2581 sage: J = DirectSumEJA(J1,J2)
2582 sage: (iota_left, iota_right) = J.inclusions()
2583 sage: (pi_left, pi_right) = J.projections()
2584 sage: pi_left*iota_left == J1.one().operator()
2586 sage: pi_right*iota_right == J2.one().operator()
2588 sage: (pi_left*iota_right).is_zero()
2590 sage: (pi_right*iota_left).is_zero()
2594 (J1
,J2
) = self
.factors()
2597 V_basis
= self
.vector_space().basis()
2598 # Need to specify the dimensions explicitly so that we don't
2599 # wind up with a zero-by-zero matrix when we want e.g. a
2600 # two-by-zero matrix (important for composing things).
2601 I1
= matrix
.column(self
.base_ring(), m
, m
+n
, V_basis
[:m
])
2602 I2
= matrix
.column(self
.base_ring(), n
, m
+n
, V_basis
[m
:])
2603 iota_left
= FiniteDimensionalEuclideanJordanAlgebraOperator(J1
,self
,I1
)
2604 iota_right
= FiniteDimensionalEuclideanJordanAlgebraOperator(J2
,self
,I2
)
2605 return (iota_left
, iota_right
)
2607 def inner_product(self
, x
, y
):
2609 The standard Cartesian inner-product.
2611 We project ``x`` and ``y`` onto our factors, and add up the
2612 inner-products from the subalgebras.
2617 sage: from mjo.eja.eja_algebra import (HadamardEJA,
2618 ....: QuaternionHermitianEJA,
2623 sage: J1 = HadamardEJA(3,QQ)
2624 sage: J2 = QuaternionHermitianEJA(2,QQ,normalize_basis=False)
2625 sage: J = DirectSumEJA(J1,J2)
2630 sage: x1.inner_product(x2)
2632 sage: y1.inner_product(y2)
2634 sage: J.one().inner_product(J.one())
2638 (pi_left
, pi_right
) = self
.projections()
2644 return (x1
.inner_product(y1
) + x2
.inner_product(y2
))
2648 random_eja
= ConcreteEuclideanJordanAlgebra
.random_instance