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
127 The multiplication table must be square with ``check_axioms=True``::
129 sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,((),()),((1,),))
130 Traceback (most recent call last):
132 ValueError: multiplication table is not square
134 The multiplication and inner-product tables must be the same
135 size (and in particular, the inner-product table must also be
136 square) with ``check_axioms=True``::
138 sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,((1,),),(()))
139 Traceback (most recent call last):
141 ValueError: multiplication and inner-product tables are
143 sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,((1,),),((1,2),))
144 Traceback (most recent call last):
146 ValueError: multiplication and inner-product tables are
151 if not field
.is_subring(RR
):
152 # Note: this does return true for the real algebraic
153 # field, the rationals, and any quadratic field where
154 # we've specified a real embedding.
155 raise ValueError("scalar field is not real")
158 # The multiplication and inner-product tables should be square
159 # if the user wants us to verify them. And we verify them as
160 # soon as possible, because we want to exploit their symmetry.
161 n
= len(multiplication_table
)
163 if not all( len(l
) == n
for l
in multiplication_table
):
164 raise ValueError("multiplication table is not square")
166 # If the multiplication table is square, we can check if
167 # the inner-product table is square by comparing it to the
168 # multiplication table's dimensions.
169 msg
= "multiplication and inner-product tables are different sizes"
170 if not len(inner_product_table
) == n
:
171 raise ValueError(msg
)
173 if not all( len(l
) == n
for l
in inner_product_table
):
174 raise ValueError(msg
)
176 if not all( multiplication_table
[j
][i
]
177 == multiplication_table
[i
][j
]
179 for j
in range(i
+1) ):
180 raise ValueError("Jordan product is not commutative")
181 if not all( inner_product_table
[j
][i
]
182 == inner_product_table
[i
][j
]
184 for j
in range(i
+1) ):
185 raise ValueError("inner-product is not commutative")
186 self
._matrix
_basis
= matrix_basis
189 category
= MagmaticAlgebras(field
).FiniteDimensional()
190 category
= category
.WithBasis().Unital()
192 fda
= super(FiniteDimensionalEuclideanJordanAlgebra
, self
)
197 self
.print_options(bracket
='')
199 # The multiplication table we're given is necessarily in terms
200 # of vectors, because we don't have an algebra yet for
201 # anything to be an element of. However, it's faster in the
202 # long run to have the multiplication table be in terms of
203 # algebra elements. We do this after calling the superclass
204 # constructor so that from_vector() knows what to do.
206 # Note: we take advantage of symmetry here, and only store
207 # the lower-triangular portion of the table.
208 self
._multiplication
_table
= [ [ self
.vector_space().zero()
209 for j
in range(i
+1) ]
214 elt
= self
.from_vector(multiplication_table
[i
][j
])
215 self
._multiplication
_table
[i
][j
] = elt
217 # Save our inner product as a matrix, since the efficiency of
218 # matrix multiplication will usually outweigh the fact that we
219 # have to store a redundant upper- or lower-triangular part.
220 # Pre-cache the fact that these are Hermitian (real symmetric,
221 # in fact) in case some e.g. matrix multiplication routine can
222 # take advantage of it.
223 self
._inner
_product
_matrix
= matrix(field
, inner_product_table
)
224 self
._inner
_product
_matrix
._cache
= {'hermitian': False}
227 if not self
._is
_jordanian
():
228 raise ValueError("Jordan identity does not hold")
229 if not self
._inner
_product
_is
_associative
():
230 raise ValueError("inner product is not associative")
232 def _element_constructor_(self
, elt
):
234 Construct an element of this algebra from its vector or matrix
237 This gets called only after the parent element _call_ method
238 fails to find a coercion for the argument.
242 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
244 ....: RealSymmetricEJA)
248 The identity in `S^n` is converted to the identity in the EJA::
250 sage: J = RealSymmetricEJA(3)
251 sage: I = matrix.identity(QQ,3)
252 sage: J(I) == J.one()
255 This skew-symmetric matrix can't be represented in the EJA::
257 sage: J = RealSymmetricEJA(3)
258 sage: A = matrix(QQ,3, lambda i,j: i-j)
260 Traceback (most recent call last):
262 ValueError: not an element of this algebra
266 Ensure that we can convert any element of the two non-matrix
267 simple algebras (whose matrix representations are columns)
268 back and forth faithfully::
270 sage: set_random_seed()
271 sage: J = HadamardEJA.random_instance()
272 sage: x = J.random_element()
273 sage: J(x.to_vector().column()) == x
275 sage: J = JordanSpinEJA.random_instance()
276 sage: x = J.random_element()
277 sage: J(x.to_vector().column()) == x
280 msg
= "not an element of this algebra"
282 # The superclass implementation of random_element()
283 # needs to be able to coerce "0" into the algebra.
285 elif elt
in self
.base_ring():
286 # Ensure that no base ring -> algebra coercion is performed
287 # by this method. There's some stupidity in sage that would
288 # otherwise propagate to this method; for example, sage thinks
289 # that the integer 3 belongs to the space of 2-by-2 matrices.
290 raise ValueError(msg
)
292 if elt
not in self
.matrix_space():
293 raise ValueError(msg
)
295 # Thanks for nothing! Matrix spaces aren't vector spaces in
296 # Sage, so we have to figure out its matrix-basis coordinates
297 # ourselves. We use the basis space's ring instead of the
298 # element's ring because the basis space might be an algebraic
299 # closure whereas the base ring of the 3-by-3 identity matrix
300 # could be QQ instead of QQbar.
301 V
= VectorSpace(self
.base_ring(), elt
.nrows()*elt
.ncols())
302 W
= V
.span_of_basis( _mat2vec(s
) for s
in self
.matrix_basis() )
305 coords
= W
.coordinate_vector(_mat2vec(elt
))
306 except ArithmeticError: # vector is not in free module
307 raise ValueError(msg
)
309 return self
.from_vector(coords
)
313 Return a string representation of ``self``.
317 sage: from mjo.eja.eja_algebra import JordanSpinEJA
321 Ensure that it says what we think it says::
323 sage: JordanSpinEJA(2, field=AA)
324 Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
325 sage: JordanSpinEJA(3, field=RDF)
326 Euclidean Jordan algebra of dimension 3 over Real Double Field
329 fmt
= "Euclidean Jordan algebra of dimension {} over {}"
330 return fmt
.format(self
.dimension(), self
.base_ring())
332 def product_on_basis(self
, i
, j
):
333 # We only stored the lower-triangular portion of the
334 # multiplication table.
336 return self
._multiplication
_table
[i
][j
]
338 return self
._multiplication
_table
[j
][i
]
340 def _is_commutative(self
):
342 Whether or not this algebra's multiplication table is commutative.
344 This method should of course always return ``True``, unless
345 this algebra was constructed with ``check_axioms=False`` and
346 passed an invalid multiplication table.
348 return all( self
.product_on_basis(i
,j
) == self
.product_on_basis(i
,j
)
349 for i
in range(self
.dimension())
350 for j
in range(self
.dimension()) )
352 def _is_jordanian(self
):
354 Whether or not this algebra's multiplication table respects the
355 Jordan identity `(x^{2})(xy) = x(x^{2}y)`.
357 We only check one arrangement of `x` and `y`, so for a
358 ``True`` result to be truly true, you should also check
359 :meth:`_is_commutative`. This method should of course always
360 return ``True``, unless this algebra was constructed with
361 ``check_axioms=False`` and passed an invalid multiplication table.
363 return all( (self
.monomial(i
)**2)*(self
.monomial(i
)*self
.monomial(j
))
365 (self
.monomial(i
))*((self
.monomial(i
)**2)*self
.monomial(j
))
366 for i
in range(self
.dimension())
367 for j
in range(self
.dimension()) )
369 def _inner_product_is_associative(self
):
371 Return whether or not this algebra's inner product `B` is
372 associative; that is, whether or not `B(xy,z) = B(x,yz)`.
374 This method should of course always return ``True``, unless
375 this algebra was constructed with ``check_axioms=False`` and
376 passed an invalid multiplication table.
379 # Used to check whether or not something is zero in an inexact
380 # ring. This number is sufficient to allow the construction of
381 # QuaternionHermitianEJA(2, RDF) with check_axioms=True.
384 for i
in range(self
.dimension()):
385 for j
in range(self
.dimension()):
386 for k
in range(self
.dimension()):
390 diff
= (x
*y
).inner_product(z
) - x
.inner_product(y
*z
)
392 if self
.base_ring().is_exact():
396 if diff
.abs() > epsilon
:
402 def characteristic_polynomial_of(self
):
404 Return the algebra's "characteristic polynomial of" function,
405 which is itself a multivariate polynomial that, when evaluated
406 at the coordinates of some algebra element, returns that
407 element's characteristic polynomial.
409 The resulting polynomial has `n+1` variables, where `n` is the
410 dimension of this algebra. The first `n` variables correspond to
411 the coordinates of an algebra element: when evaluated at the
412 coordinates of an algebra element with respect to a certain
413 basis, the result is a univariate polynomial (in the one
414 remaining variable ``t``), namely the characteristic polynomial
419 sage: from mjo.eja.eja_algebra import JordanSpinEJA, TrivialEJA
423 The characteristic polynomial in the spin algebra is given in
424 Alizadeh, Example 11.11::
426 sage: J = JordanSpinEJA(3)
427 sage: p = J.characteristic_polynomial_of(); p
428 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
429 sage: xvec = J.one().to_vector()
433 By definition, the characteristic polynomial is a monic
434 degree-zero polynomial in a rank-zero algebra. Note that
435 Cayley-Hamilton is indeed satisfied since the polynomial
436 ``1`` evaluates to the identity element of the algebra on
439 sage: J = TrivialEJA()
440 sage: J.characteristic_polynomial_of()
447 # The list of coefficient polynomials a_0, a_1, a_2, ..., a_(r-1).
448 a
= self
._charpoly
_coefficients
()
450 # We go to a bit of trouble here to reorder the
451 # indeterminates, so that it's easier to evaluate the
452 # characteristic polynomial at x's coordinates and get back
453 # something in terms of t, which is what we want.
454 S
= PolynomialRing(self
.base_ring(),'t')
458 S
= PolynomialRing(S
, R
.variable_names())
461 return (t
**r
+ sum( a
[k
]*(t
**k
) for k
in range(r
) ))
463 def coordinate_polynomial_ring(self
):
465 The multivariate polynomial ring in which this algebra's
466 :meth:`characteristic_polynomial_of` lives.
470 sage: from mjo.eja.eja_algebra import (HadamardEJA,
471 ....: RealSymmetricEJA)
475 sage: J = HadamardEJA(2)
476 sage: J.coordinate_polynomial_ring()
477 Multivariate Polynomial Ring in X1, X2...
478 sage: J = RealSymmetricEJA(3,QQ)
479 sage: J.coordinate_polynomial_ring()
480 Multivariate Polynomial Ring in X1, X2, X3, X4, X5, X6...
483 var_names
= tuple( "X%d" % z
for z
in range(1, self
.dimension()+1) )
484 return PolynomialRing(self
.base_ring(), var_names
)
486 def inner_product(self
, x
, y
):
488 The inner product associated with this Euclidean Jordan algebra.
490 Defaults to the trace inner product, but can be overridden by
491 subclasses if they are sure that the necessary properties are
496 sage: from mjo.eja.eja_algebra import (random_eja,
498 ....: BilinearFormEJA)
502 Our inner product is "associative," which means the following for
503 a symmetric bilinear form::
505 sage: set_random_seed()
506 sage: J = random_eja()
507 sage: x,y,z = J.random_elements(3)
508 sage: (x*y).inner_product(z) == y.inner_product(x*z)
513 Ensure that this is the usual inner product for the algebras
516 sage: set_random_seed()
517 sage: J = HadamardEJA.random_instance()
518 sage: x,y = J.random_elements(2)
519 sage: actual = x.inner_product(y)
520 sage: expected = x.to_vector().inner_product(y.to_vector())
521 sage: actual == expected
524 Ensure that this is one-half of the trace inner-product in a
525 BilinearFormEJA that isn't just the reals (when ``n`` isn't
526 one). This is in Faraut and Koranyi, and also my "On the
529 sage: set_random_seed()
530 sage: J = BilinearFormEJA.random_instance()
531 sage: n = J.dimension()
532 sage: x = J.random_element()
533 sage: y = J.random_element()
534 sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
537 B
= self
._inner
_product
_matrix
538 return (B
*x
.to_vector()).inner_product(y
.to_vector())
541 def is_trivial(self
):
543 Return whether or not this algebra is trivial.
545 A trivial algebra contains only the zero element.
549 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
554 sage: J = ComplexHermitianEJA(3)
560 sage: J = TrivialEJA()
565 return self
.dimension() == 0
568 def multiplication_table(self
):
570 Return a visual representation of this algebra's multiplication
571 table (on basis elements).
575 sage: from mjo.eja.eja_algebra import JordanSpinEJA
579 sage: J = JordanSpinEJA(4)
580 sage: J.multiplication_table()
581 +----++----+----+----+----+
582 | * || e0 | e1 | e2 | e3 |
583 +====++====+====+====+====+
584 | e0 || e0 | e1 | e2 | e3 |
585 +----++----+----+----+----+
586 | e1 || e1 | e0 | 0 | 0 |
587 +----++----+----+----+----+
588 | e2 || e2 | 0 | e0 | 0 |
589 +----++----+----+----+----+
590 | e3 || e3 | 0 | 0 | e0 |
591 +----++----+----+----+----+
595 M
= [ [ self
.zero() for j
in range(n
) ]
599 M
[i
][j
] = self
._multiplication
_table
[i
][j
]
603 # Prepend the left "header" column entry Can't do this in
604 # the loop because it messes up the symmetry.
605 M
[i
] = [self
.monomial(i
)] + M
[i
]
607 # Prepend the header row.
608 M
= [["*"] + list(self
.gens())] + M
609 return table(M
, header_row
=True, header_column
=True, frame
=True)
612 def matrix_basis(self
):
614 Return an (often more natural) representation of this algebras
615 basis as an ordered tuple of matrices.
617 Every finite-dimensional Euclidean Jordan Algebra is a, up to
618 Jordan isomorphism, a direct sum of five simple
619 algebras---four of which comprise Hermitian matrices. And the
620 last type of algebra can of course be thought of as `n`-by-`1`
621 column matrices (ambiguusly called column vectors) to avoid
622 special cases. As a result, matrices (and column vectors) are
623 a natural representation format for Euclidean Jordan algebra
626 But, when we construct an algebra from a basis of matrices,
627 those matrix representations are lost in favor of coordinate
628 vectors *with respect to* that basis. We could eventually
629 convert back if we tried hard enough, but having the original
630 representations handy is valuable enough that we simply store
631 them and return them from this method.
633 Why implement this for non-matrix algebras? Avoiding special
634 cases for the :class:`BilinearFormEJA` pays with simplicity in
635 its own right. But mainly, we would like to be able to assume
636 that elements of a :class:`DirectSumEJA` can be displayed
637 nicely, without having to have special classes for direct sums
638 one of whose components was a matrix algebra.
642 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
643 ....: RealSymmetricEJA)
647 sage: J = RealSymmetricEJA(2)
649 Finite family {0: e0, 1: e1, 2: e2}
650 sage: J.matrix_basis()
652 [1 0] [ 0 0.7071067811865475?] [0 0]
653 [0 0], [0.7071067811865475? 0], [0 1]
658 sage: J = JordanSpinEJA(2)
660 Finite family {0: e0, 1: e1}
661 sage: J.matrix_basis()
667 if self
._matrix
_basis
is None:
668 M
= self
.matrix_space()
669 return tuple( M(b
.to_vector()) for b
in self
.basis() )
671 return self
._matrix
_basis
674 def matrix_space(self
):
676 Return the matrix space in which this algebra's elements live, if
677 we think of them as matrices (including column vectors of the
680 Generally this will be an `n`-by-`1` column-vector space,
681 except when the algebra is trivial. There it's `n`-by-`n`
682 (where `n` is zero), to ensure that two elements of the matrix
683 space (empty matrices) can be multiplied.
685 Matrix algebras override this with something more useful.
687 if self
.is_trivial():
688 return MatrixSpace(self
.base_ring(), 0)
689 elif self
._matrix
_basis
is None or len(self
._matrix
_basis
) == 0:
690 return MatrixSpace(self
.base_ring(), self
.dimension(), 1)
692 return self
._matrix
_basis
[0].matrix_space()
698 Return the unit element of this algebra.
702 sage: from mjo.eja.eja_algebra import (HadamardEJA,
707 sage: J = HadamardEJA(5)
709 e0 + e1 + e2 + e3 + e4
713 The identity element acts like the identity::
715 sage: set_random_seed()
716 sage: J = random_eja()
717 sage: x = J.random_element()
718 sage: J.one()*x == x and x*J.one() == x
721 The matrix of the unit element's operator is the identity::
723 sage: set_random_seed()
724 sage: J = random_eja()
725 sage: actual = J.one().operator().matrix()
726 sage: expected = matrix.identity(J.base_ring(), J.dimension())
727 sage: actual == expected
730 Ensure that the cached unit element (often precomputed by
731 hand) agrees with the computed one::
733 sage: set_random_seed()
734 sage: J = random_eja()
735 sage: cached = J.one()
736 sage: J.one.clear_cache()
737 sage: J.one() == cached
741 # We can brute-force compute the matrices of the operators
742 # that correspond to the basis elements of this algebra.
743 # If some linear combination of those basis elements is the
744 # algebra identity, then the same linear combination of
745 # their matrices has to be the identity matrix.
747 # Of course, matrices aren't vectors in sage, so we have to
748 # appeal to the "long vectors" isometry.
749 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
751 # Now we use basic linear algebra to find the coefficients,
752 # of the matrices-as-vectors-linear-combination, which should
753 # work for the original algebra basis too.
754 A
= matrix(self
.base_ring(), oper_vecs
)
756 # We used the isometry on the left-hand side already, but we
757 # still need to do it for the right-hand side. Recall that we
758 # wanted something that summed to the identity matrix.
759 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
761 # Now if there's an identity element in the algebra, this
762 # should work. We solve on the left to avoid having to
763 # transpose the matrix "A".
764 return self
.from_vector(A
.solve_left(b
))
767 def peirce_decomposition(self
, c
):
769 The Peirce decomposition of this algebra relative to the
772 In the future, this can be extended to a complete system of
773 orthogonal idempotents.
777 - ``c`` -- an idempotent of this algebra.
781 A triple (J0, J5, J1) containing two subalgebras and one subspace
784 - ``J0`` -- the algebra on the eigenspace of ``c.operator()``
785 corresponding to the eigenvalue zero.
787 - ``J5`` -- the eigenspace (NOT a subalgebra) of ``c.operator()``
788 corresponding to the eigenvalue one-half.
790 - ``J1`` -- the algebra on the eigenspace of ``c.operator()``
791 corresponding to the eigenvalue one.
793 These are the only possible eigenspaces for that operator, and this
794 algebra is a direct sum of them. The spaces ``J0`` and ``J1`` are
795 orthogonal, and are subalgebras of this algebra with the appropriate
800 sage: from mjo.eja.eja_algebra import random_eja, RealSymmetricEJA
804 The canonical example comes from the symmetric matrices, which
805 decompose into diagonal and off-diagonal parts::
807 sage: J = RealSymmetricEJA(3)
808 sage: C = matrix(QQ, [ [1,0,0],
812 sage: J0,J5,J1 = J.peirce_decomposition(c)
814 Euclidean Jordan algebra of dimension 1...
816 Vector space of degree 6 and dimension 2...
818 Euclidean Jordan algebra of dimension 3...
819 sage: J0.one().to_matrix()
823 sage: orig_df = AA.options.display_format
824 sage: AA.options.display_format = 'radical'
825 sage: J.from_vector(J5.basis()[0]).to_matrix()
829 sage: J.from_vector(J5.basis()[1]).to_matrix()
833 sage: AA.options.display_format = orig_df
834 sage: J1.one().to_matrix()
841 Every algebra decomposes trivially with respect to its identity
844 sage: set_random_seed()
845 sage: J = random_eja()
846 sage: J0,J5,J1 = J.peirce_decomposition(J.one())
847 sage: J0.dimension() == 0 and J5.dimension() == 0
849 sage: J1.superalgebra() == J and J1.dimension() == J.dimension()
852 The decomposition is into eigenspaces, and its components are
853 therefore necessarily orthogonal. Moreover, the identity
854 elements in the two subalgebras are the projections onto their
855 respective subspaces of the superalgebra's identity element::
857 sage: set_random_seed()
858 sage: J = random_eja()
859 sage: x = J.random_element()
860 sage: if not J.is_trivial():
861 ....: while x.is_nilpotent():
862 ....: x = J.random_element()
863 sage: c = x.subalgebra_idempotent()
864 sage: J0,J5,J1 = J.peirce_decomposition(c)
866 sage: for (w,y,z) in zip(J0.basis(), J5.basis(), J1.basis()):
867 ....: w = w.superalgebra_element()
868 ....: y = J.from_vector(y)
869 ....: z = z.superalgebra_element()
870 ....: ipsum += w.inner_product(y).abs()
871 ....: ipsum += w.inner_product(z).abs()
872 ....: ipsum += y.inner_product(z).abs()
875 sage: J1(c) == J1.one()
877 sage: J0(J.one() - c) == J0.one()
881 if not c
.is_idempotent():
882 raise ValueError("element is not idempotent: %s" % c
)
884 from mjo
.eja
.eja_subalgebra
import FiniteDimensionalEuclideanJordanSubalgebra
886 # Default these to what they should be if they turn out to be
887 # trivial, because eigenspaces_left() won't return eigenvalues
888 # corresponding to trivial spaces (e.g. it returns only the
889 # eigenspace corresponding to lambda=1 if you take the
890 # decomposition relative to the identity element).
891 trivial
= FiniteDimensionalEuclideanJordanSubalgebra(self
, ())
892 J0
= trivial
# eigenvalue zero
893 J5
= VectorSpace(self
.base_ring(), 0) # eigenvalue one-half
894 J1
= trivial
# eigenvalue one
896 for (eigval
, eigspace
) in c
.operator().matrix().right_eigenspaces():
897 if eigval
== ~
(self
.base_ring()(2)):
900 gens
= tuple( self
.from_vector(b
) for b
in eigspace
.basis() )
901 subalg
= FiniteDimensionalEuclideanJordanSubalgebra(self
,
909 raise ValueError("unexpected eigenvalue: %s" % eigval
)
914 def random_element(self
, thorough
=False):
916 Return a random element of this algebra.
918 Our algebra superclass method only returns a linear
919 combination of at most two basis elements. We instead
920 want the vector space "random element" method that
921 returns a more diverse selection.
925 - ``thorough`` -- (boolean; default False) whether or not we
926 should generate irrational coefficients for the random
927 element when our base ring is irrational; this slows the
928 algebra operations to a crawl, but any truly random method
932 # For a general base ring... maybe we can trust this to do the
933 # right thing? Unlikely, but.
934 V
= self
.vector_space()
935 v
= V
.random_element()
937 if self
.base_ring() is AA
:
938 # The "random element" method of the algebraic reals is
939 # stupid at the moment, and only returns integers between
940 # -2 and 2, inclusive:
942 # https://trac.sagemath.org/ticket/30875
944 # Instead, we implement our own "random vector" method,
945 # and then coerce that into the algebra. We use the vector
946 # space degree here instead of the dimension because a
947 # subalgebra could (for example) be spanned by only two
948 # vectors, each with five coordinates. We need to
949 # generate all five coordinates.
951 v
*= QQbar
.random_element().real()
953 v
*= QQ
.random_element()
955 return self
.from_vector(V
.coordinate_vector(v
))
957 def random_elements(self
, count
, thorough
=False):
959 Return ``count`` random elements as a tuple.
963 - ``thorough`` -- (boolean; default False) whether or not we
964 should generate irrational coefficients for the random
965 elements when our base ring is irrational; this slows the
966 algebra operations to a crawl, but any truly random method
971 sage: from mjo.eja.eja_algebra import JordanSpinEJA
975 sage: J = JordanSpinEJA(3)
976 sage: x,y,z = J.random_elements(3)
977 sage: all( [ x in J, y in J, z in J ])
979 sage: len( J.random_elements(10) ) == 10
983 return tuple( self
.random_element(thorough
)
984 for idx
in range(count
) )
988 def _charpoly_coefficients(self
):
990 The `r` polynomial coefficients of the "characteristic polynomial
994 R
= self
.coordinate_polynomial_ring()
996 F
= R
.fraction_field()
999 # From a result in my book, these are the entries of the
1000 # basis representation of L_x.
1001 return sum( vars[k
]*self
.monomial(k
).operator().matrix()[i
,j
]
1004 L_x
= matrix(F
, n
, n
, L_x_i_j
)
1007 if self
.rank
.is_in_cache():
1009 # There's no need to pad the system with redundant
1010 # columns if we *know* they'll be redundant.
1013 # Compute an extra power in case the rank is equal to
1014 # the dimension (otherwise, we would stop at x^(r-1)).
1015 x_powers
= [ (L_x
**k
)*self
.one().to_vector()
1016 for k
in range(n
+1) ]
1017 A
= matrix
.column(F
, x_powers
[:n
])
1018 AE
= A
.extended_echelon_form()
1025 # The theory says that only the first "r" coefficients are
1026 # nonzero, and they actually live in the original polynomial
1027 # ring and not the fraction field. We negate them because
1028 # in the actual characteristic polynomial, they get moved
1029 # to the other side where x^r lives.
1030 return -A_rref
.solve_right(E
*b
).change_ring(R
)[:r
]
1035 Return the rank of this EJA.
1037 This is a cached method because we know the rank a priori for
1038 all of the algebras we can construct. Thus we can avoid the
1039 expensive ``_charpoly_coefficients()`` call unless we truly
1040 need to compute the whole characteristic polynomial.
1044 sage: from mjo.eja.eja_algebra import (HadamardEJA,
1045 ....: JordanSpinEJA,
1046 ....: RealSymmetricEJA,
1047 ....: ComplexHermitianEJA,
1048 ....: QuaternionHermitianEJA,
1053 The rank of the Jordan spin algebra is always two::
1055 sage: JordanSpinEJA(2).rank()
1057 sage: JordanSpinEJA(3).rank()
1059 sage: JordanSpinEJA(4).rank()
1062 The rank of the `n`-by-`n` Hermitian real, complex, or
1063 quaternion matrices is `n`::
1065 sage: RealSymmetricEJA(4).rank()
1067 sage: ComplexHermitianEJA(3).rank()
1069 sage: QuaternionHermitianEJA(2).rank()
1074 Ensure that every EJA that we know how to construct has a
1075 positive integer rank, unless the algebra is trivial in
1076 which case its rank will be zero::
1078 sage: set_random_seed()
1079 sage: J = random_eja()
1083 sage: r > 0 or (r == 0 and J.is_trivial())
1086 Ensure that computing the rank actually works, since the ranks
1087 of all simple algebras are known and will be cached by default::
1089 sage: set_random_seed() # long time
1090 sage: J = random_eja() # long time
1091 sage: caches = J.rank() # long time
1092 sage: J.rank.clear_cache() # long time
1093 sage: J.rank() == cached # long time
1097 return len(self
._charpoly
_coefficients
())
1100 def vector_space(self
):
1102 Return the vector space that underlies this algebra.
1106 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1110 sage: J = RealSymmetricEJA(2)
1111 sage: J.vector_space()
1112 Vector space of dimension 3 over...
1115 return self
.zero().to_vector().parent().ambient_vector_space()
1118 Element
= FiniteDimensionalEuclideanJordanAlgebraElement
1120 class RationalBasisEuclideanJordanAlgebraNg(FiniteDimensionalEuclideanJordanAlgebra
):
1122 New class for algebras whose supplied basis elements have all rational entries.
1126 sage: from mjo.eja.eja_algebra import BilinearFormEJA
1130 The supplied basis is orthonormalized by default::
1132 sage: B = matrix(QQ, [[1, 0, 0], [0, 25, -32], [0, -32, 41]])
1133 sage: J = BilinearFormEJA(B)
1134 sage: J.matrix_basis()
1147 orthonormalize
=True,
1154 vector_basis
= basis
1156 from sage
.structure
.element
import is_Matrix
1157 basis_is_matrices
= False
1161 if is_Matrix(basis
[0]):
1162 basis_is_matrices
= True
1163 vector_basis
= tuple( map(_mat2vec
,basis
) )
1164 degree
= basis
[0].nrows()**2
1166 degree
= basis
[0].degree()
1168 V
= VectorSpace(field
, degree
)
1170 # If we were asked to orthonormalize, and if the orthonormal
1171 # basis is different from the given one, then we also want to
1172 # compute multiplication and inner-product tables for the
1173 # deorthonormalized basis. These can be used later to
1174 # construct a deorthonormalized copy of this algebra over QQ
1175 # in which several operations are much faster.
1176 self
._deortho
_multiplication
_table
= None
1177 self
._deortho
_inner
_product
_table
= None
1180 from mjo
.eja
.eja_utils
import gram_schmidt
1181 vector_basis
= gram_schmidt(vector_basis
, inner_product
)
1182 W
= V
.span_of_basis( vector_basis
)
1184 # Normalize the "matrix" basis, too!
1185 basis
= vector_basis
1187 if basis_is_matrices
:
1188 from mjo
.eja
.eja_utils
import _vec2mat
1189 basis
= tuple( map(_vec2mat
,basis
) )
1191 W
= V
.span_of_basis( vector_basis
)
1193 mult_table
= [ [0 for i
in range(n
)] for j
in range(n
) ]
1194 ip_table
= [ [0 for i
in range(n
)] for j
in range(n
) ]
1196 # Note: the Jordan and inner- products are defined in terms
1197 # of the ambient basis. It's important that their arguments
1198 # are in ambient coordinates as well.
1200 for j
in range(i
+1):
1201 # ortho basis w.r.t. ambient coords
1202 q_i
= vector_basis
[i
]
1203 q_j
= vector_basis
[j
]
1205 if basis_is_matrices
:
1209 elt
= jordan_product(q_i
, q_j
)
1210 ip
= inner_product(q_i
, q_j
)
1212 if basis_is_matrices
:
1213 # do another mat2vec because the multiplication
1214 # table is in terms of vectors
1217 elt
= W
.coordinate_vector(elt
)
1218 mult_table
[i
][j
] = elt
1219 mult_table
[j
][i
] = elt
1223 if basis_is_matrices
:
1227 basis
= tuple( x
.column() for x
in basis
)
1229 super().__init
__(field
,
1234 basis
, # matrix basis
1238 class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra
):
1240 Algebras whose basis consists of vectors with rational
1241 entries. Equivalently, algebras whose multiplication tables
1242 contain only rational coefficients.
1244 When an EJA has a basis that can be made rational, we can speed up
1245 the computation of its characteristic polynomial by doing it over
1246 ``QQ``. All of the named EJA constructors that we provide fall
1250 def _charpoly_coefficients(self
):
1252 Override the parent method with something that tries to compute
1253 over a faster (non-extension) field.
1257 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1261 The base ring of the resulting polynomial coefficients is what
1262 it should be, and not the rationals (unless the algebra was
1263 already over the rationals)::
1265 sage: J = JordanSpinEJA(3)
1266 sage: J._charpoly_coefficients()
1267 (X1^2 - X2^2 - X3^2, -2*X1)
1268 sage: a0 = J._charpoly_coefficients()[0]
1270 Algebraic Real Field
1271 sage: a0.base_ring()
1272 Algebraic Real Field
1275 if self
.base_ring() is QQ
:
1276 # There's no need to construct *another* algebra over the
1277 # rationals if this one is already over the rationals.
1278 superclass
= super(RationalBasisEuclideanJordanAlgebra
, self
)
1279 return superclass
._charpoly
_coefficients
()
1282 tuple(map(lambda x
: x
.to_vector(), ls
))
1283 for ls
in self
._multiplication
_table
1286 # Do the computation over the rationals. The answer will be
1287 # the same, because our basis coordinates are (essentially)
1289 J
= FiniteDimensionalEuclideanJordanAlgebra(QQ
,
1293 a
= J
._charpoly
_coefficients
()
1294 return tuple(map(lambda x
: x
.change_ring(self
.base_ring()), a
))
1297 class ConcreteEuclideanJordanAlgebra
:
1299 A class for the Euclidean Jordan algebras that we know by name.
1301 These are the Jordan algebras whose basis, multiplication table,
1302 rank, and so on are known a priori. More to the point, they are
1303 the Euclidean Jordan algebras for which we are able to conjure up
1304 a "random instance."
1308 sage: from mjo.eja.eja_algebra import ConcreteEuclideanJordanAlgebra
1312 Our basis is normalized with respect to the algebra's inner
1313 product, unless we specify otherwise::
1315 sage: set_random_seed()
1316 sage: J = ConcreteEuclideanJordanAlgebra.random_instance()
1317 sage: all( b.norm() == 1 for b in J.gens() )
1320 Since our basis is orthonormal with respect to the algebra's inner
1321 product, and since we know that this algebra is an EJA, any
1322 left-multiplication operator's matrix will be symmetric because
1323 natural->EJA basis representation is an isometry and within the
1324 EJA the operator is self-adjoint by the Jordan axiom::
1326 sage: set_random_seed()
1327 sage: J = ConcreteEuclideanJordanAlgebra.random_instance()
1328 sage: x = J.random_element()
1329 sage: x.operator().is_self_adjoint()
1334 def _max_random_instance_size():
1336 Return an integer "size" that is an upper bound on the size of
1337 this algebra when it is used in a random test
1338 case. Unfortunately, the term "size" is ambiguous -- when
1339 dealing with `R^n` under either the Hadamard or Jordan spin
1340 product, the "size" refers to the dimension `n`. When dealing
1341 with a matrix algebra (real symmetric or complex/quaternion
1342 Hermitian), it refers to the size of the matrix, which is far
1343 less than the dimension of the underlying vector space.
1345 This method must be implemented in each subclass.
1347 raise NotImplementedError
1350 def random_instance(cls
, field
=AA
, **kwargs
):
1352 Return a random instance of this type of algebra.
1354 This method should be implemented in each subclass.
1356 from sage
.misc
.prandom
import choice
1357 eja_class
= choice(cls
.__subclasses
__())
1358 return eja_class
.random_instance(field
)
1361 class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra
):
1363 def __init__(self
, field
, basis
, normalize_basis
=True, **kwargs
):
1365 Compared to the superclass constructor, we take a basis instead of
1366 a multiplication table because the latter can be computed in terms
1367 of the former when the product is known (like it is here).
1369 # Used in this class's fast _charpoly_coefficients() override.
1370 self
._basis
_normalizers
= None
1372 # We're going to loop through this a few times, so now's a good
1373 # time to ensure that it isn't a generator expression.
1374 basis
= tuple(basis
)
1376 algebra_dim
= len(basis
)
1377 degree
= 0 # size of the matrices
1379 degree
= basis
[0].nrows()
1381 if algebra_dim
> 1 and normalize_basis
:
1382 # We'll need sqrt(2) to normalize the basis, and this
1383 # winds up in the multiplication table, so the whole
1384 # algebra needs to be over the field extension.
1385 R
= PolynomialRing(field
, 'z')
1388 if p
.is_irreducible():
1389 field
= field
.extension(p
, 'sqrt2', embedding
=RLF(2).sqrt())
1390 basis
= tuple( s
.change_ring(field
) for s
in basis
)
1391 self
._basis
_normalizers
= tuple(
1392 ~
(self
.matrix_inner_product(s
,s
).sqrt()) for s
in basis
)
1393 basis
= tuple(s
*c
for (s
,c
) in zip(basis
,self
._basis
_normalizers
))
1395 # Now compute the multiplication and inner product tables.
1396 # We have to do this *after* normalizing the basis, because
1397 # scaling affects the answers.
1398 V
= VectorSpace(field
, degree
**2)
1399 W
= V
.span_of_basis( _mat2vec(s
) for s
in basis
)
1400 mult_table
= [[W
.zero() for j
in range(algebra_dim
)]
1401 for i
in range(algebra_dim
)]
1402 ip_table
= [[field
.zero() for j
in range(algebra_dim
)]
1403 for i
in range(algebra_dim
)]
1404 for i
in range(algebra_dim
):
1405 for j
in range(algebra_dim
):
1406 mat_entry
= (basis
[i
]*basis
[j
] + basis
[j
]*basis
[i
])/2
1407 mult_table
[i
][j
] = W
.coordinate_vector(_mat2vec(mat_entry
))
1410 # HACK: ignore the error here if we don't need the
1411 # inner product (as is the case when we construct
1412 # a dummy QQ-algebra for fast charpoly coefficients.
1413 ip_table
[i
][j
] = self
.matrix_inner_product(basis
[i
],
1418 super(MatrixEuclideanJordanAlgebra
, self
).__init
__(field
,
1424 if algebra_dim
== 0:
1425 self
.one
.set_cache(self
.zero())
1427 n
= basis
[0].nrows()
1428 # The identity wrt (A,B) -> (AB + BA)/2 is independent of the
1429 # details of this algebra.
1430 self
.one
.set_cache(self(matrix
.identity(field
,n
)))
1434 def _charpoly_coefficients(self
):
1436 Override the parent method with something that tries to compute
1437 over a faster (non-extension) field.
1439 if self
._basis
_normalizers
is None or self
.base_ring() is QQ
:
1440 # We didn't normalize, or the basis we started with had
1441 # entries in a nice field already. Just compute the thing.
1442 return super(MatrixEuclideanJordanAlgebra
, self
)._charpoly
_coefficients
()
1444 basis
= ( (b
/n
) for (b
,n
) in zip(self
.matrix_basis(),
1445 self
._basis
_normalizers
) )
1447 # Do this over the rationals and convert back at the end.
1448 # Only works because we know the entries of the basis are
1449 # integers. The argument ``check_axioms=False`` is required
1450 # because the trace inner-product method for this
1451 # class is a stub and can't actually be checked.
1452 J
= MatrixEuclideanJordanAlgebra(QQ
,
1454 normalize_basis
=False,
1457 a
= J
._charpoly
_coefficients
()
1459 # Unfortunately, changing the basis does change the
1460 # coefficients of the characteristic polynomial, but since
1461 # these are really the coefficients of the "characteristic
1462 # polynomial of" function, everything is still nice and
1463 # unevaluated. It's therefore "obvious" how scaling the
1464 # basis affects the coordinate variables X1, X2, et
1465 # cetera. Scaling the first basis vector up by "n" adds a
1466 # factor of 1/n into every "X1" term, for example. So here
1467 # we simply undo the basis_normalizer scaling that we
1468 # performed earlier.
1470 # The a[0] access here is safe because trivial algebras
1471 # won't have any basis normalizers and therefore won't
1472 # make it to this "else" branch.
1473 XS
= a
[0].parent().gens()
1474 subs_dict
= { XS
[i
]: self
._basis
_normalizers
[i
]*XS
[i
]
1475 for i
in range(len(XS
)) }
1476 return tuple( a_i
.subs(subs_dict
) for a_i
in a
)
1482 Embed the matrix ``M`` into a space of real matrices.
1484 The matrix ``M`` can have entries in any field at the moment:
1485 the real numbers, complex numbers, or quaternions. And although
1486 they are not a field, we can probably support octonions at some
1487 point, too. This function returns a real matrix that "acts like"
1488 the original with respect to matrix multiplication; i.e.
1490 real_embed(M*N) = real_embed(M)*real_embed(N)
1493 raise NotImplementedError
1497 def real_unembed(M
):
1499 The inverse of :meth:`real_embed`.
1501 raise NotImplementedError
1504 def matrix_inner_product(cls
,X
,Y
):
1505 Xu
= cls
.real_unembed(X
)
1506 Yu
= cls
.real_unembed(Y
)
1507 tr
= (Xu
*Yu
).trace()
1510 # Works in QQ, AA, RDF, et cetera.
1512 except AttributeError:
1513 # A quaternion doesn't have a real() method, but does
1514 # have coefficient_tuple() method that returns the
1515 # coefficients of 1, i, j, and k -- in that order.
1516 return tr
.coefficient_tuple()[0]
1519 class RealMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1523 The identity function, for embedding real matrices into real
1529 def real_unembed(M
):
1531 The identity function, for unembedding real matrices from real
1537 class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra
,
1538 ConcreteEuclideanJordanAlgebra
):
1540 The rank-n simple EJA consisting of real symmetric n-by-n
1541 matrices, the usual symmetric Jordan product, and the trace inner
1542 product. It has dimension `(n^2 + n)/2` over the reals.
1546 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1550 sage: J = RealSymmetricEJA(2)
1551 sage: e0, e1, e2 = J.gens()
1559 In theory, our "field" can be any subfield of the reals::
1561 sage: RealSymmetricEJA(2, RDF)
1562 Euclidean Jordan algebra of dimension 3 over Real Double Field
1563 sage: RealSymmetricEJA(2, RR)
1564 Euclidean Jordan algebra of dimension 3 over Real Field with
1565 53 bits of precision
1569 The dimension of this algebra is `(n^2 + n) / 2`::
1571 sage: set_random_seed()
1572 sage: n_max = RealSymmetricEJA._max_random_instance_size()
1573 sage: n = ZZ.random_element(1, n_max)
1574 sage: J = RealSymmetricEJA(n)
1575 sage: J.dimension() == (n^2 + n)/2
1578 The Jordan multiplication is what we think it is::
1580 sage: set_random_seed()
1581 sage: J = RealSymmetricEJA.random_instance()
1582 sage: x,y = J.random_elements(2)
1583 sage: actual = (x*y).to_matrix()
1584 sage: X = x.to_matrix()
1585 sage: Y = y.to_matrix()
1586 sage: expected = (X*Y + Y*X)/2
1587 sage: actual == expected
1589 sage: J(expected) == x*y
1592 We can change the generator prefix::
1594 sage: RealSymmetricEJA(3, prefix='q').gens()
1595 (q0, q1, q2, q3, q4, q5)
1597 We can construct the (trivial) algebra of rank zero::
1599 sage: RealSymmetricEJA(0)
1600 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1604 def _denormalized_basis(cls
, n
, field
):
1606 Return a basis for the space of real symmetric n-by-n matrices.
1610 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1614 sage: set_random_seed()
1615 sage: n = ZZ.random_element(1,5)
1616 sage: B = RealSymmetricEJA._denormalized_basis(n,QQ)
1617 sage: all( M.is_symmetric() for M in B)
1621 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1625 for j
in range(i
+1):
1626 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
1630 Sij
= Eij
+ Eij
.transpose()
1636 def _max_random_instance_size():
1637 return 4 # Dimension 10
1640 def random_instance(cls
, field
=AA
, **kwargs
):
1642 Return a random instance of this type of algebra.
1644 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1645 return cls(n
, field
, **kwargs
)
1647 def __init__(self
, n
, field
=AA
, **kwargs
):
1648 basis
= self
._denormalized
_basis
(n
, field
)
1649 super(RealSymmetricEJA
, self
).__init
__(field
,
1653 self
.rank
.set_cache(n
)
1656 class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1660 Embed the n-by-n complex matrix ``M`` into the space of real
1661 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1662 bi` to the block matrix ``[[a,b],[-b,a]]``.
1666 sage: from mjo.eja.eja_algebra import \
1667 ....: ComplexMatrixEuclideanJordanAlgebra
1671 sage: F = QuadraticField(-1, 'I')
1672 sage: x1 = F(4 - 2*i)
1673 sage: x2 = F(1 + 2*i)
1676 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1677 sage: ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1686 Embedding is a homomorphism (isomorphism, in fact)::
1688 sage: set_random_seed()
1689 sage: n = ZZ.random_element(3)
1690 sage: F = QuadraticField(-1, 'I')
1691 sage: X = random_matrix(F, n)
1692 sage: Y = random_matrix(F, n)
1693 sage: Xe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X)
1694 sage: Ye = ComplexMatrixEuclideanJordanAlgebra.real_embed(Y)
1695 sage: XYe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1702 raise ValueError("the matrix 'M' must be square")
1704 # We don't need any adjoined elements...
1705 field
= M
.base_ring().base_ring()
1709 a
= z
.list()[0] # real part, I guess
1710 b
= z
.list()[1] # imag part, I guess
1711 blocks
.append(matrix(field
, 2, [[a
,b
],[-b
,a
]]))
1713 return matrix
.block(field
, n
, blocks
)
1717 def real_unembed(M
):
1719 The inverse of _embed_complex_matrix().
1723 sage: from mjo.eja.eja_algebra import \
1724 ....: ComplexMatrixEuclideanJordanAlgebra
1728 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1729 ....: [-2, 1, -4, 3],
1730 ....: [ 9, 10, 11, 12],
1731 ....: [-10, 9, -12, 11] ])
1732 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(A)
1734 [ 10*I + 9 12*I + 11]
1738 Unembedding is the inverse of embedding::
1740 sage: set_random_seed()
1741 sage: F = QuadraticField(-1, 'I')
1742 sage: M = random_matrix(F, 3)
1743 sage: Me = ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1744 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1750 raise ValueError("the matrix 'M' must be square")
1751 if not n
.mod(2).is_zero():
1752 raise ValueError("the matrix 'M' must be a complex embedding")
1754 # If "M" was normalized, its base ring might have roots
1755 # adjoined and they can stick around after unembedding.
1756 field
= M
.base_ring()
1757 R
= PolynomialRing(field
, 'z')
1760 # Sage doesn't know how to embed AA into QQbar, i.e. how
1761 # to adjoin sqrt(-1) to AA.
1764 F
= field
.extension(z
**2 + 1, 'I', embedding
=CLF(-1).sqrt())
1767 # Go top-left to bottom-right (reading order), converting every
1768 # 2-by-2 block we see to a single complex element.
1770 for k
in range(n
/2):
1771 for j
in range(n
/2):
1772 submat
= M
[2*k
:2*k
+2,2*j
:2*j
+2]
1773 if submat
[0,0] != submat
[1,1]:
1774 raise ValueError('bad on-diagonal submatrix')
1775 if submat
[0,1] != -submat
[1,0]:
1776 raise ValueError('bad off-diagonal submatrix')
1777 z
= submat
[0,0] + submat
[0,1]*i
1780 return matrix(F
, n
/2, elements
)
1784 def matrix_inner_product(cls
,X
,Y
):
1786 Compute a matrix inner product in this algebra directly from
1791 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1795 This gives the same answer as the slow, default method implemented
1796 in :class:`MatrixEuclideanJordanAlgebra`::
1798 sage: set_random_seed()
1799 sage: J = ComplexHermitianEJA.random_instance()
1800 sage: x,y = J.random_elements(2)
1801 sage: Xe = x.to_matrix()
1802 sage: Ye = y.to_matrix()
1803 sage: X = ComplexHermitianEJA.real_unembed(Xe)
1804 sage: Y = ComplexHermitianEJA.real_unembed(Ye)
1805 sage: expected = (X*Y).trace().real()
1806 sage: actual = ComplexHermitianEJA.matrix_inner_product(Xe,Ye)
1807 sage: actual == expected
1811 return RealMatrixEuclideanJordanAlgebra
.matrix_inner_product(X
,Y
)/2
1814 class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra
,
1815 ConcreteEuclideanJordanAlgebra
):
1817 The rank-n simple EJA consisting of complex Hermitian n-by-n
1818 matrices over the real numbers, the usual symmetric Jordan product,
1819 and the real-part-of-trace inner product. It has dimension `n^2` over
1824 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1828 In theory, our "field" can be any subfield of the reals::
1830 sage: ComplexHermitianEJA(2, RDF)
1831 Euclidean Jordan algebra of dimension 4 over Real Double Field
1832 sage: ComplexHermitianEJA(2, RR)
1833 Euclidean Jordan algebra of dimension 4 over Real Field with
1834 53 bits of precision
1838 The dimension of this algebra is `n^2`::
1840 sage: set_random_seed()
1841 sage: n_max = ComplexHermitianEJA._max_random_instance_size()
1842 sage: n = ZZ.random_element(1, n_max)
1843 sage: J = ComplexHermitianEJA(n)
1844 sage: J.dimension() == n^2
1847 The Jordan multiplication is what we think it is::
1849 sage: set_random_seed()
1850 sage: J = ComplexHermitianEJA.random_instance()
1851 sage: x,y = J.random_elements(2)
1852 sage: actual = (x*y).to_matrix()
1853 sage: X = x.to_matrix()
1854 sage: Y = y.to_matrix()
1855 sage: expected = (X*Y + Y*X)/2
1856 sage: actual == expected
1858 sage: J(expected) == x*y
1861 We can change the generator prefix::
1863 sage: ComplexHermitianEJA(2, prefix='z').gens()
1866 We can construct the (trivial) algebra of rank zero::
1868 sage: ComplexHermitianEJA(0)
1869 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1874 def _denormalized_basis(cls
, n
, field
):
1876 Returns a basis for the space of complex Hermitian n-by-n matrices.
1878 Why do we embed these? Basically, because all of numerical linear
1879 algebra assumes that you're working with vectors consisting of `n`
1880 entries from a field and scalars from the same field. There's no way
1881 to tell SageMath that (for example) the vectors contain complex
1882 numbers, while the scalar field is real.
1886 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1890 sage: set_random_seed()
1891 sage: n = ZZ.random_element(1,5)
1892 sage: field = QuadraticField(2, 'sqrt2')
1893 sage: B = ComplexHermitianEJA._denormalized_basis(n, field)
1894 sage: all( M.is_symmetric() for M in B)
1898 R
= PolynomialRing(field
, 'z')
1900 F
= field
.extension(z
**2 + 1, 'I')
1903 # This is like the symmetric case, but we need to be careful:
1905 # * We want conjugate-symmetry, not just symmetry.
1906 # * The diagonal will (as a result) be real.
1910 for j
in range(i
+1):
1911 Eij
= matrix(F
, n
, lambda k
,l
: k
==i
and l
==j
)
1913 Sij
= cls
.real_embed(Eij
)
1916 # The second one has a minus because it's conjugated.
1917 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
1919 Sij_imag
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
1922 # Since we embedded these, we can drop back to the "field" that we
1923 # started with instead of the complex extension "F".
1924 return ( s
.change_ring(field
) for s
in S
)
1927 def __init__(self
, n
, field
=AA
, **kwargs
):
1928 basis
= self
._denormalized
_basis
(n
,field
)
1929 super(ComplexHermitianEJA
,self
).__init
__(field
,
1933 self
.rank
.set_cache(n
)
1936 def _max_random_instance_size():
1937 return 3 # Dimension 9
1940 def random_instance(cls
, field
=AA
, **kwargs
):
1942 Return a random instance of this type of algebra.
1944 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1945 return cls(n
, field
, **kwargs
)
1947 class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1951 Embed the n-by-n quaternion matrix ``M`` into the space of real
1952 matrices of size 4n-by-4n by first sending each quaternion entry `z
1953 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
1954 c+di],[-c + di, a-bi]]`, and then embedding those into a real
1959 sage: from mjo.eja.eja_algebra import \
1960 ....: QuaternionMatrixEuclideanJordanAlgebra
1964 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1965 sage: i,j,k = Q.gens()
1966 sage: x = 1 + 2*i + 3*j + 4*k
1967 sage: M = matrix(Q, 1, [[x]])
1968 sage: QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1974 Embedding is a homomorphism (isomorphism, in fact)::
1976 sage: set_random_seed()
1977 sage: n = ZZ.random_element(2)
1978 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1979 sage: X = random_matrix(Q, n)
1980 sage: Y = random_matrix(Q, n)
1981 sage: Xe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X)
1982 sage: Ye = QuaternionMatrixEuclideanJordanAlgebra.real_embed(Y)
1983 sage: XYe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1988 quaternions
= M
.base_ring()
1991 raise ValueError("the matrix 'M' must be square")
1993 F
= QuadraticField(-1, 'I')
1998 t
= z
.coefficient_tuple()
2003 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
2004 [-c
+ d
*i
, a
- b
*i
]])
2005 realM
= ComplexMatrixEuclideanJordanAlgebra
.real_embed(cplxM
)
2006 blocks
.append(realM
)
2008 # We should have real entries by now, so use the realest field
2009 # we've got for the return value.
2010 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
2015 def real_unembed(M
):
2017 The inverse of _embed_quaternion_matrix().
2021 sage: from mjo.eja.eja_algebra import \
2022 ....: QuaternionMatrixEuclideanJordanAlgebra
2026 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
2027 ....: [-2, 1, -4, 3],
2028 ....: [-3, 4, 1, -2],
2029 ....: [-4, -3, 2, 1]])
2030 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(M)
2031 [1 + 2*i + 3*j + 4*k]
2035 Unembedding is the inverse of embedding::
2037 sage: set_random_seed()
2038 sage: Q = QuaternionAlgebra(QQ, -1, -1)
2039 sage: M = random_matrix(Q, 3)
2040 sage: Me = QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
2041 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
2047 raise ValueError("the matrix 'M' must be square")
2048 if not n
.mod(4).is_zero():
2049 raise ValueError("the matrix 'M' must be a quaternion embedding")
2051 # Use the base ring of the matrix to ensure that its entries can be
2052 # multiplied by elements of the quaternion algebra.
2053 field
= M
.base_ring()
2054 Q
= QuaternionAlgebra(field
,-1,-1)
2057 # Go top-left to bottom-right (reading order), converting every
2058 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
2061 for l
in range(n
/4):
2062 for m
in range(n
/4):
2063 submat
= ComplexMatrixEuclideanJordanAlgebra
.real_unembed(
2064 M
[4*l
:4*l
+4,4*m
:4*m
+4] )
2065 if submat
[0,0] != submat
[1,1].conjugate():
2066 raise ValueError('bad on-diagonal submatrix')
2067 if submat
[0,1] != -submat
[1,0].conjugate():
2068 raise ValueError('bad off-diagonal submatrix')
2069 z
= submat
[0,0].real()
2070 z
+= submat
[0,0].imag()*i
2071 z
+= submat
[0,1].real()*j
2072 z
+= submat
[0,1].imag()*k
2075 return matrix(Q
, n
/4, elements
)
2079 def matrix_inner_product(cls
,X
,Y
):
2081 Compute a matrix inner product in this algebra directly from
2086 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2090 This gives the same answer as the slow, default method implemented
2091 in :class:`MatrixEuclideanJordanAlgebra`::
2093 sage: set_random_seed()
2094 sage: J = QuaternionHermitianEJA.random_instance()
2095 sage: x,y = J.random_elements(2)
2096 sage: Xe = x.to_matrix()
2097 sage: Ye = y.to_matrix()
2098 sage: X = QuaternionHermitianEJA.real_unembed(Xe)
2099 sage: Y = QuaternionHermitianEJA.real_unembed(Ye)
2100 sage: expected = (X*Y).trace().coefficient_tuple()[0]
2101 sage: actual = QuaternionHermitianEJA.matrix_inner_product(Xe,Ye)
2102 sage: actual == expected
2106 return RealMatrixEuclideanJordanAlgebra
.matrix_inner_product(X
,Y
)/4
2109 class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra
,
2110 ConcreteEuclideanJordanAlgebra
):
2112 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
2113 matrices, the usual symmetric Jordan product, and the
2114 real-part-of-trace inner product. It has dimension `2n^2 - n` over
2119 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2123 In theory, our "field" can be any subfield of the reals::
2125 sage: QuaternionHermitianEJA(2, RDF)
2126 Euclidean Jordan algebra of dimension 6 over Real Double Field
2127 sage: QuaternionHermitianEJA(2, RR)
2128 Euclidean Jordan algebra of dimension 6 over Real Field with
2129 53 bits of precision
2133 The dimension of this algebra is `2*n^2 - n`::
2135 sage: set_random_seed()
2136 sage: n_max = QuaternionHermitianEJA._max_random_instance_size()
2137 sage: n = ZZ.random_element(1, n_max)
2138 sage: J = QuaternionHermitianEJA(n)
2139 sage: J.dimension() == 2*(n^2) - n
2142 The Jordan multiplication is what we think it is::
2144 sage: set_random_seed()
2145 sage: J = QuaternionHermitianEJA.random_instance()
2146 sage: x,y = J.random_elements(2)
2147 sage: actual = (x*y).to_matrix()
2148 sage: X = x.to_matrix()
2149 sage: Y = y.to_matrix()
2150 sage: expected = (X*Y + Y*X)/2
2151 sage: actual == expected
2153 sage: J(expected) == x*y
2156 We can change the generator prefix::
2158 sage: QuaternionHermitianEJA(2, prefix='a').gens()
2159 (a0, a1, a2, a3, a4, a5)
2161 We can construct the (trivial) algebra of rank zero::
2163 sage: QuaternionHermitianEJA(0)
2164 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2168 def _denormalized_basis(cls
, n
, field
):
2170 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
2172 Why do we embed these? Basically, because all of numerical
2173 linear algebra assumes that you're working with vectors consisting
2174 of `n` entries from a field and scalars from the same field. There's
2175 no way to tell SageMath that (for example) the vectors contain
2176 complex numbers, while the scalar field is real.
2180 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2184 sage: set_random_seed()
2185 sage: n = ZZ.random_element(1,5)
2186 sage: B = QuaternionHermitianEJA._denormalized_basis(n,QQ)
2187 sage: all( M.is_symmetric() for M in B )
2191 Q
= QuaternionAlgebra(QQ
,-1,-1)
2194 # This is like the symmetric case, but we need to be careful:
2196 # * We want conjugate-symmetry, not just symmetry.
2197 # * The diagonal will (as a result) be real.
2201 for j
in range(i
+1):
2202 Eij
= matrix(Q
, n
, lambda k
,l
: k
==i
and l
==j
)
2204 Sij
= cls
.real_embed(Eij
)
2207 # The second, third, and fourth ones have a minus
2208 # because they're conjugated.
2209 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
2211 Sij_I
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
2213 Sij_J
= cls
.real_embed(J
*Eij
- J
*Eij
.transpose())
2215 Sij_K
= cls
.real_embed(K
*Eij
- K
*Eij
.transpose())
2218 # Since we embedded these, we can drop back to the "field" that we
2219 # started with instead of the quaternion algebra "Q".
2220 return ( s
.change_ring(field
) for s
in S
)
2223 def __init__(self
, n
, field
=AA
, **kwargs
):
2224 basis
= self
._denormalized
_basis
(n
,field
)
2225 super(QuaternionHermitianEJA
,self
).__init
__(field
,
2229 self
.rank
.set_cache(n
)
2232 def _max_random_instance_size():
2234 The maximum rank of a random QuaternionHermitianEJA.
2236 return 2 # Dimension 6
2239 def random_instance(cls
, field
=AA
, **kwargs
):
2241 Return a random instance of this type of algebra.
2243 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2244 return cls(n
, field
, **kwargs
)
2247 class HadamardEJA(RationalBasisEuclideanJordanAlgebraNg
,
2248 ConcreteEuclideanJordanAlgebra
):
2250 Return the Euclidean Jordan Algebra corresponding to the set
2251 `R^n` under the Hadamard product.
2253 Note: this is nothing more than the Cartesian product of ``n``
2254 copies of the spin algebra. Once Cartesian product algebras
2255 are implemented, this can go.
2259 sage: from mjo.eja.eja_algebra import HadamardEJA
2263 This multiplication table can be verified by hand::
2265 sage: J = HadamardEJA(3)
2266 sage: e0,e1,e2 = J.gens()
2282 We can change the generator prefix::
2284 sage: HadamardEJA(3, prefix='r').gens()
2288 def __init__(self
, n
, field
=AA
, **kwargs
):
2289 V
= VectorSpace(field
, n
)
2292 def jordan_product(x
,y
):
2293 return V([ xi
*yi
for (xi
,yi
) in zip(x
,y
) ])
2294 def inner_product(x
,y
):
2295 return x
.inner_product(y
)
2297 super(HadamardEJA
, self
).__init
__(field
,
2302 self
.rank
.set_cache(n
)
2305 self
.one
.set_cache( self
.zero() )
2307 self
.one
.set_cache( sum(self
.gens()) )
2310 def _max_random_instance_size():
2312 The maximum dimension of a random HadamardEJA.
2317 def random_instance(cls
, field
=AA
, **kwargs
):
2319 Return a random instance of this type of algebra.
2321 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2322 return cls(n
, field
, **kwargs
)
2325 class BilinearFormEJA(RationalBasisEuclideanJordanAlgebraNg
,
2326 ConcreteEuclideanJordanAlgebra
):
2328 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2329 with the half-trace inner product and jordan product ``x*y =
2330 (<Bx,y>,y_bar>, x0*y_bar + y0*x_bar)`` where `B = 1 \times B22` is
2331 a symmetric positive-definite "bilinear form" matrix. Its
2332 dimension is the size of `B`, and it has rank two in dimensions
2333 larger than two. It reduces to the ``JordanSpinEJA`` when `B` is
2334 the identity matrix of order ``n``.
2336 We insist that the one-by-one upper-left identity block of `B` be
2337 passed in as well so that we can be passed a matrix of size zero
2338 to construct a trivial algebra.
2342 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
2343 ....: JordanSpinEJA)
2347 When no bilinear form is specified, the identity matrix is used,
2348 and the resulting algebra is the Jordan spin algebra::
2350 sage: B = matrix.identity(AA,3)
2351 sage: J0 = BilinearFormEJA(B)
2352 sage: J1 = JordanSpinEJA(3)
2353 sage: J0.multiplication_table() == J0.multiplication_table()
2356 An error is raised if the matrix `B` does not correspond to a
2357 positive-definite bilinear form::
2359 sage: B = matrix.random(QQ,2,3)
2360 sage: J = BilinearFormEJA(B)
2361 Traceback (most recent call last):
2363 ValueError: bilinear form is not positive-definite
2364 sage: B = matrix.zero(QQ,3)
2365 sage: J = BilinearFormEJA(B)
2366 Traceback (most recent call last):
2368 ValueError: bilinear form is not positive-definite
2372 We can create a zero-dimensional algebra::
2374 sage: B = matrix.identity(AA,0)
2375 sage: J = BilinearFormEJA(B)
2379 We can check the multiplication condition given in the Jordan, von
2380 Neumann, and Wigner paper (and also discussed on my "On the
2381 symmetry..." paper). Note that this relies heavily on the standard
2382 choice of basis, as does anything utilizing the bilinear form
2383 matrix. We opt not to orthonormalize the basis, because if we
2384 did, we would have to normalize the `s_{i}` in a similar manner::
2386 sage: set_random_seed()
2387 sage: n = ZZ.random_element(5)
2388 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2389 sage: B11 = matrix.identity(QQ,1)
2390 sage: B22 = M.transpose()*M
2391 sage: B = block_matrix(2,2,[ [B11,0 ],
2393 sage: J = BilinearFormEJA(B, orthonormalize=False)
2394 sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
2395 sage: V = J.vector_space()
2396 sage: sis = [ J( V([0] + (M.inverse()*ei).list()).column() )
2397 ....: for ei in eis ]
2398 sage: actual = [ sis[i]*sis[j]
2399 ....: for i in range(n-1)
2400 ....: for j in range(n-1) ]
2401 sage: expected = [ J.one() if i == j else J.zero()
2402 ....: for i in range(n-1)
2403 ....: for j in range(n-1) ]
2404 sage: actual == expected
2407 def __init__(self
, B
, field
=AA
, **kwargs
):
2408 if not B
.is_positive_definite():
2409 raise ValueError("bilinear form is not positive-definite")
2412 V
= VectorSpace(field
, n
)
2414 def inner_product(x
,y
):
2415 return (B
*x
).inner_product(y
)
2417 def jordan_product(x
,y
):
2422 z0
= inner_product(x
,y
)
2423 zbar
= y0
*xbar
+ x0
*ybar
2424 return V([z0
] + zbar
.list())
2426 super(BilinearFormEJA
, self
).__init
__(field
,
2432 # The rank of this algebra is two, unless we're in a
2433 # one-dimensional ambient space (because the rank is bounded
2434 # by the ambient dimension).
2435 self
.rank
.set_cache(min(n
,2))
2438 self
.one
.set_cache( self
.zero() )
2440 self
.one
.set_cache( self
.monomial(0) )
2443 def _max_random_instance_size():
2445 The maximum dimension of a random BilinearFormEJA.
2450 def random_instance(cls
, field
=AA
, **kwargs
):
2452 Return a random instance of this algebra.
2454 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2456 B
= matrix
.identity(field
, n
)
2457 return cls(B
, field
, **kwargs
)
2459 B11
= matrix
.identity(field
,1)
2460 M
= matrix
.random(field
, n
-1)
2461 I
= matrix
.identity(field
, n
-1)
2462 alpha
= field
.zero()
2463 while alpha
.is_zero():
2464 alpha
= field
.random_element().abs()
2465 B22
= M
.transpose()*M
+ alpha
*I
2467 from sage
.matrix
.special
import block_matrix
2468 B
= block_matrix(2,2, [ [B11
, ZZ(0) ],
2471 return cls(B
, field
, **kwargs
)
2474 class JordanSpinEJA(BilinearFormEJA
):
2476 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2477 with the usual inner product and jordan product ``x*y =
2478 (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2483 sage: from mjo.eja.eja_algebra import JordanSpinEJA
2487 This multiplication table can be verified by hand::
2489 sage: J = JordanSpinEJA(4)
2490 sage: e0,e1,e2,e3 = J.gens()
2506 We can change the generator prefix::
2508 sage: JordanSpinEJA(2, prefix='B').gens()
2513 Ensure that we have the usual inner product on `R^n`::
2515 sage: set_random_seed()
2516 sage: J = JordanSpinEJA.random_instance()
2517 sage: x,y = J.random_elements(2)
2518 sage: actual = x.inner_product(y)
2519 sage: expected = x.to_vector().inner_product(y.to_vector())
2520 sage: actual == expected
2524 def __init__(self
, n
, field
=AA
, **kwargs
):
2525 # This is a special case of the BilinearFormEJA with the identity
2526 # matrix as its bilinear form.
2527 B
= matrix
.identity(field
, n
)
2528 super(JordanSpinEJA
, self
).__init
__(B
, field
, **kwargs
)
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(FiniteDimensionalEuclideanJordanAlgebra
,
2549 ConcreteEuclideanJordanAlgebra
):
2551 The trivial Euclidean Jordan algebra consisting of only a zero element.
2555 sage: from mjo.eja.eja_algebra import TrivialEJA
2559 sage: J = TrivialEJA()
2566 sage: 7*J.one()*12*J.one()
2568 sage: J.one().inner_product(J.one())
2570 sage: J.one().norm()
2572 sage: J.one().subalgebra_generated_by()
2573 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2578 def __init__(self
, field
=AA
, **kwargs
):
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)
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(n
) ]
2650 p
= (J1
.monomial(i
)*J1
.monomial(j
)).to_vector()
2651 mult_table
[i
][j
] = V(p
.list() + [field
.zero()]*n2
)
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).
2661 super(DirectSumEJA
, self
).__init
__(field
,
2666 self
.rank
.set_cache(J1
.rank() + J2
.rank())
2671 Return the pair of this algebra's factors.
2675 sage: from mjo.eja.eja_algebra import (HadamardEJA,
2676 ....: JordanSpinEJA,
2681 sage: J1 = HadamardEJA(2,QQ)
2682 sage: J2 = JordanSpinEJA(3,QQ)
2683 sage: J = DirectSumEJA(J1,J2)
2685 (Euclidean Jordan algebra of dimension 2 over Rational Field,
2686 Euclidean Jordan algebra of dimension 3 over Rational Field)
2689 return self
._factors
2691 def projections(self
):
2693 Return a pair of projections onto this algebra's factors.
2697 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
2698 ....: ComplexHermitianEJA,
2703 sage: J1 = JordanSpinEJA(2)
2704 sage: J2 = ComplexHermitianEJA(2)
2705 sage: J = DirectSumEJA(J1,J2)
2706 sage: (pi_left, pi_right) = J.projections()
2707 sage: J.one().to_vector()
2709 sage: pi_left(J.one()).to_vector()
2711 sage: pi_right(J.one()).to_vector()
2715 (J1
,J2
) = self
.factors()
2718 V_basis
= self
.vector_space().basis()
2719 # Need to specify the dimensions explicitly so that we don't
2720 # wind up with a zero-by-zero matrix when we want e.g. a
2721 # zero-by-two matrix (important for composing things).
2722 P1
= matrix(self
.base_ring(), m
, m
+n
, V_basis
[:m
])
2723 P2
= matrix(self
.base_ring(), n
, m
+n
, V_basis
[m
:])
2724 pi_left
= FiniteDimensionalEuclideanJordanAlgebraOperator(self
,J1
,P1
)
2725 pi_right
= FiniteDimensionalEuclideanJordanAlgebraOperator(self
,J2
,P2
)
2726 return (pi_left
, pi_right
)
2728 def inclusions(self
):
2730 Return the pair of inclusion maps from our factors into us.
2734 sage: from mjo.eja.eja_algebra import (random_eja,
2735 ....: JordanSpinEJA,
2736 ....: RealSymmetricEJA,
2741 sage: J1 = JordanSpinEJA(3)
2742 sage: J2 = RealSymmetricEJA(2)
2743 sage: J = DirectSumEJA(J1,J2)
2744 sage: (iota_left, iota_right) = J.inclusions()
2745 sage: iota_left(J1.zero()) == J.zero()
2747 sage: iota_right(J2.zero()) == J.zero()
2749 sage: J1.one().to_vector()
2751 sage: iota_left(J1.one()).to_vector()
2753 sage: J2.one().to_vector()
2755 sage: iota_right(J2.one()).to_vector()
2757 sage: J.one().to_vector()
2762 Composing a projection with the corresponding inclusion should
2763 produce the identity map, and mismatching them should produce
2766 sage: set_random_seed()
2767 sage: J1 = random_eja()
2768 sage: J2 = random_eja()
2769 sage: J = DirectSumEJA(J1,J2)
2770 sage: (iota_left, iota_right) = J.inclusions()
2771 sage: (pi_left, pi_right) = J.projections()
2772 sage: pi_left*iota_left == J1.one().operator()
2774 sage: pi_right*iota_right == J2.one().operator()
2776 sage: (pi_left*iota_right).is_zero()
2778 sage: (pi_right*iota_left).is_zero()
2782 (J1
,J2
) = self
.factors()
2785 V_basis
= self
.vector_space().basis()
2786 # Need to specify the dimensions explicitly so that we don't
2787 # wind up with a zero-by-zero matrix when we want e.g. a
2788 # two-by-zero matrix (important for composing things).
2789 I1
= matrix
.column(self
.base_ring(), m
, m
+n
, V_basis
[:m
])
2790 I2
= matrix
.column(self
.base_ring(), n
, m
+n
, V_basis
[m
:])
2791 iota_left
= FiniteDimensionalEuclideanJordanAlgebraOperator(J1
,self
,I1
)
2792 iota_right
= FiniteDimensionalEuclideanJordanAlgebraOperator(J2
,self
,I2
)
2793 return (iota_left
, iota_right
)
2795 def inner_product(self
, x
, y
):
2797 The standard Cartesian inner-product.
2799 We project ``x`` and ``y`` onto our factors, and add up the
2800 inner-products from the subalgebras.
2805 sage: from mjo.eja.eja_algebra import (HadamardEJA,
2806 ....: QuaternionHermitianEJA,
2811 sage: J1 = HadamardEJA(3,QQ)
2812 sage: J2 = QuaternionHermitianEJA(2,QQ,normalize_basis=False)
2813 sage: J = DirectSumEJA(J1,J2)
2818 sage: x1.inner_product(x2)
2820 sage: y1.inner_product(y2)
2822 sage: J.one().inner_product(J.one())
2826 (pi_left
, pi_right
) = self
.projections()
2832 return (x1
.inner_product(y1
) + x2
.inner_product(y2
))
2836 random_eja
= ConcreteEuclideanJordanAlgebra
.random_instance