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 The lowest-level class for representing a Euclidean Jordan algebra.
40 def _coerce_map_from_base_ring(self
):
42 Disable the map from the base ring into the algebra.
44 Performing a nonsense conversion like this automatically
45 is counterpedagogical. The fallback is to try the usual
46 element constructor, which should also fail.
50 sage: from mjo.eja.eja_algebra import random_eja
54 sage: set_random_seed()
55 sage: J = random_eja()
57 Traceback (most recent call last):
59 ValueError: not an element of this algebra
76 * field -- the scalar field for this algebra (must be real)
78 * multiplication_table -- the multiplication table for this
79 algebra's implicit basis. Only the lower-triangular portion
80 of the table is used, since the multiplication is assumed
85 sage: from mjo.eja.eja_algebra import (
86 ....: FiniteDimensionalEuclideanJordanAlgebra,
92 By definition, Jordan multiplication commutes::
94 sage: set_random_seed()
95 sage: J = random_eja()
96 sage: x,y = J.random_elements(2)
100 An error is raised if the Jordan product is not commutative::
102 sage: JP = ((1,2),(0,0))
103 sage: IP = ((1,0),(0,1))
104 sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,JP,IP)
105 Traceback (most recent call last):
107 ValueError: Jordan product is not commutative
109 An error is raised if the inner-product is not commutative::
111 sage: JP = ((1,0),(0,1))
112 sage: IP = ((1,2),(0,0))
113 sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,JP,IP)
114 Traceback (most recent call last):
116 ValueError: inner-product is not commutative
120 The ``field`` we're given must be real with ``check_field=True``::
122 sage: JordanSpinEJA(2, QQbar)
123 Traceback (most recent call last):
125 ValueError: scalar field is not real
126 sage: JordanSpinEJA(2, QQbar, check_field=False)
127 Euclidean Jordan algebra of dimension 2 over Algebraic Field
129 The multiplication table must be square with ``check_axioms=True``::
131 sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,((),()),((1,),))
132 Traceback (most recent call last):
134 ValueError: multiplication table is not square
136 The multiplication and inner-product tables must be the same
137 size (and in particular, the inner-product table must also be
138 square) with ``check_axioms=True``::
140 sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,((1,),),(()))
141 Traceback (most recent call last):
143 ValueError: multiplication and inner-product tables are
145 sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,((1,),),((1,2),))
146 Traceback (most recent call last):
148 ValueError: multiplication and inner-product tables are
153 if not field
.is_subring(RR
):
154 # Note: this does return true for the real algebraic
155 # field, the rationals, and any quadratic field where
156 # we've specified a real embedding.
157 raise ValueError("scalar field is not real")
160 # The multiplication and inner-product tables should be square
161 # if the user wants us to verify them. And we verify them as
162 # soon as possible, because we want to exploit their symmetry.
163 n
= len(multiplication_table
)
165 if not all( len(l
) == n
for l
in multiplication_table
):
166 raise ValueError("multiplication table is not square")
168 # If the multiplication table is square, we can check if
169 # the inner-product table is square by comparing it to the
170 # multiplication table's dimensions.
171 msg
= "multiplication and inner-product tables are different sizes"
172 if not len(inner_product_table
) == n
:
173 raise ValueError(msg
)
175 if not all( len(l
) == n
for l
in inner_product_table
):
176 raise ValueError(msg
)
178 # Check commutativity of the Jordan product (symmetry of
179 # the multiplication table) and the commutativity of the
180 # inner-product (symmetry of the inner-product table)
181 # first if we're going to check them at all.. This has to
182 # be done before we define product_on_basis(), because
183 # that method assumes that self._multiplication_table is
184 # symmetric. And it has to be done before we build
185 # self._inner_product_matrix, because the process used to
186 # construct it assumes symmetry as well.
187 if not all( multiplication_table
[j
][i
]
188 == multiplication_table
[i
][j
]
190 for j
in range(i
+1) ):
191 raise ValueError("Jordan product is not commutative")
193 if not all( inner_product_table
[j
][i
]
194 == inner_product_table
[i
][j
]
196 for j
in range(i
+1) ):
197 raise ValueError("inner-product is not commutative")
199 self
._matrix
_basis
= matrix_basis
202 category
= MagmaticAlgebras(field
).FiniteDimensional()
203 category
= category
.WithBasis().Unital()
205 fda
= super(FiniteDimensionalEuclideanJordanAlgebra
, self
)
210 self
.print_options(bracket
='')
212 # The multiplication table we're given is necessarily in terms
213 # of vectors, because we don't have an algebra yet for
214 # anything to be an element of. However, it's faster in the
215 # long run to have the multiplication table be in terms of
216 # algebra elements. We do this after calling the superclass
217 # constructor so that from_vector() knows what to do.
219 # Note: we take advantage of symmetry here, and only store
220 # the lower-triangular portion of the table.
221 self
._multiplication
_table
= [ [ self
.vector_space().zero()
222 for j
in range(i
+1) ]
227 elt
= self
.from_vector(multiplication_table
[i
][j
])
228 self
._multiplication
_table
[i
][j
] = elt
230 self
._multiplication
_table
= tuple(map(tuple, self
._multiplication
_table
))
232 # Save our inner product as a matrix, since the efficiency of
233 # matrix multiplication will usually outweigh the fact that we
234 # have to store a redundant upper- or lower-triangular part.
235 # Pre-cache the fact that these are Hermitian (real symmetric,
236 # in fact) in case some e.g. matrix multiplication routine can
237 # take advantage of it.
238 ip_matrix_constructor
= lambda i
,j
: inner_product_table
[i
][j
] if j
<= i
else inner_product_table
[j
][i
]
239 self
._inner
_product
_matrix
= matrix(field
, n
, ip_matrix_constructor
)
240 self
._inner
_product
_matrix
._cache
= {'hermitian': True}
241 self
._inner
_product
_matrix
.set_immutable()
244 if not self
._is
_jordanian
():
245 raise ValueError("Jordan identity does not hold")
246 if not self
._inner
_product
_is
_associative
():
247 raise ValueError("inner product is not associative")
249 def _element_constructor_(self
, elt
):
251 Construct an element of this algebra from its vector or matrix
254 This gets called only after the parent element _call_ method
255 fails to find a coercion for the argument.
259 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
261 ....: RealSymmetricEJA)
265 The identity in `S^n` is converted to the identity in the EJA::
267 sage: J = RealSymmetricEJA(3)
268 sage: I = matrix.identity(QQ,3)
269 sage: J(I) == J.one()
272 This skew-symmetric matrix can't be represented in the EJA::
274 sage: J = RealSymmetricEJA(3)
275 sage: A = matrix(QQ,3, lambda i,j: i-j)
277 Traceback (most recent call last):
279 ValueError: not an element of this algebra
283 Ensure that we can convert any element of the two non-matrix
284 simple algebras (whose matrix representations are columns)
285 back and forth faithfully::
287 sage: set_random_seed()
288 sage: J = HadamardEJA.random_instance()
289 sage: x = J.random_element()
290 sage: J(x.to_vector().column()) == x
292 sage: J = JordanSpinEJA.random_instance()
293 sage: x = J.random_element()
294 sage: J(x.to_vector().column()) == x
297 msg
= "not an element of this algebra"
299 # The superclass implementation of random_element()
300 # needs to be able to coerce "0" into the algebra.
302 elif elt
in self
.base_ring():
303 # Ensure that no base ring -> algebra coercion is performed
304 # by this method. There's some stupidity in sage that would
305 # otherwise propagate to this method; for example, sage thinks
306 # that the integer 3 belongs to the space of 2-by-2 matrices.
307 raise ValueError(msg
)
309 if elt
not in self
.matrix_space():
310 raise ValueError(msg
)
312 # Thanks for nothing! Matrix spaces aren't vector spaces in
313 # Sage, so we have to figure out its matrix-basis coordinates
314 # ourselves. We use the basis space's ring instead of the
315 # element's ring because the basis space might be an algebraic
316 # closure whereas the base ring of the 3-by-3 identity matrix
317 # could be QQ instead of QQbar.
319 # We pass check=False because the matrix basis is "guaranteed"
320 # to be linearly independent... right? Ha ha.
321 V
= VectorSpace(self
.base_ring(), elt
.nrows()*elt
.ncols())
322 W
= V
.span_of_basis( (_mat2vec(s
) for s
in self
.matrix_basis()),
326 coords
= W
.coordinate_vector(_mat2vec(elt
))
327 except ArithmeticError: # vector is not in free module
328 raise ValueError(msg
)
330 return self
.from_vector(coords
)
334 Return a string representation of ``self``.
338 sage: from mjo.eja.eja_algebra import JordanSpinEJA
342 Ensure that it says what we think it says::
344 sage: JordanSpinEJA(2, field=AA)
345 Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
346 sage: JordanSpinEJA(3, field=RDF)
347 Euclidean Jordan algebra of dimension 3 over Real Double Field
350 fmt
= "Euclidean Jordan algebra of dimension {} over {}"
351 return fmt
.format(self
.dimension(), self
.base_ring())
353 def product_on_basis(self
, i
, j
):
354 # We only stored the lower-triangular portion of the
355 # multiplication table.
357 return self
._multiplication
_table
[i
][j
]
359 return self
._multiplication
_table
[j
][i
]
361 def _is_commutative(self
):
363 Whether or not this algebra's multiplication table is commutative.
365 This method should of course always return ``True``, unless
366 this algebra was constructed with ``check_axioms=False`` and
367 passed an invalid multiplication table.
369 return all( self
.product_on_basis(i
,j
) == self
.product_on_basis(i
,j
)
370 for i
in range(self
.dimension())
371 for j
in range(self
.dimension()) )
373 def _is_jordanian(self
):
375 Whether or not this algebra's multiplication table respects the
376 Jordan identity `(x^{2})(xy) = x(x^{2}y)`.
378 We only check one arrangement of `x` and `y`, so for a
379 ``True`` result to be truly true, you should also check
380 :meth:`_is_commutative`. This method should of course always
381 return ``True``, unless this algebra was constructed with
382 ``check_axioms=False`` and passed an invalid multiplication table.
384 return all( (self
.monomial(i
)**2)*(self
.monomial(i
)*self
.monomial(j
))
386 (self
.monomial(i
))*((self
.monomial(i
)**2)*self
.monomial(j
))
387 for i
in range(self
.dimension())
388 for j
in range(self
.dimension()) )
390 def _inner_product_is_associative(self
):
392 Return whether or not this algebra's inner product `B` is
393 associative; that is, whether or not `B(xy,z) = B(x,yz)`.
395 This method should of course always return ``True``, unless
396 this algebra was constructed with ``check_axioms=False`` and
397 passed an invalid multiplication table.
400 # Used to check whether or not something is zero in an inexact
401 # ring. This number is sufficient to allow the construction of
402 # QuaternionHermitianEJA(2, RDF) with check_axioms=True.
405 for i
in range(self
.dimension()):
406 for j
in range(self
.dimension()):
407 for k
in range(self
.dimension()):
411 diff
= (x
*y
).inner_product(z
) - x
.inner_product(y
*z
)
413 if self
.base_ring().is_exact():
417 if diff
.abs() > epsilon
:
423 def characteristic_polynomial_of(self
):
425 Return the algebra's "characteristic polynomial of" function,
426 which is itself a multivariate polynomial that, when evaluated
427 at the coordinates of some algebra element, returns that
428 element's characteristic polynomial.
430 The resulting polynomial has `n+1` variables, where `n` is the
431 dimension of this algebra. The first `n` variables correspond to
432 the coordinates of an algebra element: when evaluated at the
433 coordinates of an algebra element with respect to a certain
434 basis, the result is a univariate polynomial (in the one
435 remaining variable ``t``), namely the characteristic polynomial
440 sage: from mjo.eja.eja_algebra import JordanSpinEJA, TrivialEJA
444 The characteristic polynomial in the spin algebra is given in
445 Alizadeh, Example 11.11::
447 sage: J = JordanSpinEJA(3)
448 sage: p = J.characteristic_polynomial_of(); p
449 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
450 sage: xvec = J.one().to_vector()
454 By definition, the characteristic polynomial is a monic
455 degree-zero polynomial in a rank-zero algebra. Note that
456 Cayley-Hamilton is indeed satisfied since the polynomial
457 ``1`` evaluates to the identity element of the algebra on
460 sage: J = TrivialEJA()
461 sage: J.characteristic_polynomial_of()
468 # The list of coefficient polynomials a_0, a_1, a_2, ..., a_(r-1).
469 a
= self
._charpoly
_coefficients
()
471 # We go to a bit of trouble here to reorder the
472 # indeterminates, so that it's easier to evaluate the
473 # characteristic polynomial at x's coordinates and get back
474 # something in terms of t, which is what we want.
475 S
= PolynomialRing(self
.base_ring(),'t')
479 S
= PolynomialRing(S
, R
.variable_names())
482 return (t
**r
+ sum( a
[k
]*(t
**k
) for k
in range(r
) ))
484 def coordinate_polynomial_ring(self
):
486 The multivariate polynomial ring in which this algebra's
487 :meth:`characteristic_polynomial_of` lives.
491 sage: from mjo.eja.eja_algebra import (HadamardEJA,
492 ....: RealSymmetricEJA)
496 sage: J = HadamardEJA(2)
497 sage: J.coordinate_polynomial_ring()
498 Multivariate Polynomial Ring in X1, X2...
499 sage: J = RealSymmetricEJA(3,QQ,orthonormalize=False)
500 sage: J.coordinate_polynomial_ring()
501 Multivariate Polynomial Ring in X1, X2, X3, X4, X5, X6...
504 var_names
= tuple( "X%d" % z
for z
in range(1, self
.dimension()+1) )
505 return PolynomialRing(self
.base_ring(), var_names
)
507 def inner_product(self
, x
, y
):
509 The inner product associated with this Euclidean Jordan algebra.
511 Defaults to the trace inner product, but can be overridden by
512 subclasses if they are sure that the necessary properties are
517 sage: from mjo.eja.eja_algebra import (random_eja,
519 ....: BilinearFormEJA)
523 Our inner product is "associative," which means the following for
524 a symmetric bilinear form::
526 sage: set_random_seed()
527 sage: J = random_eja()
528 sage: x,y,z = J.random_elements(3)
529 sage: (x*y).inner_product(z) == y.inner_product(x*z)
534 Ensure that this is the usual inner product for the algebras
537 sage: set_random_seed()
538 sage: J = HadamardEJA.random_instance()
539 sage: x,y = J.random_elements(2)
540 sage: actual = x.inner_product(y)
541 sage: expected = x.to_vector().inner_product(y.to_vector())
542 sage: actual == expected
545 Ensure that this is one-half of the trace inner-product in a
546 BilinearFormEJA that isn't just the reals (when ``n`` isn't
547 one). This is in Faraut and Koranyi, and also my "On the
550 sage: set_random_seed()
551 sage: J = BilinearFormEJA.random_instance()
552 sage: n = J.dimension()
553 sage: x = J.random_element()
554 sage: y = J.random_element()
555 sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
558 B
= self
._inner
_product
_matrix
559 return (B
*x
.to_vector()).inner_product(y
.to_vector())
562 def is_trivial(self
):
564 Return whether or not this algebra is trivial.
566 A trivial algebra contains only the zero element.
570 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
575 sage: J = ComplexHermitianEJA(3)
581 sage: J = TrivialEJA()
586 return self
.dimension() == 0
589 def multiplication_table(self
):
591 Return a visual representation of this algebra's multiplication
592 table (on basis elements).
596 sage: from mjo.eja.eja_algebra import JordanSpinEJA
600 sage: J = JordanSpinEJA(4)
601 sage: J.multiplication_table()
602 +----++----+----+----+----+
603 | * || e0 | e1 | e2 | e3 |
604 +====++====+====+====+====+
605 | e0 || e0 | e1 | e2 | e3 |
606 +----++----+----+----+----+
607 | e1 || e1 | e0 | 0 | 0 |
608 +----++----+----+----+----+
609 | e2 || e2 | 0 | e0 | 0 |
610 +----++----+----+----+----+
611 | e3 || e3 | 0 | 0 | e0 |
612 +----++----+----+----+----+
616 M
= [ [ self
.zero() for j
in range(n
) ]
620 M
[i
][j
] = self
._multiplication
_table
[i
][j
]
624 # Prepend the left "header" column entry Can't do this in
625 # the loop because it messes up the symmetry.
626 M
[i
] = [self
.monomial(i
)] + M
[i
]
628 # Prepend the header row.
629 M
= [["*"] + list(self
.gens())] + M
630 return table(M
, header_row
=True, header_column
=True, frame
=True)
633 def matrix_basis(self
):
635 Return an (often more natural) representation of this algebras
636 basis as an ordered tuple of matrices.
638 Every finite-dimensional Euclidean Jordan Algebra is a, up to
639 Jordan isomorphism, a direct sum of five simple
640 algebras---four of which comprise Hermitian matrices. And the
641 last type of algebra can of course be thought of as `n`-by-`1`
642 column matrices (ambiguusly called column vectors) to avoid
643 special cases. As a result, matrices (and column vectors) are
644 a natural representation format for Euclidean Jordan algebra
647 But, when we construct an algebra from a basis of matrices,
648 those matrix representations are lost in favor of coordinate
649 vectors *with respect to* that basis. We could eventually
650 convert back if we tried hard enough, but having the original
651 representations handy is valuable enough that we simply store
652 them and return them from this method.
654 Why implement this for non-matrix algebras? Avoiding special
655 cases for the :class:`BilinearFormEJA` pays with simplicity in
656 its own right. But mainly, we would like to be able to assume
657 that elements of a :class:`DirectSumEJA` can be displayed
658 nicely, without having to have special classes for direct sums
659 one of whose components was a matrix algebra.
663 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
664 ....: RealSymmetricEJA)
668 sage: J = RealSymmetricEJA(2)
670 Finite family {0: e0, 1: e1, 2: e2}
671 sage: J.matrix_basis()
673 [1 0] [ 0 0.7071067811865475?] [0 0]
674 [0 0], [0.7071067811865475? 0], [0 1]
679 sage: J = JordanSpinEJA(2)
681 Finite family {0: e0, 1: e1}
682 sage: J.matrix_basis()
688 if self
._matrix
_basis
is None:
689 M
= self
.matrix_space()
690 return tuple( M(b
.to_vector()) for b
in self
.basis() )
692 return self
._matrix
_basis
695 def matrix_space(self
):
697 Return the matrix space in which this algebra's elements live, if
698 we think of them as matrices (including column vectors of the
701 Generally this will be an `n`-by-`1` column-vector space,
702 except when the algebra is trivial. There it's `n`-by-`n`
703 (where `n` is zero), to ensure that two elements of the matrix
704 space (empty matrices) can be multiplied.
706 Matrix algebras override this with something more useful.
708 if self
.is_trivial():
709 return MatrixSpace(self
.base_ring(), 0)
710 elif self
._matrix
_basis
is None or len(self
._matrix
_basis
) == 0:
711 return MatrixSpace(self
.base_ring(), self
.dimension(), 1)
713 return self
._matrix
_basis
[0].matrix_space()
719 Return the unit element of this algebra.
723 sage: from mjo.eja.eja_algebra import (HadamardEJA,
728 sage: J = HadamardEJA(5)
730 e0 + e1 + e2 + e3 + e4
734 The identity element acts like the identity::
736 sage: set_random_seed()
737 sage: J = random_eja()
738 sage: x = J.random_element()
739 sage: J.one()*x == x and x*J.one() == x
742 The matrix of the unit element's operator is the identity::
744 sage: set_random_seed()
745 sage: J = random_eja()
746 sage: actual = J.one().operator().matrix()
747 sage: expected = matrix.identity(J.base_ring(), J.dimension())
748 sage: actual == expected
751 Ensure that the cached unit element (often precomputed by
752 hand) agrees with the computed one::
754 sage: set_random_seed()
755 sage: J = random_eja()
756 sage: cached = J.one()
757 sage: J.one.clear_cache()
758 sage: J.one() == cached
762 # We can brute-force compute the matrices of the operators
763 # that correspond to the basis elements of this algebra.
764 # If some linear combination of those basis elements is the
765 # algebra identity, then the same linear combination of
766 # their matrices has to be the identity matrix.
768 # Of course, matrices aren't vectors in sage, so we have to
769 # appeal to the "long vectors" isometry.
770 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
772 # Now we use basic linear algebra to find the coefficients,
773 # of the matrices-as-vectors-linear-combination, which should
774 # work for the original algebra basis too.
775 A
= matrix(self
.base_ring(), oper_vecs
)
777 # We used the isometry on the left-hand side already, but we
778 # still need to do it for the right-hand side. Recall that we
779 # wanted something that summed to the identity matrix.
780 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
782 # Now if there's an identity element in the algebra, this
783 # should work. We solve on the left to avoid having to
784 # transpose the matrix "A".
785 return self
.from_vector(A
.solve_left(b
))
788 def peirce_decomposition(self
, c
):
790 The Peirce decomposition of this algebra relative to the
793 In the future, this can be extended to a complete system of
794 orthogonal idempotents.
798 - ``c`` -- an idempotent of this algebra.
802 A triple (J0, J5, J1) containing two subalgebras and one subspace
805 - ``J0`` -- the algebra on the eigenspace of ``c.operator()``
806 corresponding to the eigenvalue zero.
808 - ``J5`` -- the eigenspace (NOT a subalgebra) of ``c.operator()``
809 corresponding to the eigenvalue one-half.
811 - ``J1`` -- the algebra on the eigenspace of ``c.operator()``
812 corresponding to the eigenvalue one.
814 These are the only possible eigenspaces for that operator, and this
815 algebra is a direct sum of them. The spaces ``J0`` and ``J1`` are
816 orthogonal, and are subalgebras of this algebra with the appropriate
821 sage: from mjo.eja.eja_algebra import random_eja, RealSymmetricEJA
825 The canonical example comes from the symmetric matrices, which
826 decompose into diagonal and off-diagonal parts::
828 sage: J = RealSymmetricEJA(3)
829 sage: C = matrix(QQ, [ [1,0,0],
833 sage: J0,J5,J1 = J.peirce_decomposition(c)
835 Euclidean Jordan algebra of dimension 1...
837 Vector space of degree 6 and dimension 2...
839 Euclidean Jordan algebra of dimension 3...
840 sage: J0.one().to_matrix()
844 sage: orig_df = AA.options.display_format
845 sage: AA.options.display_format = 'radical'
846 sage: J.from_vector(J5.basis()[0]).to_matrix()
850 sage: J.from_vector(J5.basis()[1]).to_matrix()
854 sage: AA.options.display_format = orig_df
855 sage: J1.one().to_matrix()
862 Every algebra decomposes trivially with respect to its identity
865 sage: set_random_seed()
866 sage: J = random_eja()
867 sage: J0,J5,J1 = J.peirce_decomposition(J.one())
868 sage: J0.dimension() == 0 and J5.dimension() == 0
870 sage: J1.superalgebra() == J and J1.dimension() == J.dimension()
873 The decomposition is into eigenspaces, and its components are
874 therefore necessarily orthogonal. Moreover, the identity
875 elements in the two subalgebras are the projections onto their
876 respective subspaces of the superalgebra's identity element::
878 sage: set_random_seed()
879 sage: J = random_eja()
880 sage: x = J.random_element()
881 sage: if not J.is_trivial():
882 ....: while x.is_nilpotent():
883 ....: x = J.random_element()
884 sage: c = x.subalgebra_idempotent()
885 sage: J0,J5,J1 = J.peirce_decomposition(c)
887 sage: for (w,y,z) in zip(J0.basis(), J5.basis(), J1.basis()):
888 ....: w = w.superalgebra_element()
889 ....: y = J.from_vector(y)
890 ....: z = z.superalgebra_element()
891 ....: ipsum += w.inner_product(y).abs()
892 ....: ipsum += w.inner_product(z).abs()
893 ....: ipsum += y.inner_product(z).abs()
896 sage: J1(c) == J1.one()
898 sage: J0(J.one() - c) == J0.one()
902 if not c
.is_idempotent():
903 raise ValueError("element is not idempotent: %s" % c
)
905 from mjo
.eja
.eja_subalgebra
import FiniteDimensionalEuclideanJordanSubalgebra
907 # Default these to what they should be if they turn out to be
908 # trivial, because eigenspaces_left() won't return eigenvalues
909 # corresponding to trivial spaces (e.g. it returns only the
910 # eigenspace corresponding to lambda=1 if you take the
911 # decomposition relative to the identity element).
912 trivial
= FiniteDimensionalEuclideanJordanSubalgebra(self
, ())
913 J0
= trivial
# eigenvalue zero
914 J5
= VectorSpace(self
.base_ring(), 0) # eigenvalue one-half
915 J1
= trivial
# eigenvalue one
917 for (eigval
, eigspace
) in c
.operator().matrix().right_eigenspaces():
918 if eigval
== ~
(self
.base_ring()(2)):
921 gens
= tuple( self
.from_vector(b
) for b
in eigspace
.basis() )
922 subalg
= FiniteDimensionalEuclideanJordanSubalgebra(self
,
930 raise ValueError("unexpected eigenvalue: %s" % eigval
)
935 def random_element(self
, thorough
=False):
937 Return a random element of this algebra.
939 Our algebra superclass method only returns a linear
940 combination of at most two basis elements. We instead
941 want the vector space "random element" method that
942 returns a more diverse selection.
946 - ``thorough`` -- (boolean; default False) whether or not we
947 should generate irrational coefficients for the random
948 element when our base ring is irrational; this slows the
949 algebra operations to a crawl, but any truly random method
953 # For a general base ring... maybe we can trust this to do the
954 # right thing? Unlikely, but.
955 V
= self
.vector_space()
956 v
= V
.random_element()
958 if self
.base_ring() is AA
:
959 # The "random element" method of the algebraic reals is
960 # stupid at the moment, and only returns integers between
961 # -2 and 2, inclusive:
963 # https://trac.sagemath.org/ticket/30875
965 # Instead, we implement our own "random vector" method,
966 # and then coerce that into the algebra. We use the vector
967 # space degree here instead of the dimension because a
968 # subalgebra could (for example) be spanned by only two
969 # vectors, each with five coordinates. We need to
970 # generate all five coordinates.
972 v
*= QQbar
.random_element().real()
974 v
*= QQ
.random_element()
976 return self
.from_vector(V
.coordinate_vector(v
))
978 def random_elements(self
, count
, thorough
=False):
980 Return ``count`` random elements as a tuple.
984 - ``thorough`` -- (boolean; default False) whether or not we
985 should generate irrational coefficients for the random
986 elements when our base ring is irrational; this slows the
987 algebra operations to a crawl, but any truly random method
992 sage: from mjo.eja.eja_algebra import JordanSpinEJA
996 sage: J = JordanSpinEJA(3)
997 sage: x,y,z = J.random_elements(3)
998 sage: all( [ x in J, y in J, z in J ])
1000 sage: len( J.random_elements(10) ) == 10
1004 return tuple( self
.random_element(thorough
)
1005 for idx
in range(count
) )
1009 def _charpoly_coefficients(self
):
1011 The `r` polynomial coefficients of the "characteristic polynomial
1014 n
= self
.dimension()
1015 R
= self
.coordinate_polynomial_ring()
1017 F
= R
.fraction_field()
1020 # From a result in my book, these are the entries of the
1021 # basis representation of L_x.
1022 return sum( vars[k
]*self
.monomial(k
).operator().matrix()[i
,j
]
1025 L_x
= matrix(F
, n
, n
, L_x_i_j
)
1028 if self
.rank
.is_in_cache():
1030 # There's no need to pad the system with redundant
1031 # columns if we *know* they'll be redundant.
1034 # Compute an extra power in case the rank is equal to
1035 # the dimension (otherwise, we would stop at x^(r-1)).
1036 x_powers
= [ (L_x
**k
)*self
.one().to_vector()
1037 for k
in range(n
+1) ]
1038 A
= matrix
.column(F
, x_powers
[:n
])
1039 AE
= A
.extended_echelon_form()
1046 # The theory says that only the first "r" coefficients are
1047 # nonzero, and they actually live in the original polynomial
1048 # ring and not the fraction field. We negate them because
1049 # in the actual characteristic polynomial, they get moved
1050 # to the other side where x^r lives.
1051 return -A_rref
.solve_right(E
*b
).change_ring(R
)[:r
]
1056 Return the rank of this EJA.
1058 This is a cached method because we know the rank a priori for
1059 all of the algebras we can construct. Thus we can avoid the
1060 expensive ``_charpoly_coefficients()`` call unless we truly
1061 need to compute the whole characteristic polynomial.
1065 sage: from mjo.eja.eja_algebra import (HadamardEJA,
1066 ....: JordanSpinEJA,
1067 ....: RealSymmetricEJA,
1068 ....: ComplexHermitianEJA,
1069 ....: QuaternionHermitianEJA,
1074 The rank of the Jordan spin algebra is always two::
1076 sage: JordanSpinEJA(2).rank()
1078 sage: JordanSpinEJA(3).rank()
1080 sage: JordanSpinEJA(4).rank()
1083 The rank of the `n`-by-`n` Hermitian real, complex, or
1084 quaternion matrices is `n`::
1086 sage: RealSymmetricEJA(4).rank()
1088 sage: ComplexHermitianEJA(3).rank()
1090 sage: QuaternionHermitianEJA(2).rank()
1095 Ensure that every EJA that we know how to construct has a
1096 positive integer rank, unless the algebra is trivial in
1097 which case its rank will be zero::
1099 sage: set_random_seed()
1100 sage: J = random_eja()
1104 sage: r > 0 or (r == 0 and J.is_trivial())
1107 Ensure that computing the rank actually works, since the ranks
1108 of all simple algebras are known and will be cached by default::
1110 sage: set_random_seed() # long time
1111 sage: J = random_eja() # long time
1112 sage: caches = J.rank() # long time
1113 sage: J.rank.clear_cache() # long time
1114 sage: J.rank() == cached # long time
1118 return len(self
._charpoly
_coefficients
())
1121 def vector_space(self
):
1123 Return the vector space that underlies this algebra.
1127 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1131 sage: J = RealSymmetricEJA(2)
1132 sage: J.vector_space()
1133 Vector space of dimension 3 over...
1136 return self
.zero().to_vector().parent().ambient_vector_space()
1139 Element
= FiniteDimensionalEuclideanJordanAlgebraElement
1141 class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra
):
1143 New class for algebras whose supplied basis elements have all rational entries.
1147 sage: from mjo.eja.eja_algebra import BilinearFormEJA
1151 The supplied basis is orthonormalized by default::
1153 sage: B = matrix(QQ, [[1, 0, 0], [0, 25, -32], [0, -32, 41]])
1154 sage: J = BilinearFormEJA(B)
1155 sage: J.matrix_basis()
1168 orthonormalize
=True,
1175 # Abuse the check_field parameter to check that the entries of
1176 # out basis (in ambient coordinates) are in the field QQ.
1177 if not all( all(b_i
in QQ
for b_i
in b
.list()) for b
in basis
):
1178 raise TypeError("basis not rational")
1181 vector_basis
= basis
1183 from sage
.structure
.element
import is_Matrix
1184 basis_is_matrices
= False
1188 if is_Matrix(basis
[0]):
1189 basis_is_matrices
= True
1190 from mjo
.eja
.eja_utils
import _vec2mat
1191 vector_basis
= tuple( map(_mat2vec
,basis
) )
1192 degree
= basis
[0].nrows()**2
1194 degree
= basis
[0].degree()
1196 V
= VectorSpace(field
, degree
)
1198 # If we were asked to orthonormalize, and if the orthonormal
1199 # basis is different from the given one, then we also want to
1200 # compute multiplication and inner-product tables for the
1201 # deorthonormalized basis. These can be used later to
1202 # construct a deorthonormalized copy of this algebra over QQ
1203 # in which several operations are much faster.
1204 self
._rational
_algebra
= None
1207 if self
.base_ring() is not QQ
:
1208 # There's no point in constructing the extra algebra if this
1209 # one is already rational. If the original basis is rational
1210 # but normalization would make it irrational, then this whole
1211 # constructor will just fail anyway as it tries to stick an
1212 # irrational number into a rational algebra.
1214 # Note: the same Jordan and inner-products work here,
1215 # because they are necessarily defined with respect to
1216 # ambient coordinates and not any particular basis.
1217 self
._rational
_algebra
= RationalBasisEuclideanJordanAlgebra(
1222 orthonormalize
=False,
1228 # Compute the deorthonormalized tables before we orthonormalize
1229 # the given basis. The "check" parameter here guarantees that
1230 # the basis is linearly-independent.
1231 W
= V
.span_of_basis( vector_basis
, check
=check_axioms
)
1233 # Note: the Jordan and inner-products are defined in terms
1234 # of the ambient basis. It's important that their arguments
1235 # are in ambient coordinates as well.
1237 for j
in range(i
+1):
1238 # given basis w.r.t. ambient coords
1239 q_i
= vector_basis
[i
]
1240 q_j
= vector_basis
[j
]
1242 if basis_is_matrices
:
1246 elt
= jordan_product(q_i
, q_j
)
1247 ip
= inner_product(q_i
, q_j
)
1249 if basis_is_matrices
:
1250 # do another mat2vec because the multiplication
1251 # table is in terms of vectors
1254 # We overwrite the name "vector_basis" in a second, but never modify it
1255 # in place, to this effectively makes a copy of it.
1256 deortho_vector_basis
= vector_basis
1257 self
._deortho
_matrix
= None
1260 from mjo
.eja
.eja_utils
import gram_schmidt
1261 if basis_is_matrices
:
1262 vector_ip
= lambda x
,y
: inner_product(_vec2mat(x
), _vec2mat(y
))
1263 vector_basis
= gram_schmidt(vector_basis
, vector_ip
)
1265 vector_basis
= gram_schmidt(vector_basis
, inner_product
)
1267 # Normalize the "matrix" basis, too!
1268 basis
= vector_basis
1270 if basis_is_matrices
:
1271 basis
= tuple( map(_vec2mat
,basis
) )
1273 W
= V
.span_of_basis( vector_basis
, check
=check_axioms
)
1275 # Now "W" is the vector space of our algebra coordinates. The
1276 # variables "X1", "X2",... refer to the entries of vectors in
1277 # W. Thus to convert back and forth between the orthonormal
1278 # coordinates and the given ones, we need to stick the original
1280 U
= V
.span_of_basis( deortho_vector_basis
, check
=check_axioms
)
1281 self
._deortho
_matrix
= matrix( U
.coordinate_vector(q
)
1282 for q
in vector_basis
)
1284 # If the superclass constructor is going to verify the
1285 # symmetry of this table, it has better at least be
1288 mult_table
= [ [0 for j
in range(n
)] for i
in range(n
) ]
1289 ip_table
= [ [0 for j
in range(n
)] for i
in range(n
) ]
1291 mult_table
= [ [0 for j
in range(i
+1)] for i
in range(n
) ]
1292 ip_table
= [ [0 for j
in range(i
+1)] for i
in range(n
) ]
1294 # Note: the Jordan and inner-products are defined in terms
1295 # of the ambient basis. It's important that their arguments
1296 # are in ambient coordinates as well.
1298 for j
in range(i
+1):
1299 # ortho basis w.r.t. ambient coords
1300 q_i
= vector_basis
[i
]
1301 q_j
= vector_basis
[j
]
1303 if basis_is_matrices
:
1307 elt
= jordan_product(q_i
, q_j
)
1308 ip
= inner_product(q_i
, q_j
)
1310 if basis_is_matrices
:
1311 # do another mat2vec because the multiplication
1312 # table is in terms of vectors
1315 elt
= W
.coordinate_vector(elt
)
1316 mult_table
[i
][j
] = elt
1319 # The tables are square if we're verifying that they
1321 mult_table
[j
][i
] = elt
1324 if basis_is_matrices
:
1328 basis
= tuple( x
.column() for x
in basis
)
1330 super().__init
__(field
,
1335 basis
, # matrix basis
1340 def _charpoly_coefficients(self
):
1344 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
1345 ....: JordanSpinEJA)
1349 The base ring of the resulting polynomial coefficients is what
1350 it should be, and not the rationals (unless the algebra was
1351 already over the rationals)::
1353 sage: J = JordanSpinEJA(3)
1354 sage: J._charpoly_coefficients()
1355 (X1^2 - X2^2 - X3^2, -2*X1)
1356 sage: a0 = J._charpoly_coefficients()[0]
1358 Algebraic Real Field
1359 sage: a0.base_ring()
1360 Algebraic Real Field
1363 if self
.base_ring() is QQ
or self
._rational
_algebra
is None:
1364 # There's no need to construct *another* algebra over the
1365 # rationals if this one is already over the
1366 # rationals. Likewise, if we never orthonormalized our
1367 # basis, we might as well just use the given one.
1368 superclass
= super(RationalBasisEuclideanJordanAlgebra
, self
)
1369 return superclass
._charpoly
_coefficients
()
1371 # Do the computation over the rationals. The answer will be
1372 # the same, because all we've done is a change of basis.
1373 # Then, change back from QQ to our real base ring
1374 a
= ( a_i
.change_ring(self
.base_ring())
1375 for a_i
in self
._rational
_algebra
._charpoly
_coefficients
() )
1377 # Now convert the coordinate variables back to the
1378 # deorthonormalized ones.
1379 R
= self
.coordinate_polynomial_ring()
1380 from sage
.modules
.free_module_element
import vector
1381 X
= vector(R
, R
.gens())
1382 BX
= self
._deortho
_matrix
*X
1384 subs_dict
= { X[i]: BX[i] for i in range(len(X)) }
1385 return tuple( a_i
.subs(subs_dict
) for a_i
in a
)
1387 class ConcreteEuclideanJordanAlgebra(RationalBasisEuclideanJordanAlgebra
):
1389 A class for the Euclidean Jordan algebras that we know by name.
1391 These are the Jordan algebras whose basis, multiplication table,
1392 rank, and so on are known a priori. More to the point, they are
1393 the Euclidean Jordan algebras for which we are able to conjure up
1394 a "random instance."
1398 sage: from mjo.eja.eja_algebra import ConcreteEuclideanJordanAlgebra
1402 Our basis is normalized with respect to the algebra's inner
1403 product, unless we specify otherwise::
1405 sage: set_random_seed()
1406 sage: J = ConcreteEuclideanJordanAlgebra.random_instance()
1407 sage: all( b.norm() == 1 for b in J.gens() )
1410 Since our basis is orthonormal with respect to the algebra's inner
1411 product, and since we know that this algebra is an EJA, any
1412 left-multiplication operator's matrix will be symmetric because
1413 natural->EJA basis representation is an isometry and within the
1414 EJA the operator is self-adjoint by the Jordan axiom::
1416 sage: set_random_seed()
1417 sage: J = ConcreteEuclideanJordanAlgebra.random_instance()
1418 sage: x = J.random_element()
1419 sage: x.operator().is_self_adjoint()
1424 def _max_random_instance_size():
1426 Return an integer "size" that is an upper bound on the size of
1427 this algebra when it is used in a random test
1428 case. Unfortunately, the term "size" is ambiguous -- when
1429 dealing with `R^n` under either the Hadamard or Jordan spin
1430 product, the "size" refers to the dimension `n`. When dealing
1431 with a matrix algebra (real symmetric or complex/quaternion
1432 Hermitian), it refers to the size of the matrix, which is far
1433 less than the dimension of the underlying vector space.
1435 This method must be implemented in each subclass.
1437 raise NotImplementedError
1440 def random_instance(cls
, *args
, **kwargs
):
1442 Return a random instance of this type of algebra.
1444 This method should be implemented in each subclass.
1446 from sage
.misc
.prandom
import choice
1447 eja_class
= choice(cls
.__subclasses
__())
1449 # These all bubble up to the RationalBasisEuclideanJordanAlgebra
1450 # superclass constructor, so any (kw)args valid there are also
1452 return eja_class
.random_instance(*args
, **kwargs
)
1455 class MatrixEuclideanJordanAlgebra
:
1459 Embed the matrix ``M`` into a space of real matrices.
1461 The matrix ``M`` can have entries in any field at the moment:
1462 the real numbers, complex numbers, or quaternions. And although
1463 they are not a field, we can probably support octonions at some
1464 point, too. This function returns a real matrix that "acts like"
1465 the original with respect to matrix multiplication; i.e.
1467 real_embed(M*N) = real_embed(M)*real_embed(N)
1470 raise NotImplementedError
1474 def real_unembed(M
):
1476 The inverse of :meth:`real_embed`.
1478 raise NotImplementedError
1481 def jordan_product(X
,Y
):
1482 return (X
*Y
+ Y
*X
)/2
1485 def trace_inner_product(cls
,X
,Y
):
1486 Xu
= cls
.real_unembed(X
)
1487 Yu
= cls
.real_unembed(Y
)
1488 tr
= (Xu
*Yu
).trace()
1491 # Works in QQ, AA, RDF, et cetera.
1493 except AttributeError:
1494 # A quaternion doesn't have a real() method, but does
1495 # have coefficient_tuple() method that returns the
1496 # coefficients of 1, i, j, and k -- in that order.
1497 return tr
.coefficient_tuple()[0]
1500 class RealMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1504 The identity function, for embedding real matrices into real
1510 def real_unembed(M
):
1512 The identity function, for unembedding real matrices from real
1518 class RealSymmetricEJA(ConcreteEuclideanJordanAlgebra
,
1519 RealMatrixEuclideanJordanAlgebra
):
1521 The rank-n simple EJA consisting of real symmetric n-by-n
1522 matrices, the usual symmetric Jordan product, and the trace inner
1523 product. It has dimension `(n^2 + n)/2` over the reals.
1527 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1531 sage: J = RealSymmetricEJA(2)
1532 sage: e0, e1, e2 = J.gens()
1540 In theory, our "field" can be any subfield of the reals::
1542 sage: RealSymmetricEJA(2, RDF)
1543 Euclidean Jordan algebra of dimension 3 over Real Double Field
1544 sage: RealSymmetricEJA(2, RR)
1545 Euclidean Jordan algebra of dimension 3 over Real Field with
1546 53 bits of precision
1550 The dimension of this algebra is `(n^2 + n) / 2`::
1552 sage: set_random_seed()
1553 sage: n_max = RealSymmetricEJA._max_random_instance_size()
1554 sage: n = ZZ.random_element(1, n_max)
1555 sage: J = RealSymmetricEJA(n)
1556 sage: J.dimension() == (n^2 + n)/2
1559 The Jordan multiplication is what we think it is::
1561 sage: set_random_seed()
1562 sage: J = RealSymmetricEJA.random_instance()
1563 sage: x,y = J.random_elements(2)
1564 sage: actual = (x*y).to_matrix()
1565 sage: X = x.to_matrix()
1566 sage: Y = y.to_matrix()
1567 sage: expected = (X*Y + Y*X)/2
1568 sage: actual == expected
1570 sage: J(expected) == x*y
1573 We can change the generator prefix::
1575 sage: RealSymmetricEJA(3, prefix='q').gens()
1576 (q0, q1, q2, q3, q4, q5)
1578 We can construct the (trivial) algebra of rank zero::
1580 sage: RealSymmetricEJA(0)
1581 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1585 def _denormalized_basis(cls
, n
, field
):
1587 Return a basis for the space of real symmetric n-by-n matrices.
1591 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1595 sage: set_random_seed()
1596 sage: n = ZZ.random_element(1,5)
1597 sage: B = RealSymmetricEJA._denormalized_basis(n,QQ)
1598 sage: all( M.is_symmetric() for M in B)
1602 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1606 for j
in range(i
+1):
1607 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
1611 Sij
= Eij
+ Eij
.transpose()
1617 def _max_random_instance_size():
1618 return 4 # Dimension 10
1621 def random_instance(cls
, field
=AA
, **kwargs
):
1623 Return a random instance of this type of algebra.
1625 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1626 return cls(n
, field
, **kwargs
)
1628 def __init__(self
, n
, field
=AA
, **kwargs
):
1629 basis
= self
._denormalized
_basis
(n
, field
)
1630 super(RealSymmetricEJA
, self
).__init
__(field
,
1632 self
.jordan_product
,
1633 self
.trace_inner_product
,
1635 self
.rank
.set_cache(n
)
1636 self
.one
.set_cache(self(matrix
.identity(field
,n
)))
1639 class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1643 Embed the n-by-n complex matrix ``M`` into the space of real
1644 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1645 bi` to the block matrix ``[[a,b],[-b,a]]``.
1649 sage: from mjo.eja.eja_algebra import \
1650 ....: ComplexMatrixEuclideanJordanAlgebra
1654 sage: F = QuadraticField(-1, 'I')
1655 sage: x1 = F(4 - 2*i)
1656 sage: x2 = F(1 + 2*i)
1659 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1660 sage: ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1669 Embedding is a homomorphism (isomorphism, in fact)::
1671 sage: set_random_seed()
1672 sage: n = ZZ.random_element(3)
1673 sage: F = QuadraticField(-1, 'I')
1674 sage: X = random_matrix(F, n)
1675 sage: Y = random_matrix(F, n)
1676 sage: Xe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X)
1677 sage: Ye = ComplexMatrixEuclideanJordanAlgebra.real_embed(Y)
1678 sage: XYe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1685 raise ValueError("the matrix 'M' must be square")
1687 # We don't need any adjoined elements...
1688 field
= M
.base_ring().base_ring()
1692 a
= z
.list()[0] # real part, I guess
1693 b
= z
.list()[1] # imag part, I guess
1694 blocks
.append(matrix(field
, 2, [[a
,b
],[-b
,a
]]))
1696 return matrix
.block(field
, n
, blocks
)
1700 def real_unembed(M
):
1702 The inverse of _embed_complex_matrix().
1706 sage: from mjo.eja.eja_algebra import \
1707 ....: ComplexMatrixEuclideanJordanAlgebra
1711 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1712 ....: [-2, 1, -4, 3],
1713 ....: [ 9, 10, 11, 12],
1714 ....: [-10, 9, -12, 11] ])
1715 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(A)
1717 [ 10*I + 9 12*I + 11]
1721 Unembedding is the inverse of embedding::
1723 sage: set_random_seed()
1724 sage: F = QuadraticField(-1, 'I')
1725 sage: M = random_matrix(F, 3)
1726 sage: Me = ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1727 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1733 raise ValueError("the matrix 'M' must be square")
1734 if not n
.mod(2).is_zero():
1735 raise ValueError("the matrix 'M' must be a complex embedding")
1737 # If "M" was normalized, its base ring might have roots
1738 # adjoined and they can stick around after unembedding.
1739 field
= M
.base_ring()
1740 R
= PolynomialRing(field
, 'z')
1743 # Sage doesn't know how to embed AA into QQbar, i.e. how
1744 # to adjoin sqrt(-1) to AA.
1747 F
= field
.extension(z
**2 + 1, 'I', embedding
=CLF(-1).sqrt())
1750 # Go top-left to bottom-right (reading order), converting every
1751 # 2-by-2 block we see to a single complex element.
1753 for k
in range(n
/2):
1754 for j
in range(n
/2):
1755 submat
= M
[2*k
:2*k
+2,2*j
:2*j
+2]
1756 if submat
[0,0] != submat
[1,1]:
1757 raise ValueError('bad on-diagonal submatrix')
1758 if submat
[0,1] != -submat
[1,0]:
1759 raise ValueError('bad off-diagonal submatrix')
1760 z
= submat
[0,0] + submat
[0,1]*i
1763 return matrix(F
, n
/2, elements
)
1767 def trace_inner_product(cls
,X
,Y
):
1769 Compute a matrix inner product in this algebra directly from
1774 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1778 This gives the same answer as the slow, default method implemented
1779 in :class:`MatrixEuclideanJordanAlgebra`::
1781 sage: set_random_seed()
1782 sage: J = ComplexHermitianEJA.random_instance()
1783 sage: x,y = J.random_elements(2)
1784 sage: Xe = x.to_matrix()
1785 sage: Ye = y.to_matrix()
1786 sage: X = ComplexHermitianEJA.real_unembed(Xe)
1787 sage: Y = ComplexHermitianEJA.real_unembed(Ye)
1788 sage: expected = (X*Y).trace().real()
1789 sage: actual = ComplexHermitianEJA.trace_inner_product(Xe,Ye)
1790 sage: actual == expected
1794 return RealMatrixEuclideanJordanAlgebra
.trace_inner_product(X
,Y
)/2
1797 class ComplexHermitianEJA(ConcreteEuclideanJordanAlgebra
,
1798 ComplexMatrixEuclideanJordanAlgebra
):
1800 The rank-n simple EJA consisting of complex Hermitian n-by-n
1801 matrices over the real numbers, the usual symmetric Jordan product,
1802 and the real-part-of-trace inner product. It has dimension `n^2` over
1807 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1811 In theory, our "field" can be any subfield of the reals::
1813 sage: ComplexHermitianEJA(2, RDF)
1814 Euclidean Jordan algebra of dimension 4 over Real Double Field
1815 sage: ComplexHermitianEJA(2, RR)
1816 Euclidean Jordan algebra of dimension 4 over Real Field with
1817 53 bits of precision
1821 The dimension of this algebra is `n^2`::
1823 sage: set_random_seed()
1824 sage: n_max = ComplexHermitianEJA._max_random_instance_size()
1825 sage: n = ZZ.random_element(1, n_max)
1826 sage: J = ComplexHermitianEJA(n)
1827 sage: J.dimension() == n^2
1830 The Jordan multiplication is what we think it is::
1832 sage: set_random_seed()
1833 sage: J = ComplexHermitianEJA.random_instance()
1834 sage: x,y = J.random_elements(2)
1835 sage: actual = (x*y).to_matrix()
1836 sage: X = x.to_matrix()
1837 sage: Y = y.to_matrix()
1838 sage: expected = (X*Y + Y*X)/2
1839 sage: actual == expected
1841 sage: J(expected) == x*y
1844 We can change the generator prefix::
1846 sage: ComplexHermitianEJA(2, prefix='z').gens()
1849 We can construct the (trivial) algebra of rank zero::
1851 sage: ComplexHermitianEJA(0)
1852 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1857 def _denormalized_basis(cls
, n
, field
):
1859 Returns a basis for the space of complex Hermitian n-by-n matrices.
1861 Why do we embed these? Basically, because all of numerical linear
1862 algebra assumes that you're working with vectors consisting of `n`
1863 entries from a field and scalars from the same field. There's no way
1864 to tell SageMath that (for example) the vectors contain complex
1865 numbers, while the scalar field is real.
1869 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1873 sage: set_random_seed()
1874 sage: n = ZZ.random_element(1,5)
1875 sage: field = QuadraticField(2, 'sqrt2')
1876 sage: B = ComplexHermitianEJA._denormalized_basis(n, field)
1877 sage: all( M.is_symmetric() for M in B)
1881 R
= PolynomialRing(field
, 'z')
1883 F
= field
.extension(z
**2 + 1, 'I')
1886 # This is like the symmetric case, but we need to be careful:
1888 # * We want conjugate-symmetry, not just symmetry.
1889 # * The diagonal will (as a result) be real.
1893 for j
in range(i
+1):
1894 Eij
= matrix(F
, n
, lambda k
,l
: k
==i
and l
==j
)
1896 Sij
= cls
.real_embed(Eij
)
1899 # The second one has a minus because it's conjugated.
1900 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
1902 Sij_imag
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
1905 # Since we embedded these, we can drop back to the "field" that we
1906 # started with instead of the complex extension "F".
1907 return tuple( s
.change_ring(field
) for s
in S
)
1910 def __init__(self
, n
, field
=AA
, **kwargs
):
1911 basis
= self
._denormalized
_basis
(n
,field
)
1912 super(ComplexHermitianEJA
, self
).__init
__(field
,
1914 self
.jordan_product
,
1915 self
.trace_inner_product
,
1917 self
.rank
.set_cache(n
)
1918 # TODO: pre-cache the identity!
1921 def _max_random_instance_size():
1922 return 3 # Dimension 9
1925 def random_instance(cls
, field
=AA
, **kwargs
):
1927 Return a random instance of this type of algebra.
1929 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1930 return cls(n
, field
, **kwargs
)
1932 class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1936 Embed the n-by-n quaternion matrix ``M`` into the space of real
1937 matrices of size 4n-by-4n by first sending each quaternion entry `z
1938 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
1939 c+di],[-c + di, a-bi]]`, and then embedding those into a real
1944 sage: from mjo.eja.eja_algebra import \
1945 ....: QuaternionMatrixEuclideanJordanAlgebra
1949 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1950 sage: i,j,k = Q.gens()
1951 sage: x = 1 + 2*i + 3*j + 4*k
1952 sage: M = matrix(Q, 1, [[x]])
1953 sage: QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1959 Embedding is a homomorphism (isomorphism, in fact)::
1961 sage: set_random_seed()
1962 sage: n = ZZ.random_element(2)
1963 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1964 sage: X = random_matrix(Q, n)
1965 sage: Y = random_matrix(Q, n)
1966 sage: Xe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X)
1967 sage: Ye = QuaternionMatrixEuclideanJordanAlgebra.real_embed(Y)
1968 sage: XYe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1973 quaternions
= M
.base_ring()
1976 raise ValueError("the matrix 'M' must be square")
1978 F
= QuadraticField(-1, 'I')
1983 t
= z
.coefficient_tuple()
1988 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
1989 [-c
+ d
*i
, a
- b
*i
]])
1990 realM
= ComplexMatrixEuclideanJordanAlgebra
.real_embed(cplxM
)
1991 blocks
.append(realM
)
1993 # We should have real entries by now, so use the realest field
1994 # we've got for the return value.
1995 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
2000 def real_unembed(M
):
2002 The inverse of _embed_quaternion_matrix().
2006 sage: from mjo.eja.eja_algebra import \
2007 ....: QuaternionMatrixEuclideanJordanAlgebra
2011 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
2012 ....: [-2, 1, -4, 3],
2013 ....: [-3, 4, 1, -2],
2014 ....: [-4, -3, 2, 1]])
2015 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(M)
2016 [1 + 2*i + 3*j + 4*k]
2020 Unembedding is the inverse of embedding::
2022 sage: set_random_seed()
2023 sage: Q = QuaternionAlgebra(QQ, -1, -1)
2024 sage: M = random_matrix(Q, 3)
2025 sage: Me = QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
2026 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
2032 raise ValueError("the matrix 'M' must be square")
2033 if not n
.mod(4).is_zero():
2034 raise ValueError("the matrix 'M' must be a quaternion embedding")
2036 # Use the base ring of the matrix to ensure that its entries can be
2037 # multiplied by elements of the quaternion algebra.
2038 field
= M
.base_ring()
2039 Q
= QuaternionAlgebra(field
,-1,-1)
2042 # Go top-left to bottom-right (reading order), converting every
2043 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
2046 for l
in range(n
/4):
2047 for m
in range(n
/4):
2048 submat
= ComplexMatrixEuclideanJordanAlgebra
.real_unembed(
2049 M
[4*l
:4*l
+4,4*m
:4*m
+4] )
2050 if submat
[0,0] != submat
[1,1].conjugate():
2051 raise ValueError('bad on-diagonal submatrix')
2052 if submat
[0,1] != -submat
[1,0].conjugate():
2053 raise ValueError('bad off-diagonal submatrix')
2054 z
= submat
[0,0].real()
2055 z
+= submat
[0,0].imag()*i
2056 z
+= submat
[0,1].real()*j
2057 z
+= submat
[0,1].imag()*k
2060 return matrix(Q
, n
/4, elements
)
2064 def trace_inner_product(cls
,X
,Y
):
2066 Compute a matrix inner product in this algebra directly from
2071 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2075 This gives the same answer as the slow, default method implemented
2076 in :class:`MatrixEuclideanJordanAlgebra`::
2078 sage: set_random_seed()
2079 sage: J = QuaternionHermitianEJA.random_instance()
2080 sage: x,y = J.random_elements(2)
2081 sage: Xe = x.to_matrix()
2082 sage: Ye = y.to_matrix()
2083 sage: X = QuaternionHermitianEJA.real_unembed(Xe)
2084 sage: Y = QuaternionHermitianEJA.real_unembed(Ye)
2085 sage: expected = (X*Y).trace().coefficient_tuple()[0]
2086 sage: actual = QuaternionHermitianEJA.trace_inner_product(Xe,Ye)
2087 sage: actual == expected
2091 return RealMatrixEuclideanJordanAlgebra
.trace_inner_product(X
,Y
)/4
2094 class QuaternionHermitianEJA(ConcreteEuclideanJordanAlgebra
,
2095 QuaternionMatrixEuclideanJordanAlgebra
):
2097 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
2098 matrices, the usual symmetric Jordan product, and the
2099 real-part-of-trace inner product. It has dimension `2n^2 - n` over
2104 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2108 In theory, our "field" can be any subfield of the reals::
2110 sage: QuaternionHermitianEJA(2, RDF)
2111 Euclidean Jordan algebra of dimension 6 over Real Double Field
2112 sage: QuaternionHermitianEJA(2, RR)
2113 Euclidean Jordan algebra of dimension 6 over Real Field with
2114 53 bits of precision
2118 The dimension of this algebra is `2*n^2 - n`::
2120 sage: set_random_seed()
2121 sage: n_max = QuaternionHermitianEJA._max_random_instance_size()
2122 sage: n = ZZ.random_element(1, n_max)
2123 sage: J = QuaternionHermitianEJA(n)
2124 sage: J.dimension() == 2*(n^2) - n
2127 The Jordan multiplication is what we think it is::
2129 sage: set_random_seed()
2130 sage: J = QuaternionHermitianEJA.random_instance()
2131 sage: x,y = J.random_elements(2)
2132 sage: actual = (x*y).to_matrix()
2133 sage: X = x.to_matrix()
2134 sage: Y = y.to_matrix()
2135 sage: expected = (X*Y + Y*X)/2
2136 sage: actual == expected
2138 sage: J(expected) == x*y
2141 We can change the generator prefix::
2143 sage: QuaternionHermitianEJA(2, prefix='a').gens()
2144 (a0, a1, a2, a3, a4, a5)
2146 We can construct the (trivial) algebra of rank zero::
2148 sage: QuaternionHermitianEJA(0)
2149 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2153 def _denormalized_basis(cls
, n
, field
):
2155 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
2157 Why do we embed these? Basically, because all of numerical
2158 linear algebra assumes that you're working with vectors consisting
2159 of `n` entries from a field and scalars from the same field. There's
2160 no way to tell SageMath that (for example) the vectors contain
2161 complex numbers, while the scalar field is real.
2165 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2169 sage: set_random_seed()
2170 sage: n = ZZ.random_element(1,5)
2171 sage: B = QuaternionHermitianEJA._denormalized_basis(n,QQ)
2172 sage: all( M.is_symmetric() for M in B )
2176 Q
= QuaternionAlgebra(QQ
,-1,-1)
2179 # This is like the symmetric case, but we need to be careful:
2181 # * We want conjugate-symmetry, not just symmetry.
2182 # * The diagonal will (as a result) be real.
2186 for j
in range(i
+1):
2187 Eij
= matrix(Q
, n
, lambda k
,l
: k
==i
and l
==j
)
2189 Sij
= cls
.real_embed(Eij
)
2192 # The second, third, and fourth ones have a minus
2193 # because they're conjugated.
2194 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
2196 Sij_I
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
2198 Sij_J
= cls
.real_embed(J
*Eij
- J
*Eij
.transpose())
2200 Sij_K
= cls
.real_embed(K
*Eij
- K
*Eij
.transpose())
2203 # Since we embedded these, we can drop back to the "field" that we
2204 # started with instead of the quaternion algebra "Q".
2205 return tuple( s
.change_ring(field
) for s
in S
)
2208 def __init__(self
, n
, field
=AA
, **kwargs
):
2209 basis
= self
._denormalized
_basis
(n
,field
)
2210 super(QuaternionHermitianEJA
, self
).__init
__(field
,
2212 self
.jordan_product
,
2213 self
.trace_inner_product
,
2215 self
.rank
.set_cache(n
)
2216 # TODO: cache one()!
2219 def _max_random_instance_size():
2221 The maximum rank of a random QuaternionHermitianEJA.
2223 return 2 # Dimension 6
2226 def random_instance(cls
, field
=AA
, **kwargs
):
2228 Return a random instance of this type of algebra.
2230 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2231 return cls(n
, field
, **kwargs
)
2234 class HadamardEJA(ConcreteEuclideanJordanAlgebra
):
2236 Return the Euclidean Jordan Algebra corresponding to the set
2237 `R^n` under the Hadamard product.
2239 Note: this is nothing more than the Cartesian product of ``n``
2240 copies of the spin algebra. Once Cartesian product algebras
2241 are implemented, this can go.
2245 sage: from mjo.eja.eja_algebra import HadamardEJA
2249 This multiplication table can be verified by hand::
2251 sage: J = HadamardEJA(3)
2252 sage: e0,e1,e2 = J.gens()
2268 We can change the generator prefix::
2270 sage: HadamardEJA(3, prefix='r').gens()
2274 def __init__(self
, n
, field
=AA
, **kwargs
):
2275 V
= VectorSpace(field
, n
)
2278 def jordan_product(x
,y
):
2279 return V([ xi
*yi
for (xi
,yi
) in zip(x
,y
) ])
2280 def inner_product(x
,y
):
2281 return x
.inner_product(y
)
2283 # Don't orthonormalize because our basis is already
2284 # orthonormal with respect to our inner-product. But also
2285 # don't pass check_field=False here, because the user can pass
2287 super(HadamardEJA
, self
).__init
__(field
,
2291 orthonormalize
=False,
2294 self
.rank
.set_cache(n
)
2297 self
.one
.set_cache( self
.zero() )
2299 self
.one
.set_cache( sum(self
.gens()) )
2302 def _max_random_instance_size():
2304 The maximum dimension of a random HadamardEJA.
2309 def random_instance(cls
, field
=AA
, **kwargs
):
2311 Return a random instance of this type of algebra.
2313 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2314 return cls(n
, field
, **kwargs
)
2317 class BilinearFormEJA(ConcreteEuclideanJordanAlgebra
):
2319 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2320 with the half-trace inner product and jordan product ``x*y =
2321 (<Bx,y>,y_bar>, x0*y_bar + y0*x_bar)`` where `B = 1 \times B22` is
2322 a symmetric positive-definite "bilinear form" matrix. Its
2323 dimension is the size of `B`, and it has rank two in dimensions
2324 larger than two. It reduces to the ``JordanSpinEJA`` when `B` is
2325 the identity matrix of order ``n``.
2327 We insist that the one-by-one upper-left identity block of `B` be
2328 passed in as well so that we can be passed a matrix of size zero
2329 to construct a trivial algebra.
2333 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
2334 ....: JordanSpinEJA)
2338 When no bilinear form is specified, the identity matrix is used,
2339 and the resulting algebra is the Jordan spin algebra::
2341 sage: B = matrix.identity(AA,3)
2342 sage: J0 = BilinearFormEJA(B)
2343 sage: J1 = JordanSpinEJA(3)
2344 sage: J0.multiplication_table() == J0.multiplication_table()
2347 An error is raised if the matrix `B` does not correspond to a
2348 positive-definite bilinear form::
2350 sage: B = matrix.random(QQ,2,3)
2351 sage: J = BilinearFormEJA(B)
2352 Traceback (most recent call last):
2354 ValueError: bilinear form is not positive-definite
2355 sage: B = matrix.zero(QQ,3)
2356 sage: J = BilinearFormEJA(B)
2357 Traceback (most recent call last):
2359 ValueError: bilinear form is not positive-definite
2363 We can create a zero-dimensional algebra::
2365 sage: B = matrix.identity(AA,0)
2366 sage: J = BilinearFormEJA(B)
2370 We can check the multiplication condition given in the Jordan, von
2371 Neumann, and Wigner paper (and also discussed on my "On the
2372 symmetry..." paper). Note that this relies heavily on the standard
2373 choice of basis, as does anything utilizing the bilinear form
2374 matrix. We opt not to orthonormalize the basis, because if we
2375 did, we would have to normalize the `s_{i}` in a similar manner::
2377 sage: set_random_seed()
2378 sage: n = ZZ.random_element(5)
2379 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2380 sage: B11 = matrix.identity(QQ,1)
2381 sage: B22 = M.transpose()*M
2382 sage: B = block_matrix(2,2,[ [B11,0 ],
2384 sage: J = BilinearFormEJA(B, orthonormalize=False)
2385 sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
2386 sage: V = J.vector_space()
2387 sage: sis = [ J( V([0] + (M.inverse()*ei).list()).column() )
2388 ....: for ei in eis ]
2389 sage: actual = [ sis[i]*sis[j]
2390 ....: for i in range(n-1)
2391 ....: for j in range(n-1) ]
2392 sage: expected = [ J.one() if i == j else J.zero()
2393 ....: for i in range(n-1)
2394 ....: for j in range(n-1) ]
2395 sage: actual == expected
2398 def __init__(self
, B
, field
=AA
, **kwargs
):
2399 if not B
.is_positive_definite():
2400 raise ValueError("bilinear form is not positive-definite")
2403 V
= VectorSpace(field
, n
)
2405 def inner_product(x
,y
):
2406 return (B
*x
).inner_product(y
)
2408 def jordan_product(x
,y
):
2413 z0
= inner_product(x
,y
)
2414 zbar
= y0
*xbar
+ x0
*ybar
2415 return V([z0
] + zbar
.list())
2417 super(BilinearFormEJA
, self
).__init
__(field
,
2423 # The rank of this algebra is two, unless we're in a
2424 # one-dimensional ambient space (because the rank is bounded
2425 # by the ambient dimension).
2426 self
.rank
.set_cache(min(n
,2))
2429 self
.one
.set_cache( self
.zero() )
2431 self
.one
.set_cache( self
.monomial(0) )
2434 def _max_random_instance_size():
2436 The maximum dimension of a random BilinearFormEJA.
2441 def random_instance(cls
, field
=AA
, **kwargs
):
2443 Return a random instance of this algebra.
2445 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2447 B
= matrix
.identity(field
, n
)
2448 return cls(B
, field
, **kwargs
)
2450 B11
= matrix
.identity(field
,1)
2451 M
= matrix
.random(field
, n
-1)
2452 I
= matrix
.identity(field
, n
-1)
2453 alpha
= field
.zero()
2454 while alpha
.is_zero():
2455 alpha
= field
.random_element().abs()
2456 B22
= M
.transpose()*M
+ alpha
*I
2458 from sage
.matrix
.special
import block_matrix
2459 B
= block_matrix(2,2, [ [B11
, ZZ(0) ],
2462 return cls(B
, field
, **kwargs
)
2465 class JordanSpinEJA(BilinearFormEJA
):
2467 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2468 with the usual inner product and jordan product ``x*y =
2469 (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2474 sage: from mjo.eja.eja_algebra import JordanSpinEJA
2478 This multiplication table can be verified by hand::
2480 sage: J = JordanSpinEJA(4)
2481 sage: e0,e1,e2,e3 = J.gens()
2497 We can change the generator prefix::
2499 sage: JordanSpinEJA(2, prefix='B').gens()
2504 Ensure that we have the usual inner product on `R^n`::
2506 sage: set_random_seed()
2507 sage: J = JordanSpinEJA.random_instance()
2508 sage: x,y = J.random_elements(2)
2509 sage: actual = x.inner_product(y)
2510 sage: expected = x.to_vector().inner_product(y.to_vector())
2511 sage: actual == expected
2515 def __init__(self
, n
, field
=AA
, **kwargs
):
2516 # This is a special case of the BilinearFormEJA with the
2517 # identity matrix as its bilinear form.
2518 B
= matrix
.identity(field
, n
)
2520 # Don't orthonormalize because our basis is already
2521 # orthonormal with respect to our inner-product. But
2522 # also don't pass check_field=False here, because the
2523 # user can pass in a field!
2524 super(JordanSpinEJA
, self
).__init
__(B
,
2526 orthonormalize
=False,
2531 def _max_random_instance_size():
2533 The maximum dimension of a random JordanSpinEJA.
2538 def random_instance(cls
, field
=AA
, **kwargs
):
2540 Return a random instance of this type of algebra.
2542 Needed here to override the implementation for ``BilinearFormEJA``.
2544 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2545 return cls(n
, field
, **kwargs
)
2548 class TrivialEJA(ConcreteEuclideanJordanAlgebra
):
2550 The trivial Euclidean Jordan algebra consisting of only a zero element.
2554 sage: from mjo.eja.eja_algebra import TrivialEJA
2558 sage: J = TrivialEJA()
2565 sage: 7*J.one()*12*J.one()
2567 sage: J.one().inner_product(J.one())
2569 sage: J.one().norm()
2571 sage: J.one().subalgebra_generated_by()
2572 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2577 def __init__(self
, field
=AA
, **kwargs
):
2578 jordan_product
= lambda x
,y
: x
2579 inner_product
= lambda x
,y
: field(0)
2581 super(TrivialEJA
, self
).__init
__(field
,
2586 # The rank is zero using my definition, namely the dimension of the
2587 # largest subalgebra generated by any element.
2588 self
.rank
.set_cache(0)
2589 self
.one
.set_cache( self
.zero() )
2592 def random_instance(cls
, field
=AA
, **kwargs
):
2593 # We don't take a "size" argument so the superclass method is
2594 # inappropriate for us.
2595 return cls(field
, **kwargs
)
2597 class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra
):
2599 The external (orthogonal) direct sum of two other Euclidean Jordan
2600 algebras. Essentially the Cartesian product of its two factors.
2601 Every Euclidean Jordan algebra decomposes into an orthogonal
2602 direct sum of simple Euclidean Jordan algebras, so no generality
2603 is lost by providing only this construction.
2607 sage: from mjo.eja.eja_algebra import (random_eja,
2609 ....: RealSymmetricEJA,
2614 sage: J1 = HadamardEJA(2)
2615 sage: J2 = RealSymmetricEJA(3)
2616 sage: J = DirectSumEJA(J1,J2)
2624 The external direct sum construction is only valid when the two factors
2625 have the same base ring; an error is raised otherwise::
2627 sage: set_random_seed()
2628 sage: J1 = random_eja(AA)
2629 sage: J2 = random_eja(QQ,orthonormalize=False)
2630 sage: J = DirectSumEJA(J1,J2)
2631 Traceback (most recent call last):
2633 ValueError: algebras must share the same base field
2636 def __init__(self
, J1
, J2
, **kwargs
):
2637 if J1
.base_ring() != J2
.base_ring():
2638 raise ValueError("algebras must share the same base field")
2639 field
= J1
.base_ring()
2641 self
._factors
= (J1
, J2
)
2645 V
= VectorSpace(field
, n
)
2646 mult_table
= [ [ V
.zero() for j
in range(i
+1) ]
2649 for j
in range(i
+1):
2650 p
= (J1
.monomial(i
)*J1
.monomial(j
)).to_vector()
2651 mult_table
[i
][j
] = V(p
.list() + [field
.zero()]*n2
)
2654 for j
in range(i
+1):
2655 p
= (J2
.monomial(i
)*J2
.monomial(j
)).to_vector()
2656 mult_table
[n1
+i
][n1
+j
] = V([field
.zero()]*n1
+ p
.list())
2658 # TODO: build the IP table here from the two constituent IP
2659 # matrices (it'll be block diagonal, I think).
2660 ip_table
= [ [ field
.zero() for j
in range(i
+1) ]
2662 super(DirectSumEJA
, self
).__init
__(field
,
2667 self
.rank
.set_cache(J1
.rank() + J2
.rank())
2672 Return the pair of this algebra's factors.
2676 sage: from mjo.eja.eja_algebra import (HadamardEJA,
2677 ....: JordanSpinEJA,
2682 sage: J1 = HadamardEJA(2,QQ)
2683 sage: J2 = JordanSpinEJA(3,QQ)
2684 sage: J = DirectSumEJA(J1,J2)
2686 (Euclidean Jordan algebra of dimension 2 over Rational Field,
2687 Euclidean Jordan algebra of dimension 3 over Rational Field)
2690 return self
._factors
2692 def projections(self
):
2694 Return a pair of projections onto this algebra's factors.
2698 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
2699 ....: ComplexHermitianEJA,
2704 sage: J1 = JordanSpinEJA(2)
2705 sage: J2 = ComplexHermitianEJA(2)
2706 sage: J = DirectSumEJA(J1,J2)
2707 sage: (pi_left, pi_right) = J.projections()
2708 sage: J.one().to_vector()
2710 sage: pi_left(J.one()).to_vector()
2712 sage: pi_right(J.one()).to_vector()
2716 (J1
,J2
) = self
.factors()
2719 V_basis
= self
.vector_space().basis()
2720 # Need to specify the dimensions explicitly so that we don't
2721 # wind up with a zero-by-zero matrix when we want e.g. a
2722 # zero-by-two matrix (important for composing things).
2723 P1
= matrix(self
.base_ring(), m
, m
+n
, V_basis
[:m
])
2724 P2
= matrix(self
.base_ring(), n
, m
+n
, V_basis
[m
:])
2725 pi_left
= FiniteDimensionalEuclideanJordanAlgebraOperator(self
,J1
,P1
)
2726 pi_right
= FiniteDimensionalEuclideanJordanAlgebraOperator(self
,J2
,P2
)
2727 return (pi_left
, pi_right
)
2729 def inclusions(self
):
2731 Return the pair of inclusion maps from our factors into us.
2735 sage: from mjo.eja.eja_algebra import (random_eja,
2736 ....: JordanSpinEJA,
2737 ....: RealSymmetricEJA,
2742 sage: J1 = JordanSpinEJA(3)
2743 sage: J2 = RealSymmetricEJA(2)
2744 sage: J = DirectSumEJA(J1,J2)
2745 sage: (iota_left, iota_right) = J.inclusions()
2746 sage: iota_left(J1.zero()) == J.zero()
2748 sage: iota_right(J2.zero()) == J.zero()
2750 sage: J1.one().to_vector()
2752 sage: iota_left(J1.one()).to_vector()
2754 sage: J2.one().to_vector()
2756 sage: iota_right(J2.one()).to_vector()
2758 sage: J.one().to_vector()
2763 Composing a projection with the corresponding inclusion should
2764 produce the identity map, and mismatching them should produce
2767 sage: set_random_seed()
2768 sage: J1 = random_eja()
2769 sage: J2 = random_eja()
2770 sage: J = DirectSumEJA(J1,J2)
2771 sage: (iota_left, iota_right) = J.inclusions()
2772 sage: (pi_left, pi_right) = J.projections()
2773 sage: pi_left*iota_left == J1.one().operator()
2775 sage: pi_right*iota_right == J2.one().operator()
2777 sage: (pi_left*iota_right).is_zero()
2779 sage: (pi_right*iota_left).is_zero()
2783 (J1
,J2
) = self
.factors()
2786 V_basis
= self
.vector_space().basis()
2787 # Need to specify the dimensions explicitly so that we don't
2788 # wind up with a zero-by-zero matrix when we want e.g. a
2789 # two-by-zero matrix (important for composing things).
2790 I1
= matrix
.column(self
.base_ring(), m
, m
+n
, V_basis
[:m
])
2791 I2
= matrix
.column(self
.base_ring(), n
, m
+n
, V_basis
[m
:])
2792 iota_left
= FiniteDimensionalEuclideanJordanAlgebraOperator(J1
,self
,I1
)
2793 iota_right
= FiniteDimensionalEuclideanJordanAlgebraOperator(J2
,self
,I2
)
2794 return (iota_left
, iota_right
)
2796 def inner_product(self
, x
, y
):
2798 The standard Cartesian inner-product.
2800 We project ``x`` and ``y`` onto our factors, and add up the
2801 inner-products from the subalgebras.
2806 sage: from mjo.eja.eja_algebra import (HadamardEJA,
2807 ....: QuaternionHermitianEJA,
2812 sage: J1 = HadamardEJA(3,QQ)
2813 sage: J2 = QuaternionHermitianEJA(2,QQ,orthonormalize=False)
2814 sage: J = DirectSumEJA(J1,J2)
2819 sage: x1.inner_product(x2)
2821 sage: y1.inner_product(y2)
2823 sage: J.one().inner_product(J.one())
2827 (pi_left
, pi_right
) = self
.projections()
2833 return (x1
.inner_product(y1
) + x2
.inner_product(y2
))
2837 random_eja
= ConcreteEuclideanJordanAlgebra
.random_instance