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 a naturally-represented algebra element
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 ArithmeticError: vector is not in free module
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 a naturally-represented algebra element"
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
)
224 coords
= W
.coordinate_vector(_mat2vec(elt
))
225 return self
.from_vector(coords
)
229 Return a string representation of ``self``.
233 sage: from mjo.eja.eja_algebra import JordanSpinEJA
237 Ensure that it says what we think it says::
239 sage: JordanSpinEJA(2, field=AA)
240 Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
241 sage: JordanSpinEJA(3, field=RDF)
242 Euclidean Jordan algebra of dimension 3 over Real Double Field
245 fmt
= "Euclidean Jordan algebra of dimension {} over {}"
246 return fmt
.format(self
.dimension(), self
.base_ring())
248 def product_on_basis(self
, i
, j
):
249 return self
._multiplication
_table
[i
][j
]
251 def _is_commutative(self
):
253 Whether or not this algebra's multiplication table is commutative.
255 This method should of course always return ``True``, unless
256 this algebra was constructed with ``check_axioms=False`` and
257 passed an invalid multiplication table.
259 return all( self
.product_on_basis(i
,j
) == self
.product_on_basis(i
,j
)
260 for i
in range(self
.dimension())
261 for j
in range(self
.dimension()) )
263 def _is_jordanian(self
):
265 Whether or not this algebra's multiplication table respects the
266 Jordan identity `(x^{2})(xy) = x(x^{2}y)`.
268 We only check one arrangement of `x` and `y`, so for a
269 ``True`` result to be truly true, you should also check
270 :meth:`_is_commutative`. This method should of course always
271 return ``True``, unless this algebra was constructed with
272 ``check_axioms=False`` and passed an invalid multiplication table.
274 return all( (self
.monomial(i
)**2)*(self
.monomial(i
)*self
.monomial(j
))
276 (self
.monomial(i
))*((self
.monomial(i
)**2)*self
.monomial(j
))
277 for i
in range(self
.dimension())
278 for j
in range(self
.dimension()) )
280 def _inner_product_is_associative(self
):
282 Return whether or not this algebra's inner product `B` is
283 associative; that is, whether or not `B(xy,z) = B(x,yz)`.
285 This method should of course always return ``True``, unless
286 this algebra was constructed with ``check_axioms=False`` and
287 passed an invalid multiplication table.
290 # Used to check whether or not something is zero in an inexact
291 # ring. This number is sufficient to allow the construction of
292 # QuaternionHermitianEJA(2, RDF) with check_axioms=True.
295 for i
in range(self
.dimension()):
296 for j
in range(self
.dimension()):
297 for k
in range(self
.dimension()):
301 diff
= (x
*y
).inner_product(z
) - x
.inner_product(y
*z
)
303 if self
.base_ring().is_exact():
307 if diff
.abs() > epsilon
:
313 def characteristic_polynomial_of(self
):
315 Return the algebra's "characteristic polynomial of" function,
316 which is itself a multivariate polynomial that, when evaluated
317 at the coordinates of some algebra element, returns that
318 element's characteristic polynomial.
320 The resulting polynomial has `n+1` variables, where `n` is the
321 dimension of this algebra. The first `n` variables correspond to
322 the coordinates of an algebra element: when evaluated at the
323 coordinates of an algebra element with respect to a certain
324 basis, the result is a univariate polynomial (in the one
325 remaining variable ``t``), namely the characteristic polynomial
330 sage: from mjo.eja.eja_algebra import JordanSpinEJA, TrivialEJA
334 The characteristic polynomial in the spin algebra is given in
335 Alizadeh, Example 11.11::
337 sage: J = JordanSpinEJA(3)
338 sage: p = J.characteristic_polynomial_of(); p
339 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
340 sage: xvec = J.one().to_vector()
344 By definition, the characteristic polynomial is a monic
345 degree-zero polynomial in a rank-zero algebra. Note that
346 Cayley-Hamilton is indeed satisfied since the polynomial
347 ``1`` evaluates to the identity element of the algebra on
350 sage: J = TrivialEJA()
351 sage: J.characteristic_polynomial_of()
358 # The list of coefficient polynomials a_0, a_1, a_2, ..., a_(r-1).
359 a
= self
._charpoly
_coefficients
()
361 # We go to a bit of trouble here to reorder the
362 # indeterminates, so that it's easier to evaluate the
363 # characteristic polynomial at x's coordinates and get back
364 # something in terms of t, which is what we want.
365 S
= PolynomialRing(self
.base_ring(),'t')
369 S
= PolynomialRing(S
, R
.variable_names())
372 return (t
**r
+ sum( a
[k
]*(t
**k
) for k
in range(r
) ))
374 def coordinate_polynomial_ring(self
):
376 The multivariate polynomial ring in which this algebra's
377 :meth:`characteristic_polynomial_of` lives.
381 sage: from mjo.eja.eja_algebra import (HadamardEJA,
382 ....: RealSymmetricEJA)
386 sage: J = HadamardEJA(2)
387 sage: J.coordinate_polynomial_ring()
388 Multivariate Polynomial Ring in X1, X2...
389 sage: J = RealSymmetricEJA(3,QQ)
390 sage: J.coordinate_polynomial_ring()
391 Multivariate Polynomial Ring in X1, X2, X3, X4, X5, X6...
394 var_names
= tuple( "X%d" % z
for z
in range(1, self
.dimension()+1) )
395 return PolynomialRing(self
.base_ring(), var_names
)
397 def inner_product(self
, x
, y
):
399 The inner product associated with this Euclidean Jordan algebra.
401 Defaults to the trace inner product, but can be overridden by
402 subclasses if they are sure that the necessary properties are
407 sage: from mjo.eja.eja_algebra import (random_eja,
409 ....: BilinearFormEJA)
413 Our inner product is "associative," which means the following for
414 a symmetric bilinear form::
416 sage: set_random_seed()
417 sage: J = random_eja()
418 sage: x,y,z = J.random_elements(3)
419 sage: (x*y).inner_product(z) == y.inner_product(x*z)
424 Ensure that this is the usual inner product for the algebras
427 sage: set_random_seed()
428 sage: J = HadamardEJA.random_instance()
429 sage: x,y = J.random_elements(2)
430 sage: actual = x.inner_product(y)
431 sage: expected = x.to_vector().inner_product(y.to_vector())
432 sage: actual == expected
435 Ensure that this is one-half of the trace inner-product in a
436 BilinearFormEJA that isn't just the reals (when ``n`` isn't
437 one). This is in Faraut and Koranyi, and also my "On the
440 sage: set_random_seed()
441 sage: J = BilinearFormEJA.random_instance()
442 sage: n = J.dimension()
443 sage: x = J.random_element()
444 sage: y = J.random_element()
445 sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
448 B
= self
._inner
_product
_matrix
449 return (B
*x
.to_vector()).inner_product(y
.to_vector())
452 def is_trivial(self
):
454 Return whether or not this algebra is trivial.
456 A trivial algebra contains only the zero element.
460 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
465 sage: J = ComplexHermitianEJA(3)
471 sage: J = TrivialEJA()
476 return self
.dimension() == 0
479 def multiplication_table(self
):
481 Return a visual representation of this algebra's multiplication
482 table (on basis elements).
486 sage: from mjo.eja.eja_algebra import JordanSpinEJA
490 sage: J = JordanSpinEJA(4)
491 sage: J.multiplication_table()
492 +----++----+----+----+----+
493 | * || e0 | e1 | e2 | e3 |
494 +====++====+====+====+====+
495 | e0 || e0 | e1 | e2 | e3 |
496 +----++----+----+----+----+
497 | e1 || e1 | e0 | 0 | 0 |
498 +----++----+----+----+----+
499 | e2 || e2 | 0 | e0 | 0 |
500 +----++----+----+----+----+
501 | e3 || e3 | 0 | 0 | e0 |
502 +----++----+----+----+----+
505 M
= list(self
._multiplication
_table
) # copy
506 for i
in range(len(M
)):
507 # M had better be "square"
508 M
[i
] = [self
.monomial(i
)] + M
[i
]
509 M
= [["*"] + list(self
.gens())] + M
510 return table(M
, header_row
=True, header_column
=True, frame
=True)
513 def natural_basis(self
):
515 Return a more-natural representation of this algebra's basis.
517 Every finite-dimensional Euclidean Jordan Algebra is a direct
518 sum of five simple algebras, four of which comprise Hermitian
519 matrices. This method returns the original "natural" basis
520 for our underlying vector space. (Typically, the natural basis
521 is used to construct the multiplication table in the first place.)
523 Note that this will always return a matrix. The standard basis
524 in `R^n` will be returned as `n`-by-`1` column matrices.
528 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
529 ....: RealSymmetricEJA)
533 sage: J = RealSymmetricEJA(2)
535 Finite family {0: e0, 1: e1, 2: e2}
536 sage: J.natural_basis()
538 [1 0] [ 0 0.7071067811865475?] [0 0]
539 [0 0], [0.7071067811865475? 0], [0 1]
544 sage: J = JordanSpinEJA(2)
546 Finite family {0: e0, 1: e1}
547 sage: J.natural_basis()
554 if self
._natural
_basis
is None:
555 M
= self
.natural_basis_space()
556 return tuple( M(b
.to_vector()) for b
in self
.basis() )
558 return self
._natural
_basis
561 def natural_basis_space(self
):
563 Return the matrix space in which this algebra's natural basis
566 Generally this will be an `n`-by-`1` column-vector space,
567 except when the algebra is trivial. There it's `n`-by-`n`
568 (where `n` is zero), to ensure that two elements of the
569 natural basis space (empty matrices) can be multiplied.
571 if self
.is_trivial():
572 return MatrixSpace(self
.base_ring(), 0)
573 elif self
._natural
_basis
is None or len(self
._natural
_basis
) == 0:
574 return MatrixSpace(self
.base_ring(), self
.dimension(), 1)
576 return self
._natural
_basis
[0].matrix_space()
582 Return the unit element of this algebra.
586 sage: from mjo.eja.eja_algebra import (HadamardEJA,
591 sage: J = HadamardEJA(5)
593 e0 + e1 + e2 + e3 + e4
597 The identity element acts like the identity::
599 sage: set_random_seed()
600 sage: J = random_eja()
601 sage: x = J.random_element()
602 sage: J.one()*x == x and x*J.one() == x
605 The matrix of the unit element's operator is the identity::
607 sage: set_random_seed()
608 sage: J = random_eja()
609 sage: actual = J.one().operator().matrix()
610 sage: expected = matrix.identity(J.base_ring(), J.dimension())
611 sage: actual == expected
614 Ensure that the cached unit element (often precomputed by
615 hand) agrees with the computed one::
617 sage: set_random_seed()
618 sage: J = random_eja()
619 sage: cached = J.one()
620 sage: J.one.clear_cache()
621 sage: J.one() == cached
625 # We can brute-force compute the matrices of the operators
626 # that correspond to the basis elements of this algebra.
627 # If some linear combination of those basis elements is the
628 # algebra identity, then the same linear combination of
629 # their matrices has to be the identity matrix.
631 # Of course, matrices aren't vectors in sage, so we have to
632 # appeal to the "long vectors" isometry.
633 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
635 # Now we use basic linear algebra to find the coefficients,
636 # of the matrices-as-vectors-linear-combination, which should
637 # work for the original algebra basis too.
638 A
= matrix(self
.base_ring(), oper_vecs
)
640 # We used the isometry on the left-hand side already, but we
641 # still need to do it for the right-hand side. Recall that we
642 # wanted something that summed to the identity matrix.
643 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
645 # Now if there's an identity element in the algebra, this
646 # should work. We solve on the left to avoid having to
647 # transpose the matrix "A".
648 return self
.from_vector(A
.solve_left(b
))
651 def peirce_decomposition(self
, c
):
653 The Peirce decomposition of this algebra relative to the
656 In the future, this can be extended to a complete system of
657 orthogonal idempotents.
661 - ``c`` -- an idempotent of this algebra.
665 A triple (J0, J5, J1) containing two subalgebras and one subspace
668 - ``J0`` -- the algebra on the eigenspace of ``c.operator()``
669 corresponding to the eigenvalue zero.
671 - ``J5`` -- the eigenspace (NOT a subalgebra) of ``c.operator()``
672 corresponding to the eigenvalue one-half.
674 - ``J1`` -- the algebra on the eigenspace of ``c.operator()``
675 corresponding to the eigenvalue one.
677 These are the only possible eigenspaces for that operator, and this
678 algebra is a direct sum of them. The spaces ``J0`` and ``J1`` are
679 orthogonal, and are subalgebras of this algebra with the appropriate
684 sage: from mjo.eja.eja_algebra import random_eja, RealSymmetricEJA
688 The canonical example comes from the symmetric matrices, which
689 decompose into diagonal and off-diagonal parts::
691 sage: J = RealSymmetricEJA(3)
692 sage: C = matrix(QQ, [ [1,0,0],
696 sage: J0,J5,J1 = J.peirce_decomposition(c)
698 Euclidean Jordan algebra of dimension 1...
700 Vector space of degree 6 and dimension 2...
702 Euclidean Jordan algebra of dimension 3...
703 sage: J0.one().natural_representation()
707 sage: orig_df = AA.options.display_format
708 sage: AA.options.display_format = 'radical'
709 sage: J.from_vector(J5.basis()[0]).natural_representation()
713 sage: J.from_vector(J5.basis()[1]).natural_representation()
717 sage: AA.options.display_format = orig_df
718 sage: J1.one().natural_representation()
725 Every algebra decomposes trivially with respect to its identity
728 sage: set_random_seed()
729 sage: J = random_eja()
730 sage: J0,J5,J1 = J.peirce_decomposition(J.one())
731 sage: J0.dimension() == 0 and J5.dimension() == 0
733 sage: J1.superalgebra() == J and J1.dimension() == J.dimension()
736 The decomposition is into eigenspaces, and its components are
737 therefore necessarily orthogonal. Moreover, the identity
738 elements in the two subalgebras are the projections onto their
739 respective subspaces of the superalgebra's identity element::
741 sage: set_random_seed()
742 sage: J = random_eja()
743 sage: x = J.random_element()
744 sage: if not J.is_trivial():
745 ....: while x.is_nilpotent():
746 ....: x = J.random_element()
747 sage: c = x.subalgebra_idempotent()
748 sage: J0,J5,J1 = J.peirce_decomposition(c)
750 sage: for (w,y,z) in zip(J0.basis(), J5.basis(), J1.basis()):
751 ....: w = w.superalgebra_element()
752 ....: y = J.from_vector(y)
753 ....: z = z.superalgebra_element()
754 ....: ipsum += w.inner_product(y).abs()
755 ....: ipsum += w.inner_product(z).abs()
756 ....: ipsum += y.inner_product(z).abs()
759 sage: J1(c) == J1.one()
761 sage: J0(J.one() - c) == J0.one()
765 if not c
.is_idempotent():
766 raise ValueError("element is not idempotent: %s" % c
)
768 from mjo
.eja
.eja_subalgebra
import FiniteDimensionalEuclideanJordanSubalgebra
770 # Default these to what they should be if they turn out to be
771 # trivial, because eigenspaces_left() won't return eigenvalues
772 # corresponding to trivial spaces (e.g. it returns only the
773 # eigenspace corresponding to lambda=1 if you take the
774 # decomposition relative to the identity element).
775 trivial
= FiniteDimensionalEuclideanJordanSubalgebra(self
, ())
776 J0
= trivial
# eigenvalue zero
777 J5
= VectorSpace(self
.base_ring(), 0) # eigenvalue one-half
778 J1
= trivial
# eigenvalue one
780 for (eigval
, eigspace
) in c
.operator().matrix().right_eigenspaces():
781 if eigval
== ~
(self
.base_ring()(2)):
784 gens
= tuple( self
.from_vector(b
) for b
in eigspace
.basis() )
785 subalg
= FiniteDimensionalEuclideanJordanSubalgebra(self
,
793 raise ValueError("unexpected eigenvalue: %s" % eigval
)
798 def random_element(self
, thorough
=False):
800 Return a random element of this algebra.
802 Our algebra superclass method only returns a linear
803 combination of at most two basis elements. We instead
804 want the vector space "random element" method that
805 returns a more diverse selection.
809 - ``thorough`` -- (boolean; default False) whether or not we
810 should generate irrational coefficients for the random
811 element when our base ring is irrational; this slows the
812 algebra operations to a crawl, but any truly random method
816 # For a general base ring... maybe we can trust this to do the
817 # right thing? Unlikely, but.
818 V
= self
.vector_space()
819 v
= V
.random_element()
821 if self
.base_ring() is AA
:
822 # The "random element" method of the algebraic reals is
823 # stupid at the moment, and only returns integers between
824 # -2 and 2, inclusive:
826 # https://trac.sagemath.org/ticket/30875
828 # Instead, we implement our own "random vector" method,
829 # and then coerce that into the algebra. We use the vector
830 # space degree here instead of the dimension because a
831 # subalgebra could (for example) be spanned by only two
832 # vectors, each with five coordinates. We need to
833 # generate all five coordinates.
835 v
*= QQbar
.random_element().real()
837 v
*= QQ
.random_element()
839 return self
.from_vector(V
.coordinate_vector(v
))
841 def random_elements(self
, count
, thorough
=False):
843 Return ``count`` random elements as a tuple.
847 - ``thorough`` -- (boolean; default False) whether or not we
848 should generate irrational coefficients for the random
849 elements when our base ring is irrational; this slows the
850 algebra operations to a crawl, but any truly random method
855 sage: from mjo.eja.eja_algebra import JordanSpinEJA
859 sage: J = JordanSpinEJA(3)
860 sage: x,y,z = J.random_elements(3)
861 sage: all( [ x in J, y in J, z in J ])
863 sage: len( J.random_elements(10) ) == 10
867 return tuple( self
.random_element(thorough
)
868 for idx
in range(count
) )
872 def _charpoly_coefficients(self
):
874 The `r` polynomial coefficients of the "characteristic polynomial
878 R
= self
.coordinate_polynomial_ring()
880 F
= R
.fraction_field()
883 # From a result in my book, these are the entries of the
884 # basis representation of L_x.
885 return sum( vars[k
]*self
.monomial(k
).operator().matrix()[i
,j
]
888 L_x
= matrix(F
, n
, n
, L_x_i_j
)
891 if self
.rank
.is_in_cache():
893 # There's no need to pad the system with redundant
894 # columns if we *know* they'll be redundant.
897 # Compute an extra power in case the rank is equal to
898 # the dimension (otherwise, we would stop at x^(r-1)).
899 x_powers
= [ (L_x
**k
)*self
.one().to_vector()
900 for k
in range(n
+1) ]
901 A
= matrix
.column(F
, x_powers
[:n
])
902 AE
= A
.extended_echelon_form()
909 # The theory says that only the first "r" coefficients are
910 # nonzero, and they actually live in the original polynomial
911 # ring and not the fraction field. We negate them because
912 # in the actual characteristic polynomial, they get moved
913 # to the other side where x^r lives.
914 return -A_rref
.solve_right(E
*b
).change_ring(R
)[:r
]
919 Return the rank of this EJA.
921 This is a cached method because we know the rank a priori for
922 all of the algebras we can construct. Thus we can avoid the
923 expensive ``_charpoly_coefficients()`` call unless we truly
924 need to compute the whole characteristic polynomial.
928 sage: from mjo.eja.eja_algebra import (HadamardEJA,
930 ....: RealSymmetricEJA,
931 ....: ComplexHermitianEJA,
932 ....: QuaternionHermitianEJA,
937 The rank of the Jordan spin algebra is always two::
939 sage: JordanSpinEJA(2).rank()
941 sage: JordanSpinEJA(3).rank()
943 sage: JordanSpinEJA(4).rank()
946 The rank of the `n`-by-`n` Hermitian real, complex, or
947 quaternion matrices is `n`::
949 sage: RealSymmetricEJA(4).rank()
951 sage: ComplexHermitianEJA(3).rank()
953 sage: QuaternionHermitianEJA(2).rank()
958 Ensure that every EJA that we know how to construct has a
959 positive integer rank, unless the algebra is trivial in
960 which case its rank will be zero::
962 sage: set_random_seed()
963 sage: J = random_eja()
967 sage: r > 0 or (r == 0 and J.is_trivial())
970 Ensure that computing the rank actually works, since the ranks
971 of all simple algebras are known and will be cached by default::
973 sage: J = HadamardEJA(4)
974 sage: J.rank.clear_cache()
980 sage: J = JordanSpinEJA(4)
981 sage: J.rank.clear_cache()
987 sage: J = RealSymmetricEJA(3)
988 sage: J.rank.clear_cache()
994 sage: J = ComplexHermitianEJA(2)
995 sage: J.rank.clear_cache()
1001 sage: J = QuaternionHermitianEJA(2)
1002 sage: J.rank.clear_cache()
1006 return len(self
._charpoly
_coefficients
())
1009 def vector_space(self
):
1011 Return the vector space that underlies this algebra.
1015 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1019 sage: J = RealSymmetricEJA(2)
1020 sage: J.vector_space()
1021 Vector space of dimension 3 over...
1024 return self
.zero().to_vector().parent().ambient_vector_space()
1027 Element
= FiniteDimensionalEuclideanJordanAlgebraElement
1030 class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra
):
1032 Algebras whose basis consists of vectors with rational
1033 entries. Equivalently, algebras whose multiplication tables
1034 contain only rational coefficients.
1036 When an EJA has a basis that can be made rational, we can speed up
1037 the computation of its characteristic polynomial by doing it over
1038 ``QQ``. All of the named EJA constructors that we provide fall
1042 def _charpoly_coefficients(self
):
1044 Override the parent method with something that tries to compute
1045 over a faster (non-extension) field.
1049 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1053 The base ring of the resulting polynomial coefficients is what
1054 it should be, and not the rationals (unless the algebra was
1055 already over the rationals)::
1057 sage: J = JordanSpinEJA(3)
1058 sage: J._charpoly_coefficients()
1059 (X1^2 - X2^2 - X3^2, -2*X1)
1060 sage: a0 = J._charpoly_coefficients()[0]
1062 Algebraic Real Field
1063 sage: a0.base_ring()
1064 Algebraic Real Field
1067 if self
.base_ring() is QQ
:
1068 # There's no need to construct *another* algebra over the
1069 # rationals if this one is already over the rationals.
1070 superclass
= super(RationalBasisEuclideanJordanAlgebra
, self
)
1071 return superclass
._charpoly
_coefficients
()
1074 map(lambda x
: x
.to_vector(), ls
)
1075 for ls
in self
._multiplication
_table
1078 # Do the computation over the rationals. The answer will be
1079 # the same, because our basis coordinates are (essentially)
1081 J
= FiniteDimensionalEuclideanJordanAlgebra(QQ
,
1085 a
= J
._charpoly
_coefficients
()
1086 return tuple(map(lambda x
: x
.change_ring(self
.base_ring()), a
))
1089 class ConcreteEuclideanJordanAlgebra
:
1091 A class for the Euclidean Jordan algebras that we know by name.
1093 These are the Jordan algebras whose basis, multiplication table,
1094 rank, and so on are known a priori. More to the point, they are
1095 the Euclidean Jordan algebras for which we are able to conjure up
1096 a "random instance."
1100 sage: from mjo.eja.eja_algebra import ConcreteEuclideanJordanAlgebra
1104 Our natural basis is normalized with respect to the natural inner
1105 product unless we specify otherwise::
1107 sage: set_random_seed()
1108 sage: J = ConcreteEuclideanJordanAlgebra.random_instance()
1109 sage: all( b.norm() == 1 for b in J.gens() )
1112 Since our natural basis is normalized with respect to the natural
1113 inner product, and since we know that this algebra is an EJA, any
1114 left-multiplication operator's matrix will be symmetric because
1115 natural->EJA basis representation is an isometry and within the EJA
1116 the operator is self-adjoint by the Jordan axiom::
1118 sage: set_random_seed()
1119 sage: J = ConcreteEuclideanJordanAlgebra.random_instance()
1120 sage: x = J.random_element()
1121 sage: x.operator().matrix().is_symmetric()
1127 def _max_random_instance_size():
1129 Return an integer "size" that is an upper bound on the size of
1130 this algebra when it is used in a random test
1131 case. Unfortunately, the term "size" is ambiguous -- when
1132 dealing with `R^n` under either the Hadamard or Jordan spin
1133 product, the "size" refers to the dimension `n`. When dealing
1134 with a matrix algebra (real symmetric or complex/quaternion
1135 Hermitian), it refers to the size of the matrix, which is far
1136 less than the dimension of the underlying vector space.
1138 This method must be implemented in each subclass.
1140 raise NotImplementedError
1143 def random_instance(cls
, field
=AA
, **kwargs
):
1145 Return a random instance of this type of algebra.
1147 This method should be implemented in each subclass.
1149 from sage
.misc
.prandom
import choice
1150 eja_class
= choice(cls
.__subclasses
__())
1151 return eja_class
.random_instance(field
)
1154 class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra
):
1156 def __init__(self
, field
, basis
, normalize_basis
=True, **kwargs
):
1158 Compared to the superclass constructor, we take a basis instead of
1159 a multiplication table because the latter can be computed in terms
1160 of the former when the product is known (like it is here).
1162 # Used in this class's fast _charpoly_coefficients() override.
1163 self
._basis
_normalizers
= None
1165 # We're going to loop through this a few times, so now's a good
1166 # time to ensure that it isn't a generator expression.
1167 basis
= tuple(basis
)
1169 algebra_dim
= len(basis
)
1170 degree
= 0 # size of the matrices
1172 degree
= basis
[0].nrows()
1174 if algebra_dim
> 1 and normalize_basis
:
1175 # We'll need sqrt(2) to normalize the basis, and this
1176 # winds up in the multiplication table, so the whole
1177 # algebra needs to be over the field extension.
1178 R
= PolynomialRing(field
, 'z')
1181 if p
.is_irreducible():
1182 field
= field
.extension(p
, 'sqrt2', embedding
=RLF(2).sqrt())
1183 basis
= tuple( s
.change_ring(field
) for s
in basis
)
1184 self
._basis
_normalizers
= tuple(
1185 ~
(self
.natural_inner_product(s
,s
).sqrt()) for s
in basis
)
1186 basis
= tuple(s
*c
for (s
,c
) in zip(basis
,self
._basis
_normalizers
))
1188 # Now compute the multiplication and inner product tables.
1189 # We have to do this *after* normalizing the basis, because
1190 # scaling affects the answers.
1191 V
= VectorSpace(field
, degree
**2)
1192 W
= V
.span_of_basis( _mat2vec(s
) for s
in basis
)
1193 mult_table
= [[W
.zero() for j
in range(algebra_dim
)]
1194 for i
in range(algebra_dim
)]
1195 ip_table
= [[W
.zero() for j
in range(algebra_dim
)]
1196 for i
in range(algebra_dim
)]
1197 for i
in range(algebra_dim
):
1198 for j
in range(algebra_dim
):
1199 mat_entry
= (basis
[i
]*basis
[j
] + basis
[j
]*basis
[i
])/2
1200 mult_table
[i
][j
] = W
.coordinate_vector(_mat2vec(mat_entry
))
1203 # HACK: ignore the error here if we don't need the
1204 # inner product (as is the case when we construct
1205 # a dummy QQ-algebra for fast charpoly coefficients.
1206 ip_table
[i
][j
] = self
.natural_inner_product(basis
[i
],
1213 self
._inner
_product
_matrix
= matrix(field
,ip_table
)
1217 super(MatrixEuclideanJordanAlgebra
, self
).__init
__(field
,
1219 natural_basis
=basis
,
1222 if algebra_dim
== 0:
1223 self
.one
.set_cache(self
.zero())
1225 n
= basis
[0].nrows()
1226 # The identity wrt (A,B) -> (AB + BA)/2 is independent of the
1227 # details of this algebra.
1228 self
.one
.set_cache(self(matrix
.identity(field
,n
)))
1232 def _charpoly_coefficients(self
):
1234 Override the parent method with something that tries to compute
1235 over a faster (non-extension) field.
1237 if self
._basis
_normalizers
is None or self
.base_ring() is QQ
:
1238 # We didn't normalize, or the basis we started with had
1239 # entries in a nice field already. Just compute the thing.
1240 return super(MatrixEuclideanJordanAlgebra
, self
)._charpoly
_coefficients
()
1242 basis
= ( (b
/n
) for (b
,n
) in zip(self
.natural_basis(),
1243 self
._basis
_normalizers
) )
1245 # Do this over the rationals and convert back at the end.
1246 # Only works because we know the entries of the basis are
1247 # integers. The argument ``check_axioms=False`` is required
1248 # because the trace inner-product method for this
1249 # class is a stub and can't actually be checked.
1250 J
= MatrixEuclideanJordanAlgebra(QQ
,
1252 normalize_basis
=False,
1255 a
= J
._charpoly
_coefficients
()
1257 # Unfortunately, changing the basis does change the
1258 # coefficients of the characteristic polynomial, but since
1259 # these are really the coefficients of the "characteristic
1260 # polynomial of" function, everything is still nice and
1261 # unevaluated. It's therefore "obvious" how scaling the
1262 # basis affects the coordinate variables X1, X2, et
1263 # cetera. Scaling the first basis vector up by "n" adds a
1264 # factor of 1/n into every "X1" term, for example. So here
1265 # we simply undo the basis_normalizer scaling that we
1266 # performed earlier.
1268 # The a[0] access here is safe because trivial algebras
1269 # won't have any basis normalizers and therefore won't
1270 # make it to this "else" branch.
1271 XS
= a
[0].parent().gens()
1272 subs_dict
= { XS
[i
]: self
._basis
_normalizers
[i
]*XS
[i
]
1273 for i
in range(len(XS
)) }
1274 return tuple( a_i
.subs(subs_dict
) for a_i
in a
)
1280 Embed the matrix ``M`` into a space of real matrices.
1282 The matrix ``M`` can have entries in any field at the moment:
1283 the real numbers, complex numbers, or quaternions. And although
1284 they are not a field, we can probably support octonions at some
1285 point, too. This function returns a real matrix that "acts like"
1286 the original with respect to matrix multiplication; i.e.
1288 real_embed(M*N) = real_embed(M)*real_embed(N)
1291 raise NotImplementedError
1295 def real_unembed(M
):
1297 The inverse of :meth:`real_embed`.
1299 raise NotImplementedError
1302 def natural_inner_product(cls
,X
,Y
):
1303 Xu
= cls
.real_unembed(X
)
1304 Yu
= cls
.real_unembed(Y
)
1305 tr
= (Xu
*Yu
).trace()
1308 # Works in QQ, AA, RDF, et cetera.
1310 except AttributeError:
1311 # A quaternion doesn't have a real() method, but does
1312 # have coefficient_tuple() method that returns the
1313 # coefficients of 1, i, j, and k -- in that order.
1314 return tr
.coefficient_tuple()[0]
1317 class RealMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1321 The identity function, for embedding real matrices into real
1327 def real_unembed(M
):
1329 The identity function, for unembedding real matrices from real
1335 class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra
,
1336 ConcreteEuclideanJordanAlgebra
):
1338 The rank-n simple EJA consisting of real symmetric n-by-n
1339 matrices, the usual symmetric Jordan product, and the trace inner
1340 product. It has dimension `(n^2 + n)/2` over the reals.
1344 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1348 sage: J = RealSymmetricEJA(2)
1349 sage: e0, e1, e2 = J.gens()
1357 In theory, our "field" can be any subfield of the reals::
1359 sage: RealSymmetricEJA(2, RDF)
1360 Euclidean Jordan algebra of dimension 3 over Real Double Field
1361 sage: RealSymmetricEJA(2, RR)
1362 Euclidean Jordan algebra of dimension 3 over Real Field with
1363 53 bits of precision
1367 The dimension of this algebra is `(n^2 + n) / 2`::
1369 sage: set_random_seed()
1370 sage: n_max = RealSymmetricEJA._max_random_instance_size()
1371 sage: n = ZZ.random_element(1, n_max)
1372 sage: J = RealSymmetricEJA(n)
1373 sage: J.dimension() == (n^2 + n)/2
1376 The Jordan multiplication is what we think it is::
1378 sage: set_random_seed()
1379 sage: J = RealSymmetricEJA.random_instance()
1380 sage: x,y = J.random_elements(2)
1381 sage: actual = (x*y).natural_representation()
1382 sage: X = x.natural_representation()
1383 sage: Y = y.natural_representation()
1384 sage: expected = (X*Y + Y*X)/2
1385 sage: actual == expected
1387 sage: J(expected) == x*y
1390 We can change the generator prefix::
1392 sage: RealSymmetricEJA(3, prefix='q').gens()
1393 (q0, q1, q2, q3, q4, q5)
1395 We can construct the (trivial) algebra of rank zero::
1397 sage: RealSymmetricEJA(0)
1398 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1402 def _denormalized_basis(cls
, n
, field
):
1404 Return a basis for the space of real symmetric n-by-n matrices.
1408 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1412 sage: set_random_seed()
1413 sage: n = ZZ.random_element(1,5)
1414 sage: B = RealSymmetricEJA._denormalized_basis(n,QQ)
1415 sage: all( M.is_symmetric() for M in B)
1419 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1423 for j
in range(i
+1):
1424 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
1428 Sij
= Eij
+ Eij
.transpose()
1434 def _max_random_instance_size():
1435 return 4 # Dimension 10
1438 def random_instance(cls
, field
=AA
, **kwargs
):
1440 Return a random instance of this type of algebra.
1442 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1443 return cls(n
, field
, **kwargs
)
1445 def __init__(self
, n
, field
=AA
, **kwargs
):
1446 basis
= self
._denormalized
_basis
(n
, field
)
1447 super(RealSymmetricEJA
, self
).__init
__(field
,
1451 self
.rank
.set_cache(n
)
1454 class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1458 Embed the n-by-n complex matrix ``M`` into the space of real
1459 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1460 bi` to the block matrix ``[[a,b],[-b,a]]``.
1464 sage: from mjo.eja.eja_algebra import \
1465 ....: ComplexMatrixEuclideanJordanAlgebra
1469 sage: F = QuadraticField(-1, 'I')
1470 sage: x1 = F(4 - 2*i)
1471 sage: x2 = F(1 + 2*i)
1474 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1475 sage: ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1484 Embedding is a homomorphism (isomorphism, in fact)::
1486 sage: set_random_seed()
1487 sage: n = ZZ.random_element(3)
1488 sage: F = QuadraticField(-1, 'I')
1489 sage: X = random_matrix(F, n)
1490 sage: Y = random_matrix(F, n)
1491 sage: Xe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X)
1492 sage: Ye = ComplexMatrixEuclideanJordanAlgebra.real_embed(Y)
1493 sage: XYe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1500 raise ValueError("the matrix 'M' must be square")
1502 # We don't need any adjoined elements...
1503 field
= M
.base_ring().base_ring()
1507 a
= z
.list()[0] # real part, I guess
1508 b
= z
.list()[1] # imag part, I guess
1509 blocks
.append(matrix(field
, 2, [[a
,b
],[-b
,a
]]))
1511 return matrix
.block(field
, n
, blocks
)
1515 def real_unembed(M
):
1517 The inverse of _embed_complex_matrix().
1521 sage: from mjo.eja.eja_algebra import \
1522 ....: ComplexMatrixEuclideanJordanAlgebra
1526 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1527 ....: [-2, 1, -4, 3],
1528 ....: [ 9, 10, 11, 12],
1529 ....: [-10, 9, -12, 11] ])
1530 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(A)
1532 [ 10*I + 9 12*I + 11]
1536 Unembedding is the inverse of embedding::
1538 sage: set_random_seed()
1539 sage: F = QuadraticField(-1, 'I')
1540 sage: M = random_matrix(F, 3)
1541 sage: Me = ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1542 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1548 raise ValueError("the matrix 'M' must be square")
1549 if not n
.mod(2).is_zero():
1550 raise ValueError("the matrix 'M' must be a complex embedding")
1552 # If "M" was normalized, its base ring might have roots
1553 # adjoined and they can stick around after unembedding.
1554 field
= M
.base_ring()
1555 R
= PolynomialRing(field
, 'z')
1558 # Sage doesn't know how to embed AA into QQbar, i.e. how
1559 # to adjoin sqrt(-1) to AA.
1562 F
= field
.extension(z
**2 + 1, 'I', embedding
=CLF(-1).sqrt())
1565 # Go top-left to bottom-right (reading order), converting every
1566 # 2-by-2 block we see to a single complex element.
1568 for k
in range(n
/2):
1569 for j
in range(n
/2):
1570 submat
= M
[2*k
:2*k
+2,2*j
:2*j
+2]
1571 if submat
[0,0] != submat
[1,1]:
1572 raise ValueError('bad on-diagonal submatrix')
1573 if submat
[0,1] != -submat
[1,0]:
1574 raise ValueError('bad off-diagonal submatrix')
1575 z
= submat
[0,0] + submat
[0,1]*i
1578 return matrix(F
, n
/2, elements
)
1582 def natural_inner_product(cls
,X
,Y
):
1584 Compute a natural inner product in this algebra directly from
1589 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1593 This gives the same answer as the slow, default method implemented
1594 in :class:`MatrixEuclideanJordanAlgebra`::
1596 sage: set_random_seed()
1597 sage: J = ComplexHermitianEJA.random_instance()
1598 sage: x,y = J.random_elements(2)
1599 sage: Xe = x.natural_representation()
1600 sage: Ye = y.natural_representation()
1601 sage: X = ComplexHermitianEJA.real_unembed(Xe)
1602 sage: Y = ComplexHermitianEJA.real_unembed(Ye)
1603 sage: expected = (X*Y).trace().real()
1604 sage: actual = ComplexHermitianEJA.natural_inner_product(Xe,Ye)
1605 sage: actual == expected
1609 return RealMatrixEuclideanJordanAlgebra
.natural_inner_product(X
,Y
)/2
1612 class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra
,
1613 ConcreteEuclideanJordanAlgebra
):
1615 The rank-n simple EJA consisting of complex Hermitian n-by-n
1616 matrices over the real numbers, the usual symmetric Jordan product,
1617 and the real-part-of-trace inner product. It has dimension `n^2` over
1622 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1626 In theory, our "field" can be any subfield of the reals::
1628 sage: ComplexHermitianEJA(2, RDF)
1629 Euclidean Jordan algebra of dimension 4 over Real Double Field
1630 sage: ComplexHermitianEJA(2, RR)
1631 Euclidean Jordan algebra of dimension 4 over Real Field with
1632 53 bits of precision
1636 The dimension of this algebra is `n^2`::
1638 sage: set_random_seed()
1639 sage: n_max = ComplexHermitianEJA._max_random_instance_size()
1640 sage: n = ZZ.random_element(1, n_max)
1641 sage: J = ComplexHermitianEJA(n)
1642 sage: J.dimension() == n^2
1645 The Jordan multiplication is what we think it is::
1647 sage: set_random_seed()
1648 sage: J = ComplexHermitianEJA.random_instance()
1649 sage: x,y = J.random_elements(2)
1650 sage: actual = (x*y).natural_representation()
1651 sage: X = x.natural_representation()
1652 sage: Y = y.natural_representation()
1653 sage: expected = (X*Y + Y*X)/2
1654 sage: actual == expected
1656 sage: J(expected) == x*y
1659 We can change the generator prefix::
1661 sage: ComplexHermitianEJA(2, prefix='z').gens()
1664 We can construct the (trivial) algebra of rank zero::
1666 sage: ComplexHermitianEJA(0)
1667 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1672 def _denormalized_basis(cls
, n
, field
):
1674 Returns a basis for the space of complex Hermitian n-by-n matrices.
1676 Why do we embed these? Basically, because all of numerical linear
1677 algebra assumes that you're working with vectors consisting of `n`
1678 entries from a field and scalars from the same field. There's no way
1679 to tell SageMath that (for example) the vectors contain complex
1680 numbers, while the scalar field is real.
1684 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1688 sage: set_random_seed()
1689 sage: n = ZZ.random_element(1,5)
1690 sage: field = QuadraticField(2, 'sqrt2')
1691 sage: B = ComplexHermitianEJA._denormalized_basis(n, field)
1692 sage: all( M.is_symmetric() for M in B)
1696 R
= PolynomialRing(field
, 'z')
1698 F
= field
.extension(z
**2 + 1, 'I')
1701 # This is like the symmetric case, but we need to be careful:
1703 # * We want conjugate-symmetry, not just symmetry.
1704 # * The diagonal will (as a result) be real.
1708 for j
in range(i
+1):
1709 Eij
= matrix(F
, n
, lambda k
,l
: k
==i
and l
==j
)
1711 Sij
= cls
.real_embed(Eij
)
1714 # The second one has a minus because it's conjugated.
1715 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
1717 Sij_imag
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
1720 # Since we embedded these, we can drop back to the "field" that we
1721 # started with instead of the complex extension "F".
1722 return ( s
.change_ring(field
) for s
in S
)
1725 def __init__(self
, n
, field
=AA
, **kwargs
):
1726 basis
= self
._denormalized
_basis
(n
,field
)
1727 super(ComplexHermitianEJA
,self
).__init
__(field
,
1731 self
.rank
.set_cache(n
)
1734 def _max_random_instance_size():
1735 return 3 # Dimension 9
1738 def random_instance(cls
, field
=AA
, **kwargs
):
1740 Return a random instance of this type of algebra.
1742 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1743 return cls(n
, field
, **kwargs
)
1745 class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1749 Embed the n-by-n quaternion matrix ``M`` into the space of real
1750 matrices of size 4n-by-4n by first sending each quaternion entry `z
1751 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
1752 c+di],[-c + di, a-bi]]`, and then embedding those into a real
1757 sage: from mjo.eja.eja_algebra import \
1758 ....: QuaternionMatrixEuclideanJordanAlgebra
1762 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1763 sage: i,j,k = Q.gens()
1764 sage: x = 1 + 2*i + 3*j + 4*k
1765 sage: M = matrix(Q, 1, [[x]])
1766 sage: QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1772 Embedding is a homomorphism (isomorphism, in fact)::
1774 sage: set_random_seed()
1775 sage: n = ZZ.random_element(2)
1776 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1777 sage: X = random_matrix(Q, n)
1778 sage: Y = random_matrix(Q, n)
1779 sage: Xe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X)
1780 sage: Ye = QuaternionMatrixEuclideanJordanAlgebra.real_embed(Y)
1781 sage: XYe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1786 quaternions
= M
.base_ring()
1789 raise ValueError("the matrix 'M' must be square")
1791 F
= QuadraticField(-1, 'I')
1796 t
= z
.coefficient_tuple()
1801 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
1802 [-c
+ d
*i
, a
- b
*i
]])
1803 realM
= ComplexMatrixEuclideanJordanAlgebra
.real_embed(cplxM
)
1804 blocks
.append(realM
)
1806 # We should have real entries by now, so use the realest field
1807 # we've got for the return value.
1808 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
1813 def real_unembed(M
):
1815 The inverse of _embed_quaternion_matrix().
1819 sage: from mjo.eja.eja_algebra import \
1820 ....: QuaternionMatrixEuclideanJordanAlgebra
1824 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
1825 ....: [-2, 1, -4, 3],
1826 ....: [-3, 4, 1, -2],
1827 ....: [-4, -3, 2, 1]])
1828 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(M)
1829 [1 + 2*i + 3*j + 4*k]
1833 Unembedding is the inverse of embedding::
1835 sage: set_random_seed()
1836 sage: Q = QuaternionAlgebra(QQ, -1, -1)
1837 sage: M = random_matrix(Q, 3)
1838 sage: Me = QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1839 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1845 raise ValueError("the matrix 'M' must be square")
1846 if not n
.mod(4).is_zero():
1847 raise ValueError("the matrix 'M' must be a quaternion embedding")
1849 # Use the base ring of the matrix to ensure that its entries can be
1850 # multiplied by elements of the quaternion algebra.
1851 field
= M
.base_ring()
1852 Q
= QuaternionAlgebra(field
,-1,-1)
1855 # Go top-left to bottom-right (reading order), converting every
1856 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
1859 for l
in range(n
/4):
1860 for m
in range(n
/4):
1861 submat
= ComplexMatrixEuclideanJordanAlgebra
.real_unembed(
1862 M
[4*l
:4*l
+4,4*m
:4*m
+4] )
1863 if submat
[0,0] != submat
[1,1].conjugate():
1864 raise ValueError('bad on-diagonal submatrix')
1865 if submat
[0,1] != -submat
[1,0].conjugate():
1866 raise ValueError('bad off-diagonal submatrix')
1867 z
= submat
[0,0].real()
1868 z
+= submat
[0,0].imag()*i
1869 z
+= submat
[0,1].real()*j
1870 z
+= submat
[0,1].imag()*k
1873 return matrix(Q
, n
/4, elements
)
1877 def natural_inner_product(cls
,X
,Y
):
1879 Compute a natural inner product in this algebra directly from
1884 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1888 This gives the same answer as the slow, default method implemented
1889 in :class:`MatrixEuclideanJordanAlgebra`::
1891 sage: set_random_seed()
1892 sage: J = QuaternionHermitianEJA.random_instance()
1893 sage: x,y = J.random_elements(2)
1894 sage: Xe = x.natural_representation()
1895 sage: Ye = y.natural_representation()
1896 sage: X = QuaternionHermitianEJA.real_unembed(Xe)
1897 sage: Y = QuaternionHermitianEJA.real_unembed(Ye)
1898 sage: expected = (X*Y).trace().coefficient_tuple()[0]
1899 sage: actual = QuaternionHermitianEJA.natural_inner_product(Xe,Ye)
1900 sage: actual == expected
1904 return RealMatrixEuclideanJordanAlgebra
.natural_inner_product(X
,Y
)/4
1907 class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra
,
1908 ConcreteEuclideanJordanAlgebra
):
1910 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
1911 matrices, the usual symmetric Jordan product, and the
1912 real-part-of-trace inner product. It has dimension `2n^2 - n` over
1917 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1921 In theory, our "field" can be any subfield of the reals::
1923 sage: QuaternionHermitianEJA(2, RDF)
1924 Euclidean Jordan algebra of dimension 6 over Real Double Field
1925 sage: QuaternionHermitianEJA(2, RR)
1926 Euclidean Jordan algebra of dimension 6 over Real Field with
1927 53 bits of precision
1931 The dimension of this algebra is `2*n^2 - n`::
1933 sage: set_random_seed()
1934 sage: n_max = QuaternionHermitianEJA._max_random_instance_size()
1935 sage: n = ZZ.random_element(1, n_max)
1936 sage: J = QuaternionHermitianEJA(n)
1937 sage: J.dimension() == 2*(n^2) - n
1940 The Jordan multiplication is what we think it is::
1942 sage: set_random_seed()
1943 sage: J = QuaternionHermitianEJA.random_instance()
1944 sage: x,y = J.random_elements(2)
1945 sage: actual = (x*y).natural_representation()
1946 sage: X = x.natural_representation()
1947 sage: Y = y.natural_representation()
1948 sage: expected = (X*Y + Y*X)/2
1949 sage: actual == expected
1951 sage: J(expected) == x*y
1954 We can change the generator prefix::
1956 sage: QuaternionHermitianEJA(2, prefix='a').gens()
1957 (a0, a1, a2, a3, a4, a5)
1959 We can construct the (trivial) algebra of rank zero::
1961 sage: QuaternionHermitianEJA(0)
1962 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1966 def _denormalized_basis(cls
, n
, field
):
1968 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
1970 Why do we embed these? Basically, because all of numerical
1971 linear algebra assumes that you're working with vectors consisting
1972 of `n` entries from a field and scalars from the same field. There's
1973 no way to tell SageMath that (for example) the vectors contain
1974 complex numbers, while the scalar field is real.
1978 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1982 sage: set_random_seed()
1983 sage: n = ZZ.random_element(1,5)
1984 sage: B = QuaternionHermitianEJA._denormalized_basis(n,QQ)
1985 sage: all( M.is_symmetric() for M in B )
1989 Q
= QuaternionAlgebra(QQ
,-1,-1)
1992 # This is like the symmetric case, but we need to be careful:
1994 # * We want conjugate-symmetry, not just symmetry.
1995 # * The diagonal will (as a result) be real.
1999 for j
in range(i
+1):
2000 Eij
= matrix(Q
, n
, lambda k
,l
: k
==i
and l
==j
)
2002 Sij
= cls
.real_embed(Eij
)
2005 # The second, third, and fourth ones have a minus
2006 # because they're conjugated.
2007 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
2009 Sij_I
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
2011 Sij_J
= cls
.real_embed(J
*Eij
- J
*Eij
.transpose())
2013 Sij_K
= cls
.real_embed(K
*Eij
- K
*Eij
.transpose())
2016 # Since we embedded these, we can drop back to the "field" that we
2017 # started with instead of the quaternion algebra "Q".
2018 return ( s
.change_ring(field
) for s
in S
)
2021 def __init__(self
, n
, field
=AA
, **kwargs
):
2022 basis
= self
._denormalized
_basis
(n
,field
)
2023 super(QuaternionHermitianEJA
,self
).__init
__(field
,
2027 self
.rank
.set_cache(n
)
2030 def _max_random_instance_size():
2032 The maximum rank of a random QuaternionHermitianEJA.
2034 return 2 # Dimension 6
2037 def random_instance(cls
, field
=AA
, **kwargs
):
2039 Return a random instance of this type of algebra.
2041 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2042 return cls(n
, field
, **kwargs
)
2045 class HadamardEJA(RationalBasisEuclideanJordanAlgebra
,
2046 ConcreteEuclideanJordanAlgebra
):
2048 Return the Euclidean Jordan Algebra corresponding to the set
2049 `R^n` under the Hadamard product.
2051 Note: this is nothing more than the Cartesian product of ``n``
2052 copies of the spin algebra. Once Cartesian product algebras
2053 are implemented, this can go.
2057 sage: from mjo.eja.eja_algebra import HadamardEJA
2061 This multiplication table can be verified by hand::
2063 sage: J = HadamardEJA(3)
2064 sage: e0,e1,e2 = J.gens()
2080 We can change the generator prefix::
2082 sage: HadamardEJA(3, prefix='r').gens()
2086 def __init__(self
, n
, field
=AA
, **kwargs
):
2087 V
= VectorSpace(field
, n
)
2088 mult_table
= [ [ V
.gen(i
)*(i
== j
) for j
in range(n
) ]
2091 # Inner products are real numbers and not algebra
2092 # elements, so once we turn the algebra element
2093 # into a vector in inner_product(), we never go
2094 # back. As a result -- contrary to what we do with
2095 # self._multiplication_table -- we store the inner
2096 # product table as a plain old matrix and not as
2097 # an algebra operator.
2098 ip_table
= matrix
.identity(field
,n
)
2099 self
._inner
_product
_matrix
= ip_table
2101 super(HadamardEJA
, self
).__init
__(field
,
2105 self
.rank
.set_cache(n
)
2108 self
.one
.set_cache( self
.zero() )
2110 self
.one
.set_cache( sum(self
.gens()) )
2113 def _max_random_instance_size():
2115 The maximum dimension of a random HadamardEJA.
2120 def random_instance(cls
, field
=AA
, **kwargs
):
2122 Return a random instance of this type of algebra.
2124 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2125 return cls(n
, field
, **kwargs
)
2128 class BilinearFormEJA(RationalBasisEuclideanJordanAlgebra
,
2129 ConcreteEuclideanJordanAlgebra
):
2131 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2132 with the half-trace inner product and jordan product ``x*y =
2133 (<Bx,y>,y_bar>, x0*y_bar + y0*x_bar)`` where `B = 1 \times B22` is
2134 a symmetric positive-definite "bilinear form" matrix. Its
2135 dimension is the size of `B`, and it has rank two in dimensions
2136 larger than two. It reduces to the ``JordanSpinEJA`` when `B` is
2137 the identity matrix of order ``n``.
2139 We insist that the one-by-one upper-left identity block of `B` be
2140 passed in as well so that we can be passed a matrix of size zero
2141 to construct a trivial algebra.
2145 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
2146 ....: JordanSpinEJA)
2150 When no bilinear form is specified, the identity matrix is used,
2151 and the resulting algebra is the Jordan spin algebra::
2153 sage: B = matrix.identity(AA,3)
2154 sage: J0 = BilinearFormEJA(B)
2155 sage: J1 = JordanSpinEJA(3)
2156 sage: J0.multiplication_table() == J0.multiplication_table()
2159 An error is raised if the matrix `B` does not correspond to a
2160 positive-definite bilinear form::
2162 sage: B = matrix.random(QQ,2,3)
2163 sage: J = BilinearFormEJA(B)
2164 Traceback (most recent call last):
2166 ValueError: bilinear form is not positive-definite
2167 sage: B = matrix.zero(QQ,3)
2168 sage: J = BilinearFormEJA(B)
2169 Traceback (most recent call last):
2171 ValueError: bilinear form is not positive-definite
2175 We can create a zero-dimensional algebra::
2177 sage: B = matrix.identity(AA,0)
2178 sage: J = BilinearFormEJA(B)
2182 We can check the multiplication condition given in the Jordan, von
2183 Neumann, and Wigner paper (and also discussed on my "On the
2184 symmetry..." paper). Note that this relies heavily on the standard
2185 choice of basis, as does anything utilizing the bilinear form matrix::
2187 sage: set_random_seed()
2188 sage: n = ZZ.random_element(5)
2189 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2190 sage: B11 = matrix.identity(QQ,1)
2191 sage: B22 = M.transpose()*M
2192 sage: B = block_matrix(2,2,[ [B11,0 ],
2194 sage: J = BilinearFormEJA(B)
2195 sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
2196 sage: V = J.vector_space()
2197 sage: sis = [ J.from_vector(V([0] + (M.inverse()*ei).list()))
2198 ....: for ei in eis ]
2199 sage: actual = [ sis[i]*sis[j]
2200 ....: for i in range(n-1)
2201 ....: for j in range(n-1) ]
2202 sage: expected = [ J.one() if i == j else J.zero()
2203 ....: for i in range(n-1)
2204 ....: for j in range(n-1) ]
2205 sage: actual == expected
2208 def __init__(self
, B
, field
=AA
, **kwargs
):
2211 if not B
.is_positive_definite():
2212 raise ValueError("bilinear form is not positive-definite")
2214 V
= VectorSpace(field
, n
)
2215 mult_table
= [[V
.zero() for j
in range(n
)] for i
in range(n
)]
2224 z0
= (B
*x
).inner_product(y
)
2225 zbar
= y0
*xbar
+ x0
*ybar
2226 z
= V([z0
] + zbar
.list())
2227 mult_table
[i
][j
] = z
2229 # Inner products are real numbers and not algebra
2230 # elements, so once we turn the algebra element
2231 # into a vector in inner_product(), we never go
2232 # back. As a result -- contrary to what we do with
2233 # self._multiplication_table -- we store the inner
2234 # product table as a plain old matrix and not as
2235 # an algebra operator.
2237 self
._inner
_product
_matrix
= ip_table
2239 super(BilinearFormEJA
, self
).__init
__(field
,
2244 # The rank of this algebra is two, unless we're in a
2245 # one-dimensional ambient space (because the rank is bounded
2246 # by the ambient dimension).
2247 self
.rank
.set_cache(min(n
,2))
2250 self
.one
.set_cache( self
.zero() )
2252 self
.one
.set_cache( self
.monomial(0) )
2255 def _max_random_instance_size():
2257 The maximum dimension of a random BilinearFormEJA.
2262 def random_instance(cls
, field
=AA
, **kwargs
):
2264 Return a random instance of this algebra.
2266 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2268 B
= matrix
.identity(field
, n
)
2269 return cls(B
, field
, **kwargs
)
2271 B11
= matrix
.identity(field
,1)
2272 M
= matrix
.random(field
, n
-1)
2273 I
= matrix
.identity(field
, n
-1)
2274 alpha
= field
.zero()
2275 while alpha
.is_zero():
2276 alpha
= field
.random_element().abs()
2277 B22
= M
.transpose()*M
+ alpha
*I
2279 from sage
.matrix
.special
import block_matrix
2280 B
= block_matrix(2,2, [ [B11
, ZZ(0) ],
2283 return cls(B
, field
, **kwargs
)
2286 class JordanSpinEJA(BilinearFormEJA
):
2288 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2289 with the usual inner product and jordan product ``x*y =
2290 (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2295 sage: from mjo.eja.eja_algebra import JordanSpinEJA
2299 This multiplication table can be verified by hand::
2301 sage: J = JordanSpinEJA(4)
2302 sage: e0,e1,e2,e3 = J.gens()
2318 We can change the generator prefix::
2320 sage: JordanSpinEJA(2, prefix='B').gens()
2325 Ensure that we have the usual inner product on `R^n`::
2327 sage: set_random_seed()
2328 sage: J = JordanSpinEJA.random_instance()
2329 sage: x,y = J.random_elements(2)
2330 sage: actual = x.inner_product(y)
2331 sage: expected = x.to_vector().inner_product(y.to_vector())
2332 sage: actual == expected
2336 def __init__(self
, n
, field
=AA
, **kwargs
):
2337 # This is a special case of the BilinearFormEJA with the identity
2338 # matrix as its bilinear form.
2339 B
= matrix
.identity(field
, n
)
2340 super(JordanSpinEJA
, self
).__init
__(B
, field
, **kwargs
)
2343 def _max_random_instance_size():
2345 The maximum dimension of a random JordanSpinEJA.
2350 def random_instance(cls
, field
=AA
, **kwargs
):
2352 Return a random instance of this type of algebra.
2354 Needed here to override the implementation for ``BilinearFormEJA``.
2356 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2357 return cls(n
, field
, **kwargs
)
2360 class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra
,
2361 ConcreteEuclideanJordanAlgebra
):
2363 The trivial Euclidean Jordan algebra consisting of only a zero element.
2367 sage: from mjo.eja.eja_algebra import TrivialEJA
2371 sage: J = TrivialEJA()
2378 sage: 7*J.one()*12*J.one()
2380 sage: J.one().inner_product(J.one())
2382 sage: J.one().norm()
2384 sage: J.one().subalgebra_generated_by()
2385 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2390 def __init__(self
, field
=AA
, **kwargs
):
2392 self
._inner
_product
_matrix
= matrix(field
,0)
2393 super(TrivialEJA
, self
).__init
__(field
,
2397 # The rank is zero using my definition, namely the dimension of the
2398 # largest subalgebra generated by any element.
2399 self
.rank
.set_cache(0)
2400 self
.one
.set_cache( self
.zero() )
2403 def random_instance(cls
, field
=AA
, **kwargs
):
2404 # We don't take a "size" argument so the superclass method is
2405 # inappropriate for us.
2406 return cls(field
, **kwargs
)
2408 class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra
):
2410 The external (orthogonal) direct sum of two other Euclidean Jordan
2411 algebras. Essentially the Cartesian product of its two factors.
2412 Every Euclidean Jordan algebra decomposes into an orthogonal
2413 direct sum of simple Euclidean Jordan algebras, so no generality
2414 is lost by providing only this construction.
2418 sage: from mjo.eja.eja_algebra import (random_eja,
2420 ....: RealSymmetricEJA,
2425 sage: J1 = HadamardEJA(2)
2426 sage: J2 = RealSymmetricEJA(3)
2427 sage: J = DirectSumEJA(J1,J2)
2435 The external direct sum construction is only valid when the two factors
2436 have the same base ring; an error is raised otherwise::
2438 sage: set_random_seed()
2439 sage: J1 = random_eja(AA)
2440 sage: J2 = random_eja(QQ)
2441 sage: J = DirectSumEJA(J1,J2)
2442 Traceback (most recent call last):
2444 ValueError: algebras must share the same base field
2447 def __init__(self
, J1
, J2
, **kwargs
):
2448 if J1
.base_ring() != J2
.base_ring():
2449 raise ValueError("algebras must share the same base field")
2450 field
= J1
.base_ring()
2452 self
._factors
= (J1
, J2
)
2456 V
= VectorSpace(field
, n
)
2457 mult_table
= [ [ V
.zero() for j
in range(n
) ]
2461 p
= (J1
.monomial(i
)*J1
.monomial(j
)).to_vector()
2462 mult_table
[i
][j
] = V(p
.list() + [field
.zero()]*n2
)
2466 p
= (J2
.monomial(i
)*J2
.monomial(j
)).to_vector()
2467 mult_table
[n1
+i
][n1
+j
] = V([field
.zero()]*n1
+ p
.list())
2469 super(DirectSumEJA
, self
).__init
__(field
,
2473 self
.rank
.set_cache(J1
.rank() + J2
.rank())
2478 Return the pair of this algebra's factors.
2482 sage: from mjo.eja.eja_algebra import (HadamardEJA,
2483 ....: JordanSpinEJA,
2488 sage: J1 = HadamardEJA(2,QQ)
2489 sage: J2 = JordanSpinEJA(3,QQ)
2490 sage: J = DirectSumEJA(J1,J2)
2492 (Euclidean Jordan algebra of dimension 2 over Rational Field,
2493 Euclidean Jordan algebra of dimension 3 over Rational Field)
2496 return self
._factors
2498 def projections(self
):
2500 Return a pair of projections onto this algebra's factors.
2504 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
2505 ....: ComplexHermitianEJA,
2510 sage: J1 = JordanSpinEJA(2)
2511 sage: J2 = ComplexHermitianEJA(2)
2512 sage: J = DirectSumEJA(J1,J2)
2513 sage: (pi_left, pi_right) = J.projections()
2514 sage: J.one().to_vector()
2516 sage: pi_left(J.one()).to_vector()
2518 sage: pi_right(J.one()).to_vector()
2522 (J1
,J2
) = self
.factors()
2525 V_basis
= self
.vector_space().basis()
2526 # Need to specify the dimensions explicitly so that we don't
2527 # wind up with a zero-by-zero matrix when we want e.g. a
2528 # zero-by-two matrix (important for composing things).
2529 P1
= matrix(self
.base_ring(), m
, m
+n
, V_basis
[:m
])
2530 P2
= matrix(self
.base_ring(), n
, m
+n
, V_basis
[m
:])
2531 pi_left
= FiniteDimensionalEuclideanJordanAlgebraOperator(self
,J1
,P1
)
2532 pi_right
= FiniteDimensionalEuclideanJordanAlgebraOperator(self
,J2
,P2
)
2533 return (pi_left
, pi_right
)
2535 def inclusions(self
):
2537 Return the pair of inclusion maps from our factors into us.
2541 sage: from mjo.eja.eja_algebra import (random_eja,
2542 ....: JordanSpinEJA,
2543 ....: RealSymmetricEJA,
2548 sage: J1 = JordanSpinEJA(3)
2549 sage: J2 = RealSymmetricEJA(2)
2550 sage: J = DirectSumEJA(J1,J2)
2551 sage: (iota_left, iota_right) = J.inclusions()
2552 sage: iota_left(J1.zero()) == J.zero()
2554 sage: iota_right(J2.zero()) == J.zero()
2556 sage: J1.one().to_vector()
2558 sage: iota_left(J1.one()).to_vector()
2560 sage: J2.one().to_vector()
2562 sage: iota_right(J2.one()).to_vector()
2564 sage: J.one().to_vector()
2569 Composing a projection with the corresponding inclusion should
2570 produce the identity map, and mismatching them should produce
2573 sage: set_random_seed()
2574 sage: J1 = random_eja()
2575 sage: J2 = random_eja()
2576 sage: J = DirectSumEJA(J1,J2)
2577 sage: (iota_left, iota_right) = J.inclusions()
2578 sage: (pi_left, pi_right) = J.projections()
2579 sage: pi_left*iota_left == J1.one().operator()
2581 sage: pi_right*iota_right == J2.one().operator()
2583 sage: (pi_left*iota_right).is_zero()
2585 sage: (pi_right*iota_left).is_zero()
2589 (J1
,J2
) = self
.factors()
2592 V_basis
= self
.vector_space().basis()
2593 # Need to specify the dimensions explicitly so that we don't
2594 # wind up with a zero-by-zero matrix when we want e.g. a
2595 # two-by-zero matrix (important for composing things).
2596 I1
= matrix
.column(self
.base_ring(), m
, m
+n
, V_basis
[:m
])
2597 I2
= matrix
.column(self
.base_ring(), n
, m
+n
, V_basis
[m
:])
2598 iota_left
= FiniteDimensionalEuclideanJordanAlgebraOperator(J1
,self
,I1
)
2599 iota_right
= FiniteDimensionalEuclideanJordanAlgebraOperator(J2
,self
,I2
)
2600 return (iota_left
, iota_right
)
2602 def inner_product(self
, x
, y
):
2604 The standard Cartesian inner-product.
2606 We project ``x`` and ``y`` onto our factors, and add up the
2607 inner-products from the subalgebras.
2612 sage: from mjo.eja.eja_algebra import (HadamardEJA,
2613 ....: QuaternionHermitianEJA,
2618 sage: J1 = HadamardEJA(3,QQ)
2619 sage: J2 = QuaternionHermitianEJA(2,QQ,normalize_basis=False)
2620 sage: J = DirectSumEJA(J1,J2)
2625 sage: x1.inner_product(x2)
2627 sage: y1.inner_product(y2)
2629 sage: J.one().inner_product(J.one())
2633 (pi_left
, pi_right
) = self
.projections()
2639 return (x1
.inner_product(y1
) + x2
.inner_product(y2
))
2643 random_eja
= ConcreteEuclideanJordanAlgebra
.random_instance