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.
8 from itertools
import repeat
10 from sage
.algebras
.quatalg
.quaternion_algebra
import QuaternionAlgebra
11 from sage
.categories
.magmatic_algebras
import MagmaticAlgebras
12 from sage
.combinat
.free_module
import CombinatorialFreeModule
13 from sage
.matrix
.constructor
import matrix
14 from sage
.matrix
.matrix_space
import MatrixSpace
15 from sage
.misc
.cachefunc
import cached_method
16 from sage
.misc
.lazy_import
import lazy_import
17 from sage
.misc
.prandom
import choice
18 from sage
.misc
.table
import table
19 from sage
.modules
.free_module
import FreeModule
, VectorSpace
20 from sage
.rings
.all
import (ZZ
, QQ
, AA
, QQbar
, RR
, RLF
, CLF
,
23 from mjo
.eja
.eja_element
import FiniteDimensionalEuclideanJordanAlgebraElement
24 lazy_import('mjo.eja.eja_subalgebra',
25 'FiniteDimensionalEuclideanJordanSubalgebra')
26 from mjo
.eja
.eja_utils
import _mat2vec
28 class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule
):
30 def _coerce_map_from_base_ring(self
):
32 Disable the map from the base ring into the algebra.
34 Performing a nonsense conversion like this automatically
35 is counterpedagogical. The fallback is to try the usual
36 element constructor, which should also fail.
40 sage: from mjo.eja.eja_algebra import random_eja
44 sage: set_random_seed()
45 sage: J = random_eja()
47 Traceback (most recent call last):
49 ValueError: not a naturally-represented algebra element
65 sage: from mjo.eja.eja_algebra import (
66 ....: FiniteDimensionalEuclideanJordanAlgebra,
72 By definition, Jordan multiplication commutes::
74 sage: set_random_seed()
75 sage: J = random_eja()
76 sage: x,y = J.random_elements(2)
82 The ``field`` we're given must be real with ``check_field=True``::
84 sage: JordanSpinEJA(2,QQbar)
85 Traceback (most recent call last):
87 ValueError: scalar field is not real
89 The multiplication table must be square with ``check_axioms=True``::
91 sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,((),()))
92 Traceback (most recent call last):
94 ValueError: multiplication table is not square
98 if not field
.is_subring(RR
):
99 # Note: this does return true for the real algebraic
100 # field, the rationals, and any quadratic field where
101 # we've specified a real embedding.
102 raise ValueError("scalar field is not real")
104 # The multiplication table had better be square
107 if not all( len(l
) == n
for l
in mult_table
):
108 raise ValueError("multiplication table is not square")
110 self
._natural
_basis
= natural_basis
113 category
= MagmaticAlgebras(field
).FiniteDimensional()
114 category
= category
.WithBasis().Unital()
116 fda
= super(FiniteDimensionalEuclideanJordanAlgebra
, self
)
121 self
.print_options(bracket
='')
123 # The multiplication table we're given is necessarily in terms
124 # of vectors, because we don't have an algebra yet for
125 # anything to be an element of. However, it's faster in the
126 # long run to have the multiplication table be in terms of
127 # algebra elements. We do this after calling the superclass
128 # constructor so that from_vector() knows what to do.
129 self
._multiplication
_table
= [
130 list(map(lambda x
: self
.from_vector(x
), ls
))
135 if not self
._is
_commutative
():
136 raise ValueError("algebra is not commutative")
137 if not self
._is
_jordanian
():
138 raise ValueError("Jordan identity does not hold")
139 if not self
._inner
_product
_is
_associative
():
140 raise ValueError("inner product is not associative")
142 def _element_constructor_(self
, elt
):
144 Construct an element of this algebra from its natural
147 This gets called only after the parent element _call_ method
148 fails to find a coercion for the argument.
152 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
154 ....: RealSymmetricEJA)
158 The identity in `S^n` is converted to the identity in the EJA::
160 sage: J = RealSymmetricEJA(3)
161 sage: I = matrix.identity(QQ,3)
162 sage: J(I) == J.one()
165 This skew-symmetric matrix can't be represented in the EJA::
167 sage: J = RealSymmetricEJA(3)
168 sage: A = matrix(QQ,3, lambda i,j: i-j)
170 Traceback (most recent call last):
172 ArithmeticError: vector is not in free module
176 Ensure that we can convert any element of the two non-matrix
177 simple algebras (whose natural representations are their usual
178 vector representations) back and forth faithfully::
180 sage: set_random_seed()
181 sage: J = HadamardEJA.random_instance()
182 sage: x = J.random_element()
183 sage: J(x.to_vector().column()) == x
185 sage: J = JordanSpinEJA.random_instance()
186 sage: x = J.random_element()
187 sage: J(x.to_vector().column()) == x
191 msg
= "not a naturally-represented algebra element"
193 # The superclass implementation of random_element()
194 # needs to be able to coerce "0" into the algebra.
196 elif elt
in self
.base_ring():
197 # Ensure that no base ring -> algebra coercion is performed
198 # by this method. There's some stupidity in sage that would
199 # otherwise propagate to this method; for example, sage thinks
200 # that the integer 3 belongs to the space of 2-by-2 matrices.
201 raise ValueError(msg
)
203 natural_basis
= self
.natural_basis()
204 basis_space
= natural_basis
[0].matrix_space()
205 if elt
not in basis_space
:
206 raise ValueError(msg
)
208 # Thanks for nothing! Matrix spaces aren't vector spaces in
209 # Sage, so we have to figure out its natural-basis coordinates
210 # ourselves. We use the basis space's ring instead of the
211 # element's ring because the basis space might be an algebraic
212 # closure whereas the base ring of the 3-by-3 identity matrix
213 # could be QQ instead of QQbar.
214 V
= VectorSpace(basis_space
.base_ring(), elt
.nrows()*elt
.ncols())
215 W
= V
.span_of_basis( _mat2vec(s
) for s
in natural_basis
)
216 coords
= W
.coordinate_vector(_mat2vec(elt
))
217 return self
.from_vector(coords
)
220 def _max_random_instance_size():
222 Return an integer "size" that is an upper bound on the size of
223 this algebra when it is used in a random test
224 case. Unfortunately, the term "size" is quite vague -- when
225 dealing with `R^n` under either the Hadamard or Jordan spin
226 product, the "size" refers to the dimension `n`. When dealing
227 with a matrix algebra (real symmetric or complex/quaternion
228 Hermitian), it refers to the size of the matrix, which is
229 far less than the dimension of the underlying vector space.
231 We default to five in this class, which is safe in `R^n`. The
232 matrix algebra subclasses (or any class where the "size" is
233 interpreted to be far less than the dimension) should override
234 with a smaller number.
236 raise NotImplementedError
240 Return a string representation of ``self``.
244 sage: from mjo.eja.eja_algebra import JordanSpinEJA
248 Ensure that it says what we think it says::
250 sage: JordanSpinEJA(2, field=AA)
251 Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
252 sage: JordanSpinEJA(3, field=RDF)
253 Euclidean Jordan algebra of dimension 3 over Real Double Field
256 fmt
= "Euclidean Jordan algebra of dimension {} over {}"
257 return fmt
.format(self
.dimension(), self
.base_ring())
259 def product_on_basis(self
, i
, j
):
260 return self
._multiplication
_table
[i
][j
]
262 def _is_commutative(self
):
264 Whether or not this algebra's multiplication table is commutative.
266 This method should of course always return ``True``, unless
267 this algebra was constructed with ``check_axioms=False`` and
268 passed an invalid multiplication table.
270 return all( self
.product_on_basis(i
,j
) == self
.product_on_basis(i
,j
)
271 for i
in range(self
.dimension())
272 for j
in range(self
.dimension()) )
274 def _is_jordanian(self
):
276 Whether or not this algebra's multiplication table respects the
277 Jordan identity `(x^{2})(xy) = x(x^{2}y)`.
279 We only check one arrangement of `x` and `y`, so for a
280 ``True`` result to be truly true, you should also check
281 :meth:`_is_commutative`. This method should of course always
282 return ``True``, unless this algebra was constructed with
283 ``check_axioms=False`` and passed an invalid multiplication table.
285 return all( (self
.monomial(i
)**2)*(self
.monomial(i
)*self
.monomial(j
))
287 (self
.monomial(i
))*((self
.monomial(i
)**2)*self
.monomial(j
))
288 for i
in range(self
.dimension())
289 for j
in range(self
.dimension()) )
291 def _inner_product_is_associative(self
):
293 Return whether or not this algebra's inner product `B` is
294 associative; that is, whether or not `B(xy,z) = B(x,yz)`.
296 This method should of course always return ``True``, unless
297 this algebra was constructed with ``check_axioms=False`` and
298 passed an invalid multiplication table.
301 # Used to check whether or not something is zero in an inexact
302 # ring. This number is sufficient to allow the construction of
303 # QuaternionHermitianEJA(2, RDF) with check_axioms=True.
306 for i
in range(self
.dimension()):
307 for j
in range(self
.dimension()):
308 for k
in range(self
.dimension()):
312 diff
= (x
*y
).inner_product(z
) - x
.inner_product(y
*z
)
314 if self
.base_ring().is_exact():
318 if diff
.abs() > epsilon
:
324 def characteristic_polynomial_of(self
):
326 Return the algebra's "characteristic polynomial of" function,
327 which is itself a multivariate polynomial that, when evaluated
328 at the coordinates of some algebra element, returns that
329 element's characteristic polynomial.
331 The resulting polynomial has `n+1` variables, where `n` is the
332 dimension of this algebra. The first `n` variables correspond to
333 the coordinates of an algebra element: when evaluated at the
334 coordinates of an algebra element with respect to a certain
335 basis, the result is a univariate polynomial (in the one
336 remaining variable ``t``), namely the characteristic polynomial
341 sage: from mjo.eja.eja_algebra import JordanSpinEJA, TrivialEJA
345 The characteristic polynomial in the spin algebra is given in
346 Alizadeh, Example 11.11::
348 sage: J = JordanSpinEJA(3)
349 sage: p = J.characteristic_polynomial_of(); p
350 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
351 sage: xvec = J.one().to_vector()
355 By definition, the characteristic polynomial is a monic
356 degree-zero polynomial in a rank-zero algebra. Note that
357 Cayley-Hamilton is indeed satisfied since the polynomial
358 ``1`` evaluates to the identity element of the algebra on
361 sage: J = TrivialEJA()
362 sage: J.characteristic_polynomial_of()
369 # The list of coefficient polynomials a_0, a_1, a_2, ..., a_(r-1).
370 a
= self
._charpoly
_coefficients
()
372 # We go to a bit of trouble here to reorder the
373 # indeterminates, so that it's easier to evaluate the
374 # characteristic polynomial at x's coordinates and get back
375 # something in terms of t, which is what we want.
376 S
= PolynomialRing(self
.base_ring(),'t')
380 S
= PolynomialRing(S
, R
.variable_names())
383 return (t
**r
+ sum( a
[k
]*(t
**k
) for k
in range(r
) ))
386 def inner_product(self
, x
, y
):
388 The inner product associated with this Euclidean Jordan algebra.
390 Defaults to the trace inner product, but can be overridden by
391 subclasses if they are sure that the necessary properties are
396 sage: from mjo.eja.eja_algebra import random_eja
400 Our inner product is "associative," which means the following for
401 a symmetric bilinear form::
403 sage: set_random_seed()
404 sage: J = random_eja()
405 sage: x,y,z = J.random_elements(3)
406 sage: (x*y).inner_product(z) == y.inner_product(x*z)
410 X
= x
.natural_representation()
411 Y
= y
.natural_representation()
412 return self
.natural_inner_product(X
,Y
)
415 def is_trivial(self
):
417 Return whether or not this algebra is trivial.
419 A trivial algebra contains only the zero element.
423 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
428 sage: J = ComplexHermitianEJA(3)
434 sage: J = TrivialEJA()
439 return self
.dimension() == 0
442 def multiplication_table(self
):
444 Return a visual representation of this algebra's multiplication
445 table (on basis elements).
449 sage: from mjo.eja.eja_algebra import JordanSpinEJA
453 sage: J = JordanSpinEJA(4)
454 sage: J.multiplication_table()
455 +----++----+----+----+----+
456 | * || e0 | e1 | e2 | e3 |
457 +====++====+====+====+====+
458 | e0 || e0 | e1 | e2 | e3 |
459 +----++----+----+----+----+
460 | e1 || e1 | e0 | 0 | 0 |
461 +----++----+----+----+----+
462 | e2 || e2 | 0 | e0 | 0 |
463 +----++----+----+----+----+
464 | e3 || e3 | 0 | 0 | e0 |
465 +----++----+----+----+----+
468 M
= list(self
._multiplication
_table
) # copy
469 for i
in range(len(M
)):
470 # M had better be "square"
471 M
[i
] = [self
.monomial(i
)] + M
[i
]
472 M
= [["*"] + list(self
.gens())] + M
473 return table(M
, header_row
=True, header_column
=True, frame
=True)
476 def natural_basis(self
):
478 Return a more-natural representation of this algebra's basis.
480 Every finite-dimensional Euclidean Jordan Algebra is a direct
481 sum of five simple algebras, four of which comprise Hermitian
482 matrices. This method returns the original "natural" basis
483 for our underlying vector space. (Typically, the natural basis
484 is used to construct the multiplication table in the first place.)
486 Note that this will always return a matrix. The standard basis
487 in `R^n` will be returned as `n`-by-`1` column matrices.
491 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
492 ....: RealSymmetricEJA)
496 sage: J = RealSymmetricEJA(2)
498 Finite family {0: e0, 1: e1, 2: e2}
499 sage: J.natural_basis()
501 [1 0] [ 0 0.7071067811865475?] [0 0]
502 [0 0], [0.7071067811865475? 0], [0 1]
507 sage: J = JordanSpinEJA(2)
509 Finite family {0: e0, 1: e1}
510 sage: J.natural_basis()
517 if self
._natural
_basis
is None:
518 M
= self
.natural_basis_space()
519 return tuple( M(b
.to_vector()) for b
in self
.basis() )
521 return self
._natural
_basis
524 def natural_basis_space(self
):
526 Return the matrix space in which this algebra's natural basis
529 Generally this will be an `n`-by-`1` column-vector space,
530 except when the algebra is trivial. There it's `n`-by-`n`
531 (where `n` is zero), to ensure that two elements of the
532 natural basis space (empty matrices) can be multiplied.
534 if self
.is_trivial():
535 return MatrixSpace(self
.base_ring(), 0)
536 elif self
._natural
_basis
is None or len(self
._natural
_basis
) == 0:
537 return MatrixSpace(self
.base_ring(), self
.dimension(), 1)
539 return self
._natural
_basis
[0].matrix_space()
543 def natural_inner_product(X
,Y
):
545 Compute the inner product of two naturally-represented elements.
547 For example in the real symmetric matrix EJA, this will compute
548 the trace inner-product of two n-by-n symmetric matrices. The
549 default should work for the real cartesian product EJA, the
550 Jordan spin EJA, and the real symmetric matrices. The others
551 will have to be overridden.
553 return (X
.conjugate_transpose()*Y
).trace()
559 Return the unit element of this algebra.
563 sage: from mjo.eja.eja_algebra import (HadamardEJA,
568 sage: J = HadamardEJA(5)
570 e0 + e1 + e2 + e3 + e4
574 The identity element acts like the identity::
576 sage: set_random_seed()
577 sage: J = random_eja()
578 sage: x = J.random_element()
579 sage: J.one()*x == x and x*J.one() == x
582 The matrix of the unit element's operator is the identity::
584 sage: set_random_seed()
585 sage: J = random_eja()
586 sage: actual = J.one().operator().matrix()
587 sage: expected = matrix.identity(J.base_ring(), J.dimension())
588 sage: actual == expected
592 # We can brute-force compute the matrices of the operators
593 # that correspond to the basis elements of this algebra.
594 # If some linear combination of those basis elements is the
595 # algebra identity, then the same linear combination of
596 # their matrices has to be the identity matrix.
598 # Of course, matrices aren't vectors in sage, so we have to
599 # appeal to the "long vectors" isometry.
600 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
602 # Now we use basic linear algebra to find the coefficients,
603 # of the matrices-as-vectors-linear-combination, which should
604 # work for the original algebra basis too.
605 A
= matrix(self
.base_ring(), oper_vecs
)
607 # We used the isometry on the left-hand side already, but we
608 # still need to do it for the right-hand side. Recall that we
609 # wanted something that summed to the identity matrix.
610 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
612 # Now if there's an identity element in the algebra, this
613 # should work. We solve on the left to avoid having to
614 # transpose the matrix "A".
615 return self
.from_vector(A
.solve_left(b
))
618 def peirce_decomposition(self
, c
):
620 The Peirce decomposition of this algebra relative to the
623 In the future, this can be extended to a complete system of
624 orthogonal idempotents.
628 - ``c`` -- an idempotent of this algebra.
632 A triple (J0, J5, J1) containing two subalgebras and one subspace
635 - ``J0`` -- the algebra on the eigenspace of ``c.operator()``
636 corresponding to the eigenvalue zero.
638 - ``J5`` -- the eigenspace (NOT a subalgebra) of ``c.operator()``
639 corresponding to the eigenvalue one-half.
641 - ``J1`` -- the algebra on the eigenspace of ``c.operator()``
642 corresponding to the eigenvalue one.
644 These are the only possible eigenspaces for that operator, and this
645 algebra is a direct sum of them. The spaces ``J0`` and ``J1`` are
646 orthogonal, and are subalgebras of this algebra with the appropriate
651 sage: from mjo.eja.eja_algebra import random_eja, RealSymmetricEJA
655 The canonical example comes from the symmetric matrices, which
656 decompose into diagonal and off-diagonal parts::
658 sage: J = RealSymmetricEJA(3)
659 sage: C = matrix(QQ, [ [1,0,0],
663 sage: J0,J5,J1 = J.peirce_decomposition(c)
665 Euclidean Jordan algebra of dimension 1...
667 Vector space of degree 6 and dimension 2...
669 Euclidean Jordan algebra of dimension 3...
670 sage: J0.one().natural_representation()
674 sage: orig_df = AA.options.display_format
675 sage: AA.options.display_format = 'radical'
676 sage: J.from_vector(J5.basis()[0]).natural_representation()
680 sage: J.from_vector(J5.basis()[1]).natural_representation()
684 sage: AA.options.display_format = orig_df
685 sage: J1.one().natural_representation()
692 Every algebra decomposes trivially with respect to its identity
695 sage: set_random_seed()
696 sage: J = random_eja()
697 sage: J0,J5,J1 = J.peirce_decomposition(J.one())
698 sage: J0.dimension() == 0 and J5.dimension() == 0
700 sage: J1.superalgebra() == J and J1.dimension() == J.dimension()
703 The decomposition is into eigenspaces, and its components are
704 therefore necessarily orthogonal. Moreover, the identity
705 elements in the two subalgebras are the projections onto their
706 respective subspaces of the superalgebra's identity element::
708 sage: set_random_seed()
709 sage: J = random_eja()
710 sage: x = J.random_element()
711 sage: if not J.is_trivial():
712 ....: while x.is_nilpotent():
713 ....: x = J.random_element()
714 sage: c = x.subalgebra_idempotent()
715 sage: J0,J5,J1 = J.peirce_decomposition(c)
717 sage: for (w,y,z) in zip(J0.basis(), J5.basis(), J1.basis()):
718 ....: w = w.superalgebra_element()
719 ....: y = J.from_vector(y)
720 ....: z = z.superalgebra_element()
721 ....: ipsum += w.inner_product(y).abs()
722 ....: ipsum += w.inner_product(z).abs()
723 ....: ipsum += y.inner_product(z).abs()
726 sage: J1(c) == J1.one()
728 sage: J0(J.one() - c) == J0.one()
732 if not c
.is_idempotent():
733 raise ValueError("element is not idempotent: %s" % c
)
735 # Default these to what they should be if they turn out to be
736 # trivial, because eigenspaces_left() won't return eigenvalues
737 # corresponding to trivial spaces (e.g. it returns only the
738 # eigenspace corresponding to lambda=1 if you take the
739 # decomposition relative to the identity element).
740 trivial
= FiniteDimensionalEuclideanJordanSubalgebra(self
, ())
741 J0
= trivial
# eigenvalue zero
742 J5
= VectorSpace(self
.base_ring(), 0) # eigenvalue one-half
743 J1
= trivial
# eigenvalue one
745 for (eigval
, eigspace
) in c
.operator().matrix().right_eigenspaces():
746 if eigval
== ~
(self
.base_ring()(2)):
749 gens
= tuple( self
.from_vector(b
) for b
in eigspace
.basis() )
750 subalg
= FiniteDimensionalEuclideanJordanSubalgebra(self
,
758 raise ValueError("unexpected eigenvalue: %s" % eigval
)
763 def random_element(self
, thorough
=False):
765 Return a random element of this algebra.
767 Our algebra superclass method only returns a linear
768 combination of at most two basis elements. We instead
769 want the vector space "random element" method that
770 returns a more diverse selection.
774 - ``thorough`` -- (boolean; default False) whether or not we
775 should generate irrational coefficients for the random
776 element when our base ring is irrational; this slows the
777 algebra operations to a crawl, but any truly random method
781 # For a general base ring... maybe we can trust this to do the
782 # right thing? Unlikely, but.
783 V
= self
.vector_space()
784 v
= V
.random_element()
786 if self
.base_ring() is AA
:
787 # The "random element" method of the algebraic reals is
788 # stupid at the moment, and only returns integers between
789 # -2 and 2, inclusive:
791 # https://trac.sagemath.org/ticket/30875
793 # Instead, we implement our own "random vector" method,
794 # and then coerce that into the algebra. We use the vector
795 # space degree here instead of the dimension because a
796 # subalgebra could (for example) be spanned by only two
797 # vectors, each with five coordinates. We need to
798 # generate all five coordinates.
800 v
*= QQbar
.random_element().real()
802 v
*= QQ
.random_element()
804 return self
.from_vector(V
.coordinate_vector(v
))
806 def random_elements(self
, count
, thorough
=False):
808 Return ``count`` random elements as a tuple.
812 - ``thorough`` -- (boolean; default False) whether or not we
813 should generate irrational coefficients for the random
814 elements when our base ring is irrational; this slows the
815 algebra operations to a crawl, but any truly random method
820 sage: from mjo.eja.eja_algebra import JordanSpinEJA
824 sage: J = JordanSpinEJA(3)
825 sage: x,y,z = J.random_elements(3)
826 sage: all( [ x in J, y in J, z in J ])
828 sage: len( J.random_elements(10) ) == 10
832 return tuple( self
.random_element(thorough
)
833 for idx
in range(count
) )
836 def random_instance(cls
, field
=AA
, **kwargs
):
838 Return a random instance of this type of algebra.
840 Beware, this will crash for "most instances" because the
841 constructor below looks wrong.
843 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
844 return cls(n
, field
, **kwargs
)
847 def _charpoly_coefficients(self
):
849 The `r` polynomial coefficients of the "characteristic polynomial
853 var_names
= [ "X" + str(z
) for z
in range(1,n
+1) ]
854 R
= PolynomialRing(self
.base_ring(), var_names
)
856 F
= R
.fraction_field()
859 # From a result in my book, these are the entries of the
860 # basis representation of L_x.
861 return sum( vars[k
]*self
.monomial(k
).operator().matrix()[i
,j
]
864 L_x
= matrix(F
, n
, n
, L_x_i_j
)
867 if self
.rank
.is_in_cache():
869 # There's no need to pad the system with redundant
870 # columns if we *know* they'll be redundant.
873 # Compute an extra power in case the rank is equal to
874 # the dimension (otherwise, we would stop at x^(r-1)).
875 x_powers
= [ (L_x
**k
)*self
.one().to_vector()
876 for k
in range(n
+1) ]
877 A
= matrix
.column(F
, x_powers
[:n
])
878 AE
= A
.extended_echelon_form()
885 # The theory says that only the first "r" coefficients are
886 # nonzero, and they actually live in the original polynomial
887 # ring and not the fraction field. We negate them because
888 # in the actual characteristic polynomial, they get moved
889 # to the other side where x^r lives.
890 return -A_rref
.solve_right(E
*b
).change_ring(R
)[:r
]
895 Return the rank of this EJA.
897 This is a cached method because we know the rank a priori for
898 all of the algebras we can construct. Thus we can avoid the
899 expensive ``_charpoly_coefficients()`` call unless we truly
900 need to compute the whole characteristic polynomial.
904 sage: from mjo.eja.eja_algebra import (HadamardEJA,
906 ....: RealSymmetricEJA,
907 ....: ComplexHermitianEJA,
908 ....: QuaternionHermitianEJA,
913 The rank of the Jordan spin algebra is always two::
915 sage: JordanSpinEJA(2).rank()
917 sage: JordanSpinEJA(3).rank()
919 sage: JordanSpinEJA(4).rank()
922 The rank of the `n`-by-`n` Hermitian real, complex, or
923 quaternion matrices is `n`::
925 sage: RealSymmetricEJA(4).rank()
927 sage: ComplexHermitianEJA(3).rank()
929 sage: QuaternionHermitianEJA(2).rank()
934 Ensure that every EJA that we know how to construct has a
935 positive integer rank, unless the algebra is trivial in
936 which case its rank will be zero::
938 sage: set_random_seed()
939 sage: J = random_eja()
943 sage: r > 0 or (r == 0 and J.is_trivial())
946 Ensure that computing the rank actually works, since the ranks
947 of all simple algebras are known and will be cached by default::
949 sage: J = HadamardEJA(4)
950 sage: J.rank.clear_cache()
956 sage: J = JordanSpinEJA(4)
957 sage: J.rank.clear_cache()
963 sage: J = RealSymmetricEJA(3)
964 sage: J.rank.clear_cache()
970 sage: J = ComplexHermitianEJA(2)
971 sage: J.rank.clear_cache()
977 sage: J = QuaternionHermitianEJA(2)
978 sage: J.rank.clear_cache()
982 return len(self
._charpoly
_coefficients
())
985 def vector_space(self
):
987 Return the vector space that underlies this algebra.
991 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
995 sage: J = RealSymmetricEJA(2)
996 sage: J.vector_space()
997 Vector space of dimension 3 over...
1000 return self
.zero().to_vector().parent().ambient_vector_space()
1003 Element
= FiniteDimensionalEuclideanJordanAlgebraElement
1007 def random_eja(field
=AA
):
1009 Return a "random" finite-dimensional Euclidean Jordan Algebra.
1013 sage: from mjo.eja.eja_algebra import random_eja
1018 Euclidean Jordan algebra of dimension...
1021 classname
= choice([TrivialEJA
,
1025 ComplexHermitianEJA
,
1026 QuaternionHermitianEJA
])
1027 return classname
.random_instance(field
=field
)
1032 class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra
):
1034 Algebras whose basis consists of vectors with rational
1035 entries. Equivalently, algebras whose multiplication tables
1036 contain only rational coefficients.
1038 When an EJA has a basis that can be made rational, we can speed up
1039 the computation of its characteristic polynomial by doing it over
1040 ``QQ``. All of the named EJA constructors that we provide fall
1044 def _charpoly_coefficients(self
):
1046 Override the parent method with something that tries to compute
1047 over a faster (non-extension) field.
1051 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1055 The base ring of the resulting polynomial coefficients is what
1056 it should be, and not the rationals (unless the algebra was
1057 already over the rationals)::
1059 sage: J = JordanSpinEJA(3)
1060 sage: J._charpoly_coefficients()
1061 (X1^2 - X2^2 - X3^2, -2*X1)
1062 sage: a0 = J._charpoly_coefficients()[0]
1064 Algebraic Real Field
1065 sage: a0.base_ring()
1066 Algebraic Real Field
1069 if self
.base_ring() is QQ
:
1070 # There's no need to construct *another* algebra over the
1071 # rationals if this one is already over the rationals.
1072 superclass
= super(RationalBasisEuclideanJordanAlgebra
, self
)
1073 return superclass
._charpoly
_coefficients
()
1076 map(lambda x
: x
.to_vector(), ls
)
1077 for ls
in self
._multiplication
_table
1080 # Do the computation over the rationals. The answer will be
1081 # the same, because our basis coordinates are (essentially)
1083 J
= FiniteDimensionalEuclideanJordanAlgebra(QQ
,
1087 a
= J
._charpoly
_coefficients
()
1088 return tuple(map(lambda x
: x
.change_ring(self
.base_ring()), a
))
1091 class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra
):
1093 def _max_random_instance_size():
1094 # Play it safe, since this will be squared and the underlying
1095 # field can have dimension 4 (quaternions) too.
1098 def __init__(self
, field
, basis
, normalize_basis
=True, **kwargs
):
1100 Compared to the superclass constructor, we take a basis instead of
1101 a multiplication table because the latter can be computed in terms
1102 of the former when the product is known (like it is here).
1104 # Used in this class's fast _charpoly_coefficients() override.
1105 self
._basis
_normalizers
= None
1107 # We're going to loop through this a few times, so now's a good
1108 # time to ensure that it isn't a generator expression.
1109 basis
= tuple(basis
)
1111 if len(basis
) > 1 and normalize_basis
:
1112 # We'll need sqrt(2) to normalize the basis, and this
1113 # winds up in the multiplication table, so the whole
1114 # algebra needs to be over the field extension.
1115 R
= PolynomialRing(field
, 'z')
1118 if p
.is_irreducible():
1119 field
= field
.extension(p
, 'sqrt2', embedding
=RLF(2).sqrt())
1120 basis
= tuple( s
.change_ring(field
) for s
in basis
)
1121 self
._basis
_normalizers
= tuple(
1122 ~
(self
.natural_inner_product(s
,s
).sqrt()) for s
in basis
)
1123 basis
= tuple(s
*c
for (s
,c
) in zip(basis
,self
._basis
_normalizers
))
1125 Qs
= self
.multiplication_table_from_matrix_basis(basis
)
1127 super(MatrixEuclideanJordanAlgebra
, self
).__init
__(field
,
1129 natural_basis
=basis
,
1134 def _charpoly_coefficients(self
):
1136 Override the parent method with something that tries to compute
1137 over a faster (non-extension) field.
1139 if self
._basis
_normalizers
is None or self
.base_ring() is QQ
:
1140 # We didn't normalize, or the basis we started with had
1141 # entries in a nice field already. Just compute the thing.
1142 return super(MatrixEuclideanJordanAlgebra
, self
)._charpoly
_coefficients
()
1144 basis
= ( (b
/n
) for (b
,n
) in zip(self
.natural_basis(),
1145 self
._basis
_normalizers
) )
1147 # Do this over the rationals and convert back at the end.
1148 # Only works because we know the entries of the basis are
1149 # integers. The argument ``check_axioms=False`` is required
1150 # because the trace inner-product method for this
1151 # class is a stub and can't actually be checked.
1152 J
= MatrixEuclideanJordanAlgebra(QQ
,
1154 normalize_basis
=False,
1157 a
= J
._charpoly
_coefficients
()
1159 # Unfortunately, changing the basis does change the
1160 # coefficients of the characteristic polynomial, but since
1161 # these are really the coefficients of the "characteristic
1162 # polynomial of" function, everything is still nice and
1163 # unevaluated. It's therefore "obvious" how scaling the
1164 # basis affects the coordinate variables X1, X2, et
1165 # cetera. Scaling the first basis vector up by "n" adds a
1166 # factor of 1/n into every "X1" term, for example. So here
1167 # we simply undo the basis_normalizer scaling that we
1168 # performed earlier.
1170 # The a[0] access here is safe because trivial algebras
1171 # won't have any basis normalizers and therefore won't
1172 # make it to this "else" branch.
1173 XS
= a
[0].parent().gens()
1174 subs_dict
= { XS
[i
]: self
._basis
_normalizers
[i
]*XS
[i
]
1175 for i
in range(len(XS
)) }
1176 return tuple( a_i
.subs(subs_dict
) for a_i
in a
)
1180 def multiplication_table_from_matrix_basis(basis
):
1182 At least three of the five simple Euclidean Jordan algebras have the
1183 symmetric multiplication (A,B) |-> (AB + BA)/2, where the
1184 multiplication on the right is matrix multiplication. Given a basis
1185 for the underlying matrix space, this function returns a
1186 multiplication table (obtained by looping through the basis
1187 elements) for an algebra of those matrices.
1189 # In S^2, for example, we nominally have four coordinates even
1190 # though the space is of dimension three only. The vector space V
1191 # is supposed to hold the entire long vector, and the subspace W
1192 # of V will be spanned by the vectors that arise from symmetric
1193 # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
1197 field
= basis
[0].base_ring()
1198 dimension
= basis
[0].nrows()
1200 V
= VectorSpace(field
, dimension
**2)
1201 W
= V
.span_of_basis( _mat2vec(s
) for s
in basis
)
1203 mult_table
= [[W
.zero() for j
in range(n
)] for i
in range(n
)]
1206 mat_entry
= (basis
[i
]*basis
[j
] + basis
[j
]*basis
[i
])/2
1207 mult_table
[i
][j
] = W
.coordinate_vector(_mat2vec(mat_entry
))
1215 Embed the matrix ``M`` into a space of real matrices.
1217 The matrix ``M`` can have entries in any field at the moment:
1218 the real numbers, complex numbers, or quaternions. And although
1219 they are not a field, we can probably support octonions at some
1220 point, too. This function returns a real matrix that "acts like"
1221 the original with respect to matrix multiplication; i.e.
1223 real_embed(M*N) = real_embed(M)*real_embed(N)
1226 raise NotImplementedError
1230 def real_unembed(M
):
1232 The inverse of :meth:`real_embed`.
1234 raise NotImplementedError
1238 def natural_inner_product(cls
,X
,Y
):
1239 Xu
= cls
.real_unembed(X
)
1240 Yu
= cls
.real_unembed(Y
)
1241 tr
= (Xu
*Yu
).trace()
1244 # Works in QQ, AA, RDF, et cetera.
1246 except AttributeError:
1247 # A quaternion doesn't have a real() method, but does
1248 # have coefficient_tuple() method that returns the
1249 # coefficients of 1, i, j, and k -- in that order.
1250 return tr
.coefficient_tuple()[0]
1253 class RealMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1257 The identity function, for embedding real matrices into real
1263 def real_unembed(M
):
1265 The identity function, for unembedding real matrices from real
1271 class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra
):
1273 The rank-n simple EJA consisting of real symmetric n-by-n
1274 matrices, the usual symmetric Jordan product, and the trace inner
1275 product. It has dimension `(n^2 + n)/2` over the reals.
1279 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1283 sage: J = RealSymmetricEJA(2)
1284 sage: e0, e1, e2 = J.gens()
1292 In theory, our "field" can be any subfield of the reals::
1294 sage: RealSymmetricEJA(2, RDF)
1295 Euclidean Jordan algebra of dimension 3 over Real Double Field
1296 sage: RealSymmetricEJA(2, RR)
1297 Euclidean Jordan algebra of dimension 3 over Real Field with
1298 53 bits of precision
1302 The dimension of this algebra is `(n^2 + n) / 2`::
1304 sage: set_random_seed()
1305 sage: n_max = RealSymmetricEJA._max_random_instance_size()
1306 sage: n = ZZ.random_element(1, n_max)
1307 sage: J = RealSymmetricEJA(n)
1308 sage: J.dimension() == (n^2 + n)/2
1311 The Jordan multiplication is what we think it is::
1313 sage: set_random_seed()
1314 sage: J = RealSymmetricEJA.random_instance()
1315 sage: x,y = J.random_elements(2)
1316 sage: actual = (x*y).natural_representation()
1317 sage: X = x.natural_representation()
1318 sage: Y = y.natural_representation()
1319 sage: expected = (X*Y + Y*X)/2
1320 sage: actual == expected
1322 sage: J(expected) == x*y
1325 We can change the generator prefix::
1327 sage: RealSymmetricEJA(3, prefix='q').gens()
1328 (q0, q1, q2, q3, q4, q5)
1330 Our natural basis is normalized with respect to the natural inner
1331 product unless we specify otherwise::
1333 sage: set_random_seed()
1334 sage: J = RealSymmetricEJA.random_instance()
1335 sage: all( b.norm() == 1 for b in J.gens() )
1338 Since our natural basis is normalized with respect to the natural
1339 inner product, and since we know that this algebra is an EJA, any
1340 left-multiplication operator's matrix will be symmetric because
1341 natural->EJA basis representation is an isometry and within the EJA
1342 the operator is self-adjoint by the Jordan axiom::
1344 sage: set_random_seed()
1345 sage: x = RealSymmetricEJA.random_instance().random_element()
1346 sage: x.operator().matrix().is_symmetric()
1349 We can construct the (trivial) algebra of rank zero::
1351 sage: RealSymmetricEJA(0)
1352 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1356 def _denormalized_basis(cls
, n
, field
):
1358 Return a basis for the space of real symmetric n-by-n matrices.
1362 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1366 sage: set_random_seed()
1367 sage: n = ZZ.random_element(1,5)
1368 sage: B = RealSymmetricEJA._denormalized_basis(n,QQ)
1369 sage: all( M.is_symmetric() for M in B)
1373 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1377 for j
in range(i
+1):
1378 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
1382 Sij
= Eij
+ Eij
.transpose()
1388 def _max_random_instance_size():
1389 return 4 # Dimension 10
1392 def __init__(self
, n
, field
=AA
, **kwargs
):
1393 basis
= self
._denormalized
_basis
(n
, field
)
1394 super(RealSymmetricEJA
, self
).__init
__(field
,
1398 self
.rank
.set_cache(n
)
1401 class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1405 Embed the n-by-n complex matrix ``M`` into the space of real
1406 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1407 bi` to the block matrix ``[[a,b],[-b,a]]``.
1411 sage: from mjo.eja.eja_algebra import \
1412 ....: ComplexMatrixEuclideanJordanAlgebra
1416 sage: F = QuadraticField(-1, 'I')
1417 sage: x1 = F(4 - 2*i)
1418 sage: x2 = F(1 + 2*i)
1421 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1422 sage: ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1431 Embedding is a homomorphism (isomorphism, in fact)::
1433 sage: set_random_seed()
1434 sage: n_max = ComplexMatrixEuclideanJordanAlgebra._max_random_instance_size()
1435 sage: n = ZZ.random_element(n_max)
1436 sage: F = QuadraticField(-1, 'I')
1437 sage: X = random_matrix(F, n)
1438 sage: Y = random_matrix(F, n)
1439 sage: Xe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X)
1440 sage: Ye = ComplexMatrixEuclideanJordanAlgebra.real_embed(Y)
1441 sage: XYe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1448 raise ValueError("the matrix 'M' must be square")
1450 # We don't need any adjoined elements...
1451 field
= M
.base_ring().base_ring()
1455 a
= z
.list()[0] # real part, I guess
1456 b
= z
.list()[1] # imag part, I guess
1457 blocks
.append(matrix(field
, 2, [[a
,b
],[-b
,a
]]))
1459 return matrix
.block(field
, n
, blocks
)
1463 def real_unembed(M
):
1465 The inverse of _embed_complex_matrix().
1469 sage: from mjo.eja.eja_algebra import \
1470 ....: ComplexMatrixEuclideanJordanAlgebra
1474 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1475 ....: [-2, 1, -4, 3],
1476 ....: [ 9, 10, 11, 12],
1477 ....: [-10, 9, -12, 11] ])
1478 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(A)
1480 [ 10*I + 9 12*I + 11]
1484 Unembedding is the inverse of embedding::
1486 sage: set_random_seed()
1487 sage: F = QuadraticField(-1, 'I')
1488 sage: M = random_matrix(F, 3)
1489 sage: Me = ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1490 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1496 raise ValueError("the matrix 'M' must be square")
1497 if not n
.mod(2).is_zero():
1498 raise ValueError("the matrix 'M' must be a complex embedding")
1500 # If "M" was normalized, its base ring might have roots
1501 # adjoined and they can stick around after unembedding.
1502 field
= M
.base_ring()
1503 R
= PolynomialRing(field
, 'z')
1506 # Sage doesn't know how to embed AA into QQbar, i.e. how
1507 # to adjoin sqrt(-1) to AA.
1510 F
= field
.extension(z
**2 + 1, 'I', embedding
=CLF(-1).sqrt())
1513 # Go top-left to bottom-right (reading order), converting every
1514 # 2-by-2 block we see to a single complex element.
1516 for k
in range(n
/2):
1517 for j
in range(n
/2):
1518 submat
= M
[2*k
:2*k
+2,2*j
:2*j
+2]
1519 if submat
[0,0] != submat
[1,1]:
1520 raise ValueError('bad on-diagonal submatrix')
1521 if submat
[0,1] != -submat
[1,0]:
1522 raise ValueError('bad off-diagonal submatrix')
1523 z
= submat
[0,0] + submat
[0,1]*i
1526 return matrix(F
, n
/2, elements
)
1530 def natural_inner_product(cls
,X
,Y
):
1532 Compute a natural inner product in this algebra directly from
1537 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1541 This gives the same answer as the slow, default method implemented
1542 in :class:`MatrixEuclideanJordanAlgebra`::
1544 sage: set_random_seed()
1545 sage: J = ComplexHermitianEJA.random_instance()
1546 sage: x,y = J.random_elements(2)
1547 sage: Xe = x.natural_representation()
1548 sage: Ye = y.natural_representation()
1549 sage: X = ComplexHermitianEJA.real_unembed(Xe)
1550 sage: Y = ComplexHermitianEJA.real_unembed(Ye)
1551 sage: expected = (X*Y).trace().real()
1552 sage: actual = ComplexHermitianEJA.natural_inner_product(Xe,Ye)
1553 sage: actual == expected
1557 return RealMatrixEuclideanJordanAlgebra
.natural_inner_product(X
,Y
)/2
1560 class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra
):
1562 The rank-n simple EJA consisting of complex Hermitian n-by-n
1563 matrices over the real numbers, the usual symmetric Jordan product,
1564 and the real-part-of-trace inner product. It has dimension `n^2` over
1569 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1573 In theory, our "field" can be any subfield of the reals::
1575 sage: ComplexHermitianEJA(2, RDF)
1576 Euclidean Jordan algebra of dimension 4 over Real Double Field
1577 sage: ComplexHermitianEJA(2, RR)
1578 Euclidean Jordan algebra of dimension 4 over Real Field with
1579 53 bits of precision
1583 The dimension of this algebra is `n^2`::
1585 sage: set_random_seed()
1586 sage: n_max = ComplexHermitianEJA._max_random_instance_size()
1587 sage: n = ZZ.random_element(1, n_max)
1588 sage: J = ComplexHermitianEJA(n)
1589 sage: J.dimension() == n^2
1592 The Jordan multiplication is what we think it is::
1594 sage: set_random_seed()
1595 sage: J = ComplexHermitianEJA.random_instance()
1596 sage: x,y = J.random_elements(2)
1597 sage: actual = (x*y).natural_representation()
1598 sage: X = x.natural_representation()
1599 sage: Y = y.natural_representation()
1600 sage: expected = (X*Y + Y*X)/2
1601 sage: actual == expected
1603 sage: J(expected) == x*y
1606 We can change the generator prefix::
1608 sage: ComplexHermitianEJA(2, prefix='z').gens()
1611 Our natural basis is normalized with respect to the natural inner
1612 product unless we specify otherwise::
1614 sage: set_random_seed()
1615 sage: J = ComplexHermitianEJA.random_instance()
1616 sage: all( b.norm() == 1 for b in J.gens() )
1619 Since our natural basis is normalized with respect to the natural
1620 inner product, and since we know that this algebra is an EJA, any
1621 left-multiplication operator's matrix will be symmetric because
1622 natural->EJA basis representation is an isometry and within the EJA
1623 the operator is self-adjoint by the Jordan axiom::
1625 sage: set_random_seed()
1626 sage: x = ComplexHermitianEJA.random_instance().random_element()
1627 sage: x.operator().matrix().is_symmetric()
1630 We can construct the (trivial) algebra of rank zero::
1632 sage: ComplexHermitianEJA(0)
1633 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1638 def _denormalized_basis(cls
, n
, field
):
1640 Returns a basis for the space of complex Hermitian n-by-n matrices.
1642 Why do we embed these? Basically, because all of numerical linear
1643 algebra assumes that you're working with vectors consisting of `n`
1644 entries from a field and scalars from the same field. There's no way
1645 to tell SageMath that (for example) the vectors contain complex
1646 numbers, while the scalar field is real.
1650 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1654 sage: set_random_seed()
1655 sage: n = ZZ.random_element(1,5)
1656 sage: field = QuadraticField(2, 'sqrt2')
1657 sage: B = ComplexHermitianEJA._denormalized_basis(n, field)
1658 sage: all( M.is_symmetric() for M in B)
1662 R
= PolynomialRing(field
, 'z')
1664 F
= field
.extension(z
**2 + 1, 'I')
1667 # This is like the symmetric case, but we need to be careful:
1669 # * We want conjugate-symmetry, not just symmetry.
1670 # * The diagonal will (as a result) be real.
1674 for j
in range(i
+1):
1675 Eij
= matrix(F
, n
, lambda k
,l
: k
==i
and l
==j
)
1677 Sij
= cls
.real_embed(Eij
)
1680 # The second one has a minus because it's conjugated.
1681 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
1683 Sij_imag
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
1686 # Since we embedded these, we can drop back to the "field" that we
1687 # started with instead of the complex extension "F".
1688 return ( s
.change_ring(field
) for s
in S
)
1691 def __init__(self
, n
, field
=AA
, **kwargs
):
1692 basis
= self
._denormalized
_basis
(n
,field
)
1693 super(ComplexHermitianEJA
,self
).__init
__(field
,
1697 self
.rank
.set_cache(n
)
1700 class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1704 Embed the n-by-n quaternion matrix ``M`` into the space of real
1705 matrices of size 4n-by-4n by first sending each quaternion entry `z
1706 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
1707 c+di],[-c + di, a-bi]]`, and then embedding those into a real
1712 sage: from mjo.eja.eja_algebra import \
1713 ....: QuaternionMatrixEuclideanJordanAlgebra
1717 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1718 sage: i,j,k = Q.gens()
1719 sage: x = 1 + 2*i + 3*j + 4*k
1720 sage: M = matrix(Q, 1, [[x]])
1721 sage: QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1727 Embedding is a homomorphism (isomorphism, in fact)::
1729 sage: set_random_seed()
1730 sage: n_max = QuaternionMatrixEuclideanJordanAlgebra._max_random_instance_size()
1731 sage: n = ZZ.random_element(n_max)
1732 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1733 sage: X = random_matrix(Q, n)
1734 sage: Y = random_matrix(Q, n)
1735 sage: Xe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X)
1736 sage: Ye = QuaternionMatrixEuclideanJordanAlgebra.real_embed(Y)
1737 sage: XYe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1742 quaternions
= M
.base_ring()
1745 raise ValueError("the matrix 'M' must be square")
1747 F
= QuadraticField(-1, 'I')
1752 t
= z
.coefficient_tuple()
1757 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
1758 [-c
+ d
*i
, a
- b
*i
]])
1759 realM
= ComplexMatrixEuclideanJordanAlgebra
.real_embed(cplxM
)
1760 blocks
.append(realM
)
1762 # We should have real entries by now, so use the realest field
1763 # we've got for the return value.
1764 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
1769 def real_unembed(M
):
1771 The inverse of _embed_quaternion_matrix().
1775 sage: from mjo.eja.eja_algebra import \
1776 ....: QuaternionMatrixEuclideanJordanAlgebra
1780 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
1781 ....: [-2, 1, -4, 3],
1782 ....: [-3, 4, 1, -2],
1783 ....: [-4, -3, 2, 1]])
1784 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(M)
1785 [1 + 2*i + 3*j + 4*k]
1789 Unembedding is the inverse of embedding::
1791 sage: set_random_seed()
1792 sage: Q = QuaternionAlgebra(QQ, -1, -1)
1793 sage: M = random_matrix(Q, 3)
1794 sage: Me = QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1795 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1801 raise ValueError("the matrix 'M' must be square")
1802 if not n
.mod(4).is_zero():
1803 raise ValueError("the matrix 'M' must be a quaternion embedding")
1805 # Use the base ring of the matrix to ensure that its entries can be
1806 # multiplied by elements of the quaternion algebra.
1807 field
= M
.base_ring()
1808 Q
= QuaternionAlgebra(field
,-1,-1)
1811 # Go top-left to bottom-right (reading order), converting every
1812 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
1815 for l
in range(n
/4):
1816 for m
in range(n
/4):
1817 submat
= ComplexMatrixEuclideanJordanAlgebra
.real_unembed(
1818 M
[4*l
:4*l
+4,4*m
:4*m
+4] )
1819 if submat
[0,0] != submat
[1,1].conjugate():
1820 raise ValueError('bad on-diagonal submatrix')
1821 if submat
[0,1] != -submat
[1,0].conjugate():
1822 raise ValueError('bad off-diagonal submatrix')
1823 z
= submat
[0,0].real()
1824 z
+= submat
[0,0].imag()*i
1825 z
+= submat
[0,1].real()*j
1826 z
+= submat
[0,1].imag()*k
1829 return matrix(Q
, n
/4, elements
)
1833 def natural_inner_product(cls
,X
,Y
):
1835 Compute a natural inner product in this algebra directly from
1840 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1844 This gives the same answer as the slow, default method implemented
1845 in :class:`MatrixEuclideanJordanAlgebra`::
1847 sage: set_random_seed()
1848 sage: J = QuaternionHermitianEJA.random_instance()
1849 sage: x,y = J.random_elements(2)
1850 sage: Xe = x.natural_representation()
1851 sage: Ye = y.natural_representation()
1852 sage: X = QuaternionHermitianEJA.real_unembed(Xe)
1853 sage: Y = QuaternionHermitianEJA.real_unembed(Ye)
1854 sage: expected = (X*Y).trace().coefficient_tuple()[0]
1855 sage: actual = QuaternionHermitianEJA.natural_inner_product(Xe,Ye)
1856 sage: actual == expected
1860 return RealMatrixEuclideanJordanAlgebra
.natural_inner_product(X
,Y
)/4
1863 class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra
):
1865 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
1866 matrices, the usual symmetric Jordan product, and the
1867 real-part-of-trace inner product. It has dimension `2n^2 - n` over
1872 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1876 In theory, our "field" can be any subfield of the reals::
1878 sage: QuaternionHermitianEJA(2, RDF)
1879 Euclidean Jordan algebra of dimension 6 over Real Double Field
1880 sage: QuaternionHermitianEJA(2, RR)
1881 Euclidean Jordan algebra of dimension 6 over Real Field with
1882 53 bits of precision
1886 The dimension of this algebra is `2*n^2 - n`::
1888 sage: set_random_seed()
1889 sage: n_max = QuaternionHermitianEJA._max_random_instance_size()
1890 sage: n = ZZ.random_element(1, n_max)
1891 sage: J = QuaternionHermitianEJA(n)
1892 sage: J.dimension() == 2*(n^2) - n
1895 The Jordan multiplication is what we think it is::
1897 sage: set_random_seed()
1898 sage: J = QuaternionHermitianEJA.random_instance()
1899 sage: x,y = J.random_elements(2)
1900 sage: actual = (x*y).natural_representation()
1901 sage: X = x.natural_representation()
1902 sage: Y = y.natural_representation()
1903 sage: expected = (X*Y + Y*X)/2
1904 sage: actual == expected
1906 sage: J(expected) == x*y
1909 We can change the generator prefix::
1911 sage: QuaternionHermitianEJA(2, prefix='a').gens()
1912 (a0, a1, a2, a3, a4, a5)
1914 Our natural basis is normalized with respect to the natural inner
1915 product unless we specify otherwise::
1917 sage: set_random_seed()
1918 sage: J = QuaternionHermitianEJA.random_instance()
1919 sage: all( b.norm() == 1 for b in J.gens() )
1922 Since our natural basis is normalized with respect to the natural
1923 inner product, and since we know that this algebra is an EJA, any
1924 left-multiplication operator's matrix will be symmetric because
1925 natural->EJA basis representation is an isometry and within the EJA
1926 the operator is self-adjoint by the Jordan axiom::
1928 sage: set_random_seed()
1929 sage: x = QuaternionHermitianEJA.random_instance().random_element()
1930 sage: x.operator().matrix().is_symmetric()
1933 We can construct the (trivial) algebra of rank zero::
1935 sage: QuaternionHermitianEJA(0)
1936 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1940 def _denormalized_basis(cls
, n
, field
):
1942 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
1944 Why do we embed these? Basically, because all of numerical
1945 linear algebra assumes that you're working with vectors consisting
1946 of `n` entries from a field and scalars from the same field. There's
1947 no way to tell SageMath that (for example) the vectors contain
1948 complex numbers, while the scalar field is real.
1952 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1956 sage: set_random_seed()
1957 sage: n = ZZ.random_element(1,5)
1958 sage: B = QuaternionHermitianEJA._denormalized_basis(n,QQ)
1959 sage: all( M.is_symmetric() for M in B )
1963 Q
= QuaternionAlgebra(QQ
,-1,-1)
1966 # This is like the symmetric case, but we need to be careful:
1968 # * We want conjugate-symmetry, not just symmetry.
1969 # * The diagonal will (as a result) be real.
1973 for j
in range(i
+1):
1974 Eij
= matrix(Q
, n
, lambda k
,l
: k
==i
and l
==j
)
1976 Sij
= cls
.real_embed(Eij
)
1979 # The second, third, and fourth ones have a minus
1980 # because they're conjugated.
1981 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
1983 Sij_I
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
1985 Sij_J
= cls
.real_embed(J
*Eij
- J
*Eij
.transpose())
1987 Sij_K
= cls
.real_embed(K
*Eij
- K
*Eij
.transpose())
1990 # Since we embedded these, we can drop back to the "field" that we
1991 # started with instead of the quaternion algebra "Q".
1992 return ( s
.change_ring(field
) for s
in S
)
1995 def __init__(self
, n
, field
=AA
, **kwargs
):
1996 basis
= self
._denormalized
_basis
(n
,field
)
1997 super(QuaternionHermitianEJA
,self
).__init
__(field
,
2001 self
.rank
.set_cache(n
)
2004 class HadamardEJA(RationalBasisEuclideanJordanAlgebra
):
2006 Return the Euclidean Jordan Algebra corresponding to the set
2007 `R^n` under the Hadamard product.
2009 Note: this is nothing more than the Cartesian product of ``n``
2010 copies of the spin algebra. Once Cartesian product algebras
2011 are implemented, this can go.
2015 sage: from mjo.eja.eja_algebra import HadamardEJA
2019 This multiplication table can be verified by hand::
2021 sage: J = HadamardEJA(3)
2022 sage: e0,e1,e2 = J.gens()
2038 We can change the generator prefix::
2040 sage: HadamardEJA(3, prefix='r').gens()
2044 def __init__(self
, n
, field
=AA
, **kwargs
):
2045 V
= VectorSpace(field
, n
)
2046 mult_table
= [ [ V
.gen(i
)*(i
== j
) for j
in range(n
) ]
2049 super(HadamardEJA
, self
).__init
__(field
,
2053 self
.rank
.set_cache(n
)
2056 def _max_random_instance_size():
2059 def inner_product(self
, x
, y
):
2061 Faster to reimplement than to use natural representations.
2065 sage: from mjo.eja.eja_algebra import HadamardEJA
2069 Ensure that this is the usual inner product for the algebras
2072 sage: set_random_seed()
2073 sage: J = HadamardEJA.random_instance()
2074 sage: x,y = J.random_elements(2)
2075 sage: X = x.natural_representation()
2076 sage: Y = y.natural_representation()
2077 sage: x.inner_product(y) == J.natural_inner_product(X,Y)
2081 return x
.to_vector().inner_product(y
.to_vector())
2084 class BilinearFormEJA(RationalBasisEuclideanJordanAlgebra
):
2086 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2087 with the half-trace inner product and jordan product ``x*y =
2088 (<Bx,y>,y_bar>, x0*y_bar + y0*x_bar)`` where `B = 1 \times B22` is
2089 a symmetric positive-definite "bilinear form" matrix. Its
2090 dimension is the size of `B`, and it has rank two in dimensions
2091 larger than two. It reduces to the ``JordanSpinEJA`` when `B` is
2092 the identity matrix of order ``n``.
2094 We insist that the one-by-one upper-left identity block of `B` be
2095 passed in as well so that we can be passed a matrix of size zero
2096 to construct a trivial algebra.
2100 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
2101 ....: JordanSpinEJA)
2105 When no bilinear form is specified, the identity matrix is used,
2106 and the resulting algebra is the Jordan spin algebra::
2108 sage: B = matrix.identity(AA,3)
2109 sage: J0 = BilinearFormEJA(B)
2110 sage: J1 = JordanSpinEJA(3)
2111 sage: J0.multiplication_table() == J0.multiplication_table()
2116 We can create a zero-dimensional algebra::
2118 sage: B = matrix.identity(AA,0)
2119 sage: J = BilinearFormEJA(B)
2123 We can check the multiplication condition given in the Jordan, von
2124 Neumann, and Wigner paper (and also discussed on my "On the
2125 symmetry..." paper). Note that this relies heavily on the standard
2126 choice of basis, as does anything utilizing the bilinear form matrix::
2128 sage: set_random_seed()
2129 sage: n = ZZ.random_element(5)
2130 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2131 sage: B11 = matrix.identity(QQ,1)
2132 sage: B22 = M.transpose()*M
2133 sage: B = block_matrix(2,2,[ [B11,0 ],
2135 sage: J = BilinearFormEJA(B)
2136 sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
2137 sage: V = J.vector_space()
2138 sage: sis = [ J.from_vector(V([0] + (M.inverse()*ei).list()))
2139 ....: for ei in eis ]
2140 sage: actual = [ sis[i]*sis[j]
2141 ....: for i in range(n-1)
2142 ....: for j in range(n-1) ]
2143 sage: expected = [ J.one() if i == j else J.zero()
2144 ....: for i in range(n-1)
2145 ....: for j in range(n-1) ]
2146 sage: actual == expected
2149 def __init__(self
, B
, field
=AA
, **kwargs
):
2153 if not B
.is_positive_definite():
2154 raise TypeError("matrix B is not positive-definite")
2156 V
= VectorSpace(field
, n
)
2157 mult_table
= [[V
.zero() for j
in range(n
)] for i
in range(n
)]
2166 z0
= (B
*x
).inner_product(y
)
2167 zbar
= y0
*xbar
+ x0
*ybar
2168 z
= V([z0
] + zbar
.list())
2169 mult_table
[i
][j
] = z
2171 # The rank of this algebra is two, unless we're in a
2172 # one-dimensional ambient space (because the rank is bounded
2173 # by the ambient dimension).
2174 super(BilinearFormEJA
, self
).__init
__(field
,
2178 self
.rank
.set_cache(min(n
,2))
2181 def _max_random_instance_size():
2185 def random_instance(cls
, field
=AA
, **kwargs
):
2187 Return a random instance of this algebra.
2189 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2191 # Special case needed since we use (n-1) below.
2192 B
= matrix
.identity(field
, 0)
2193 return cls(B
, field
, **kwargs
)
2195 B11
= matrix
.identity(field
,1)
2196 M
= matrix
.random(field
, n
-1)
2197 I
= matrix
.identity(field
, n
-1)
2198 alpha
= field
.zero()
2199 while alpha
.is_zero():
2200 alpha
= field
.random_element().abs()
2201 B22
= M
.transpose()*M
+ alpha
*I
2203 from sage
.matrix
.special
import block_matrix
2204 B
= block_matrix(2,2, [ [B11
, ZZ(0) ],
2207 return cls(B
, field
, **kwargs
)
2209 def inner_product(self
, x
, y
):
2211 Half of the trace inner product.
2213 This is defined so that the special case of the Jordan spin
2214 algebra gets the usual inner product.
2218 sage: from mjo.eja.eja_algebra import BilinearFormEJA
2222 Ensure that this is one-half of the trace inner-product when
2223 the algebra isn't just the reals (when ``n`` isn't one). This
2224 is in Faraut and Koranyi, and also my "On the symmetry..."
2227 sage: set_random_seed()
2228 sage: J = BilinearFormEJA.random_instance()
2229 sage: n = J.dimension()
2230 sage: x = J.random_element()
2231 sage: y = J.random_element()
2232 sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
2236 return (self
._B
*x
.to_vector()).inner_product(y
.to_vector())
2239 class JordanSpinEJA(BilinearFormEJA
):
2241 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2242 with the usual inner product and jordan product ``x*y =
2243 (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2248 sage: from mjo.eja.eja_algebra import JordanSpinEJA
2252 This multiplication table can be verified by hand::
2254 sage: J = JordanSpinEJA(4)
2255 sage: e0,e1,e2,e3 = J.gens()
2271 We can change the generator prefix::
2273 sage: JordanSpinEJA(2, prefix='B').gens()
2278 Ensure that we have the usual inner product on `R^n`::
2280 sage: set_random_seed()
2281 sage: J = JordanSpinEJA.random_instance()
2282 sage: x,y = J.random_elements(2)
2283 sage: X = x.natural_representation()
2284 sage: Y = y.natural_representation()
2285 sage: x.inner_product(y) == J.natural_inner_product(X,Y)
2289 def __init__(self
, n
, field
=AA
, **kwargs
):
2290 # This is a special case of the BilinearFormEJA with the identity
2291 # matrix as its bilinear form.
2292 B
= matrix
.identity(field
, n
)
2293 super(JordanSpinEJA
, self
).__init
__(B
, field
, **kwargs
)
2296 def random_instance(cls
, field
=AA
, **kwargs
):
2298 Return a random instance of this type of algebra.
2300 Needed here to override the implementation for ``BilinearFormEJA``.
2302 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2303 return cls(n
, field
, **kwargs
)
2306 class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra
):
2308 The trivial Euclidean Jordan algebra consisting of only a zero element.
2312 sage: from mjo.eja.eja_algebra import TrivialEJA
2316 sage: J = TrivialEJA()
2323 sage: 7*J.one()*12*J.one()
2325 sage: J.one().inner_product(J.one())
2327 sage: J.one().norm()
2329 sage: J.one().subalgebra_generated_by()
2330 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2335 def __init__(self
, field
=AA
, **kwargs
):
2337 super(TrivialEJA
, self
).__init
__(field
,
2341 # The rank is zero using my definition, namely the dimension of the
2342 # largest subalgebra generated by any element.
2343 self
.rank
.set_cache(0)
2346 def random_instance(cls
, field
=AA
, **kwargs
):
2347 # We don't take a "size" argument so the superclass method is
2348 # inappropriate for us.
2349 return cls(field
, **kwargs
)
2351 class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra
):
2353 The external (orthogonal) direct sum of two other Euclidean Jordan
2354 algebras. Essentially the Cartesian product of its two factors.
2355 Every Euclidean Jordan algebra decomposes into an orthogonal
2356 direct sum of simple Euclidean Jordan algebras, so no generality
2357 is lost by providing only this construction.
2361 sage: from mjo.eja.eja_algebra import (HadamardEJA,
2362 ....: RealSymmetricEJA,
2367 sage: J1 = HadamardEJA(2)
2368 sage: J2 = RealSymmetricEJA(3)
2369 sage: J = DirectSumEJA(J1,J2)
2376 def __init__(self
, J1
, J2
, field
=AA
, **kwargs
):
2377 self
._factors
= (J1
, J2
)
2381 V
= VectorSpace(field
, n
)
2382 mult_table
= [ [ V
.zero() for j
in range(n
) ]
2386 p
= (J1
.monomial(i
)*J1
.monomial(j
)).to_vector()
2387 mult_table
[i
][j
] = V(p
.list() + [field
.zero()]*n2
)
2391 p
= (J2
.monomial(i
)*J2
.monomial(j
)).to_vector()
2392 mult_table
[n1
+i
][n1
+j
] = V([field
.zero()]*n1
+ p
.list())
2394 super(DirectSumEJA
, self
).__init
__(field
,
2398 self
.rank
.set_cache(J1
.rank() + J2
.rank())
2403 Return the pair of this algebra's factors.
2407 sage: from mjo.eja.eja_algebra import (HadamardEJA,
2408 ....: JordanSpinEJA,
2413 sage: J1 = HadamardEJA(2,QQ)
2414 sage: J2 = JordanSpinEJA(3,QQ)
2415 sage: J = DirectSumEJA(J1,J2)
2417 (Euclidean Jordan algebra of dimension 2 over Rational Field,
2418 Euclidean Jordan algebra of dimension 3 over Rational Field)
2421 return self
._factors
2423 def projections(self
):
2425 Return a pair of projections onto this algebra's factors.
2429 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
2430 ....: ComplexHermitianEJA,
2435 sage: J1 = JordanSpinEJA(2)
2436 sage: J2 = ComplexHermitianEJA(2)
2437 sage: J = DirectSumEJA(J1,J2)
2438 sage: (pi_left, pi_right) = J.projections()
2439 sage: J.one().to_vector()
2441 sage: pi_left(J.one()).to_vector()
2443 sage: pi_right(J.one()).to_vector()
2447 (J1
,J2
) = self
.factors()
2449 pi_left
= lambda x
: J1
.from_vector(x
.to_vector()[:n
])
2450 pi_right
= lambda x
: J2
.from_vector(x
.to_vector()[n
:])
2451 return (pi_left
, pi_right
)
2453 def inclusions(self
):
2455 Return the pair of inclusion maps from our factors into us.
2459 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
2460 ....: RealSymmetricEJA,
2465 sage: J1 = JordanSpinEJA(3)
2466 sage: J2 = RealSymmetricEJA(2)
2467 sage: J = DirectSumEJA(J1,J2)
2468 sage: (iota_left, iota_right) = J.inclusions()
2469 sage: iota_left(J1.zero()) == J.zero()
2471 sage: iota_right(J2.zero()) == J.zero()
2473 sage: J1.one().to_vector()
2475 sage: iota_left(J1.one()).to_vector()
2477 sage: J2.one().to_vector()
2479 sage: iota_right(J2.one()).to_vector()
2481 sage: J.one().to_vector()
2485 (J1
,J2
) = self
.factors()
2487 V_basis
= self
.vector_space().basis()
2488 I1
= matrix
.column(self
.base_ring(), V_basis
[:n
])
2489 I2
= matrix
.column(self
.base_ring(), V_basis
[n
:])
2490 iota_left
= lambda x
: self
.from_vector(I1
*x
.to_vector())
2491 iota_right
= lambda x
: self
.from_vector(I2
*+x
.to_vector())
2492 return (iota_left
, iota_right
)
2494 def inner_product(self
, x
, y
):
2496 The standard Cartesian inner-product.
2498 We project ``x`` and ``y`` onto our factors, and add up the
2499 inner-products from the subalgebras.
2504 sage: from mjo.eja.eja_algebra import (HadamardEJA,
2505 ....: QuaternionHermitianEJA,
2510 sage: J1 = HadamardEJA(3)
2511 sage: J2 = QuaternionHermitianEJA(2,QQ,normalize_basis=False)
2512 sage: J = DirectSumEJA(J1,J2)
2517 sage: x1.inner_product(x2)
2519 sage: y1.inner_product(y2)
2521 sage: J.one().inner_product(J.one())
2525 (pi_left
, pi_right
) = self
.projections()
2531 return (x1
.inner_product(y1
) + x2
.inner_product(y2
))