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 basis 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
.column(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 should work.
613 coeffs
= A
.solve_right(b
)
614 return self
.linear_combination(zip(self
.gens(), coeffs
))
617 def peirce_decomposition(self
, c
):
619 The Peirce decomposition of this algebra relative to the
622 In the future, this can be extended to a complete system of
623 orthogonal idempotents.
627 - ``c`` -- an idempotent of this algebra.
631 A triple (J0, J5, J1) containing two subalgebras and one subspace
634 - ``J0`` -- the algebra on the eigenspace of ``c.operator()``
635 corresponding to the eigenvalue zero.
637 - ``J5`` -- the eigenspace (NOT a subalgebra) of ``c.operator()``
638 corresponding to the eigenvalue one-half.
640 - ``J1`` -- the algebra on the eigenspace of ``c.operator()``
641 corresponding to the eigenvalue one.
643 These are the only possible eigenspaces for that operator, and this
644 algebra is a direct sum of them. The spaces ``J0`` and ``J1`` are
645 orthogonal, and are subalgebras of this algebra with the appropriate
650 sage: from mjo.eja.eja_algebra import random_eja, RealSymmetricEJA
654 The canonical example comes from the symmetric matrices, which
655 decompose into diagonal and off-diagonal parts::
657 sage: J = RealSymmetricEJA(3)
658 sage: C = matrix(QQ, [ [1,0,0],
662 sage: J0,J5,J1 = J.peirce_decomposition(c)
664 Euclidean Jordan algebra of dimension 1...
666 Vector space of degree 6 and dimension 2...
668 Euclidean Jordan algebra of dimension 3...
669 sage: J0.one().natural_representation()
673 sage: orig_df = AA.options.display_format
674 sage: AA.options.display_format = 'radical'
675 sage: J.from_vector(J5.basis()[0]).natural_representation()
679 sage: J.from_vector(J5.basis()[1]).natural_representation()
683 sage: AA.options.display_format = orig_df
684 sage: J1.one().natural_representation()
691 Every algebra decomposes trivially with respect to its identity
694 sage: set_random_seed()
695 sage: J = random_eja()
696 sage: J0,J5,J1 = J.peirce_decomposition(J.one())
697 sage: J0.dimension() == 0 and J5.dimension() == 0
699 sage: J1.superalgebra() == J and J1.dimension() == J.dimension()
702 The decomposition is into eigenspaces, and its components are
703 therefore necessarily orthogonal. Moreover, the identity
704 elements in the two subalgebras are the projections onto their
705 respective subspaces of the superalgebra's identity element::
707 sage: set_random_seed()
708 sage: J = random_eja()
709 sage: x = J.random_element()
710 sage: if not J.is_trivial():
711 ....: while x.is_nilpotent():
712 ....: x = J.random_element()
713 sage: c = x.subalgebra_idempotent()
714 sage: J0,J5,J1 = J.peirce_decomposition(c)
716 sage: for (w,y,z) in zip(J0.basis(), J5.basis(), J1.basis()):
717 ....: w = w.superalgebra_element()
718 ....: y = J.from_vector(y)
719 ....: z = z.superalgebra_element()
720 ....: ipsum += w.inner_product(y).abs()
721 ....: ipsum += w.inner_product(z).abs()
722 ....: ipsum += y.inner_product(z).abs()
725 sage: J1(c) == J1.one()
727 sage: J0(J.one() - c) == J0.one()
731 if not c
.is_idempotent():
732 raise ValueError("element is not idempotent: %s" % c
)
734 # Default these to what they should be if they turn out to be
735 # trivial, because eigenspaces_left() won't return eigenvalues
736 # corresponding to trivial spaces (e.g. it returns only the
737 # eigenspace corresponding to lambda=1 if you take the
738 # decomposition relative to the identity element).
739 trivial
= FiniteDimensionalEuclideanJordanSubalgebra(self
, ())
740 J0
= trivial
# eigenvalue zero
741 J5
= VectorSpace(self
.base_ring(), 0) # eigenvalue one-half
742 J1
= trivial
# eigenvalue one
744 for (eigval
, eigspace
) in c
.operator().matrix().right_eigenspaces():
745 if eigval
== ~
(self
.base_ring()(2)):
748 gens
= tuple( self
.from_vector(b
) for b
in eigspace
.basis() )
749 subalg
= FiniteDimensionalEuclideanJordanSubalgebra(self
,
757 raise ValueError("unexpected eigenvalue: %s" % eigval
)
762 def random_element(self
, thorough
=False):
764 Return a random element of this algebra.
766 Our algebra superclass method only returns a linear
767 combination of at most two basis elements. We instead
768 want the vector space "random element" method that
769 returns a more diverse selection.
773 - ``thorough`` -- (boolean; default False) whether or not we
774 should generate irrational coefficients for the random
775 element when our base ring is irrational; this slows the
776 algebra operations to a crawl, but any truly random method
780 # For a general base ring... maybe we can trust this to do the
781 # right thing? Unlikely, but.
782 V
= self
.vector_space()
783 v
= V
.random_element()
785 if self
.base_ring() is AA
:
786 # The "random element" method of the algebraic reals is
787 # stupid at the moment, and only returns integers between
788 # -2 and 2, inclusive:
790 # https://trac.sagemath.org/ticket/30875
792 # Instead, we implement our own "random vector" method,
793 # and then coerce that into the algebra. We use the vector
794 # space degree here instead of the dimension because a
795 # subalgebra could (for example) be spanned by only two
796 # vectors, each with five coordinates. We need to
797 # generate all five coordinates.
799 v
*= QQbar
.random_element().real()
801 v
*= QQ
.random_element()
803 return self
.from_vector(V
.coordinate_vector(v
))
805 def random_elements(self
, count
, thorough
=False):
807 Return ``count`` random elements as a tuple.
811 - ``thorough`` -- (boolean; default False) whether or not we
812 should generate irrational coefficients for the random
813 elements when our base ring is irrational; this slows the
814 algebra operations to a crawl, but any truly random method
819 sage: from mjo.eja.eja_algebra import JordanSpinEJA
823 sage: J = JordanSpinEJA(3)
824 sage: x,y,z = J.random_elements(3)
825 sage: all( [ x in J, y in J, z in J ])
827 sage: len( J.random_elements(10) ) == 10
831 return tuple( self
.random_element(thorough
)
832 for idx
in range(count
) )
835 def random_instance(cls
, field
=AA
, **kwargs
):
837 Return a random instance of this type of algebra.
839 Beware, this will crash for "most instances" because the
840 constructor below looks wrong.
842 if cls
is TrivialEJA
:
843 # The TrivialEJA class doesn't take an "n" argument because
847 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
848 return cls(n
, field
, **kwargs
)
851 def _charpoly_coefficients(self
):
853 The `r` polynomial coefficients of the "characteristic polynomial
857 var_names
= [ "X" + str(z
) for z
in range(1,n
+1) ]
858 R
= PolynomialRing(self
.base_ring(), var_names
)
860 F
= R
.fraction_field()
863 # From a result in my book, these are the entries of the
864 # basis representation of L_x.
865 return sum( vars[k
]*self
.monomial(k
).operator().matrix()[i
,j
]
868 L_x
= matrix(F
, n
, n
, L_x_i_j
)
871 if self
.rank
.is_in_cache():
873 # There's no need to pad the system with redundant
874 # columns if we *know* they'll be redundant.
877 # Compute an extra power in case the rank is equal to
878 # the dimension (otherwise, we would stop at x^(r-1)).
879 x_powers
= [ (L_x
**k
)*self
.one().to_vector()
880 for k
in range(n
+1) ]
881 A
= matrix
.column(F
, x_powers
[:n
])
882 AE
= A
.extended_echelon_form()
889 # The theory says that only the first "r" coefficients are
890 # nonzero, and they actually live in the original polynomial
891 # ring and not the fraction field. We negate them because
892 # in the actual characteristic polynomial, they get moved
893 # to the other side where x^r lives.
894 return -A_rref
.solve_right(E
*b
).change_ring(R
)[:r
]
899 Return the rank of this EJA.
901 This is a cached method because we know the rank a priori for
902 all of the algebras we can construct. Thus we can avoid the
903 expensive ``_charpoly_coefficients()`` call unless we truly
904 need to compute the whole characteristic polynomial.
908 sage: from mjo.eja.eja_algebra import (HadamardEJA,
910 ....: RealSymmetricEJA,
911 ....: ComplexHermitianEJA,
912 ....: QuaternionHermitianEJA,
917 The rank of the Jordan spin algebra is always two::
919 sage: JordanSpinEJA(2).rank()
921 sage: JordanSpinEJA(3).rank()
923 sage: JordanSpinEJA(4).rank()
926 The rank of the `n`-by-`n` Hermitian real, complex, or
927 quaternion matrices is `n`::
929 sage: RealSymmetricEJA(4).rank()
931 sage: ComplexHermitianEJA(3).rank()
933 sage: QuaternionHermitianEJA(2).rank()
938 Ensure that every EJA that we know how to construct has a
939 positive integer rank, unless the algebra is trivial in
940 which case its rank will be zero::
942 sage: set_random_seed()
943 sage: J = random_eja()
947 sage: r > 0 or (r == 0 and J.is_trivial())
950 Ensure that computing the rank actually works, since the ranks
951 of all simple algebras are known and will be cached by default::
953 sage: J = HadamardEJA(4)
954 sage: J.rank.clear_cache()
960 sage: J = JordanSpinEJA(4)
961 sage: J.rank.clear_cache()
967 sage: J = RealSymmetricEJA(3)
968 sage: J.rank.clear_cache()
974 sage: J = ComplexHermitianEJA(2)
975 sage: J.rank.clear_cache()
981 sage: J = QuaternionHermitianEJA(2)
982 sage: J.rank.clear_cache()
986 return len(self
._charpoly
_coefficients
())
989 def vector_space(self
):
991 Return the vector space that underlies this algebra.
995 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
999 sage: J = RealSymmetricEJA(2)
1000 sage: J.vector_space()
1001 Vector space of dimension 3 over...
1004 return self
.zero().to_vector().parent().ambient_vector_space()
1007 Element
= FiniteDimensionalEuclideanJordanAlgebraElement
1011 def random_eja(field
=AA
):
1013 Return a "random" finite-dimensional Euclidean Jordan Algebra.
1017 sage: from mjo.eja.eja_algebra import random_eja
1022 Euclidean Jordan algebra of dimension...
1025 classname
= choice([TrivialEJA
,
1029 ComplexHermitianEJA
,
1030 QuaternionHermitianEJA
])
1031 return classname
.random_instance(field
=field
)
1036 class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra
):
1038 Algebras whose basis consists of vectors with rational
1039 entries. Equivalently, algebras whose multiplication tables
1040 contain only rational coefficients.
1042 When an EJA has a basis that can be made rational, we can speed up
1043 the computation of its characteristic polynomial by doing it over
1044 ``QQ``. All of the named EJA constructors that we provide fall
1048 def _charpoly_coefficients(self
):
1050 Override the parent method with something that tries to compute
1051 over a faster (non-extension) field.
1055 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1059 The base ring of the resulting polynomial coefficients is what
1060 it should be, and not the rationals (unless the algebra was
1061 already over the rationals)::
1063 sage: J = JordanSpinEJA(3)
1064 sage: J._charpoly_coefficients()
1065 (X1^2 - X2^2 - X3^2, -2*X1)
1066 sage: a0 = J._charpoly_coefficients()[0]
1068 Algebraic Real Field
1069 sage: a0.base_ring()
1070 Algebraic Real Field
1073 if self
.base_ring() is QQ
:
1074 # There's no need to construct *another* algebra over the
1075 # rationals if this one is already over the rationals.
1076 superclass
= super(RationalBasisEuclideanJordanAlgebra
, self
)
1077 return superclass
._charpoly
_coefficients
()
1080 map(lambda x
: x
.to_vector(), ls
)
1081 for ls
in self
._multiplication
_table
1084 # Do the computation over the rationals. The answer will be
1085 # the same, because our basis coordinates are (essentially)
1087 J
= FiniteDimensionalEuclideanJordanAlgebra(QQ
,
1091 a
= J
._charpoly
_coefficients
()
1092 return tuple(map(lambda x
: x
.change_ring(self
.base_ring()), a
))
1095 class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra
):
1097 def _max_random_instance_size():
1098 # Play it safe, since this will be squared and the underlying
1099 # field can have dimension 4 (quaternions) too.
1102 def __init__(self
, field
, basis
, normalize_basis
=True, **kwargs
):
1104 Compared to the superclass constructor, we take a basis instead of
1105 a multiplication table because the latter can be computed in terms
1106 of the former when the product is known (like it is here).
1108 # Used in this class's fast _charpoly_coefficients() override.
1109 self
._basis
_normalizers
= None
1111 # We're going to loop through this a few times, so now's a good
1112 # time to ensure that it isn't a generator expression.
1113 basis
= tuple(basis
)
1115 if len(basis
) > 1 and normalize_basis
:
1116 # We'll need sqrt(2) to normalize the basis, and this
1117 # winds up in the multiplication table, so the whole
1118 # algebra needs to be over the field extension.
1119 R
= PolynomialRing(field
, 'z')
1122 if p
.is_irreducible():
1123 field
= field
.extension(p
, 'sqrt2', embedding
=RLF(2).sqrt())
1124 basis
= tuple( s
.change_ring(field
) for s
in basis
)
1125 self
._basis
_normalizers
= tuple(
1126 ~
(self
.natural_inner_product(s
,s
).sqrt()) for s
in basis
)
1127 basis
= tuple(s
*c
for (s
,c
) in zip(basis
,self
._basis
_normalizers
))
1129 Qs
= self
.multiplication_table_from_matrix_basis(basis
)
1131 super(MatrixEuclideanJordanAlgebra
, self
).__init
__(field
,
1133 natural_basis
=basis
,
1138 def _charpoly_coefficients(self
):
1140 Override the parent method with something that tries to compute
1141 over a faster (non-extension) field.
1143 if self
._basis
_normalizers
is None or self
.base_ring() is QQ
:
1144 # We didn't normalize, or the basis we started with had
1145 # entries in a nice field already. Just compute the thing.
1146 return super(MatrixEuclideanJordanAlgebra
, self
)._charpoly
_coefficients
()
1148 basis
= ( (b
/n
) for (b
,n
) in zip(self
.natural_basis(),
1149 self
._basis
_normalizers
) )
1151 # Do this over the rationals and convert back at the end.
1152 # Only works because we know the entries of the basis are
1153 # integers. The argument ``check_axioms=False`` is required
1154 # because the trace inner-product method for this
1155 # class is a stub and can't actually be checked.
1156 J
= MatrixEuclideanJordanAlgebra(QQ
,
1158 normalize_basis
=False,
1161 a
= J
._charpoly
_coefficients
()
1163 # Unfortunately, changing the basis does change the
1164 # coefficients of the characteristic polynomial, but since
1165 # these are really the coefficients of the "characteristic
1166 # polynomial of" function, everything is still nice and
1167 # unevaluated. It's therefore "obvious" how scaling the
1168 # basis affects the coordinate variables X1, X2, et
1169 # cetera. Scaling the first basis vector up by "n" adds a
1170 # factor of 1/n into every "X1" term, for example. So here
1171 # we simply undo the basis_normalizer scaling that we
1172 # performed earlier.
1174 # The a[0] access here is safe because trivial algebras
1175 # won't have any basis normalizers and therefore won't
1176 # make it to this "else" branch.
1177 XS
= a
[0].parent().gens()
1178 subs_dict
= { XS
[i
]: self
._basis
_normalizers
[i
]*XS
[i
]
1179 for i
in range(len(XS
)) }
1180 return tuple( a_i
.subs(subs_dict
) for a_i
in a
)
1184 def multiplication_table_from_matrix_basis(basis
):
1186 At least three of the five simple Euclidean Jordan algebras have the
1187 symmetric multiplication (A,B) |-> (AB + BA)/2, where the
1188 multiplication on the right is matrix multiplication. Given a basis
1189 for the underlying matrix space, this function returns a
1190 multiplication table (obtained by looping through the basis
1191 elements) for an algebra of those matrices.
1193 # In S^2, for example, we nominally have four coordinates even
1194 # though the space is of dimension three only. The vector space V
1195 # is supposed to hold the entire long vector, and the subspace W
1196 # of V will be spanned by the vectors that arise from symmetric
1197 # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
1201 field
= basis
[0].base_ring()
1202 dimension
= basis
[0].nrows()
1204 V
= VectorSpace(field
, dimension
**2)
1205 W
= V
.span_of_basis( _mat2vec(s
) for s
in basis
)
1207 mult_table
= [[W
.zero() for j
in range(n
)] for i
in range(n
)]
1210 mat_entry
= (basis
[i
]*basis
[j
] + basis
[j
]*basis
[i
])/2
1211 mult_table
[i
][j
] = W
.coordinate_vector(_mat2vec(mat_entry
))
1219 Embed the matrix ``M`` into a space of real matrices.
1221 The matrix ``M`` can have entries in any field at the moment:
1222 the real numbers, complex numbers, or quaternions. And although
1223 they are not a field, we can probably support octonions at some
1224 point, too. This function returns a real matrix that "acts like"
1225 the original with respect to matrix multiplication; i.e.
1227 real_embed(M*N) = real_embed(M)*real_embed(N)
1230 raise NotImplementedError
1234 def real_unembed(M
):
1236 The inverse of :meth:`real_embed`.
1238 raise NotImplementedError
1242 def natural_inner_product(cls
,X
,Y
):
1243 Xu
= cls
.real_unembed(X
)
1244 Yu
= cls
.real_unembed(Y
)
1245 tr
= (Xu
*Yu
).trace()
1248 # Works in QQ, AA, RDF, et cetera.
1250 except AttributeError:
1251 # A quaternion doesn't have a real() method, but does
1252 # have coefficient_tuple() method that returns the
1253 # coefficients of 1, i, j, and k -- in that order.
1254 return tr
.coefficient_tuple()[0]
1257 class RealMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1261 The identity function, for embedding real matrices into real
1267 def real_unembed(M
):
1269 The identity function, for unembedding real matrices from real
1275 class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra
):
1277 The rank-n simple EJA consisting of real symmetric n-by-n
1278 matrices, the usual symmetric Jordan product, and the trace inner
1279 product. It has dimension `(n^2 + n)/2` over the reals.
1283 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1287 sage: J = RealSymmetricEJA(2)
1288 sage: e0, e1, e2 = J.gens()
1296 In theory, our "field" can be any subfield of the reals::
1298 sage: RealSymmetricEJA(2, RDF)
1299 Euclidean Jordan algebra of dimension 3 over Real Double Field
1300 sage: RealSymmetricEJA(2, RR)
1301 Euclidean Jordan algebra of dimension 3 over Real Field with
1302 53 bits of precision
1306 The dimension of this algebra is `(n^2 + n) / 2`::
1308 sage: set_random_seed()
1309 sage: n_max = RealSymmetricEJA._max_random_instance_size()
1310 sage: n = ZZ.random_element(1, n_max)
1311 sage: J = RealSymmetricEJA(n)
1312 sage: J.dimension() == (n^2 + n)/2
1315 The Jordan multiplication is what we think it is::
1317 sage: set_random_seed()
1318 sage: J = RealSymmetricEJA.random_instance()
1319 sage: x,y = J.random_elements(2)
1320 sage: actual = (x*y).natural_representation()
1321 sage: X = x.natural_representation()
1322 sage: Y = y.natural_representation()
1323 sage: expected = (X*Y + Y*X)/2
1324 sage: actual == expected
1326 sage: J(expected) == x*y
1329 We can change the generator prefix::
1331 sage: RealSymmetricEJA(3, prefix='q').gens()
1332 (q0, q1, q2, q3, q4, q5)
1334 Our natural basis is normalized with respect to the natural inner
1335 product unless we specify otherwise::
1337 sage: set_random_seed()
1338 sage: J = RealSymmetricEJA.random_instance()
1339 sage: all( b.norm() == 1 for b in J.gens() )
1342 Since our natural basis is normalized with respect to the natural
1343 inner product, and since we know that this algebra is an EJA, any
1344 left-multiplication operator's matrix will be symmetric because
1345 natural->EJA basis representation is an isometry and within the EJA
1346 the operator is self-adjoint by the Jordan axiom::
1348 sage: set_random_seed()
1349 sage: x = RealSymmetricEJA.random_instance().random_element()
1350 sage: x.operator().matrix().is_symmetric()
1353 We can construct the (trivial) algebra of rank zero::
1355 sage: RealSymmetricEJA(0)
1356 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1360 def _denormalized_basis(cls
, n
, field
):
1362 Return a basis for the space of real symmetric n-by-n matrices.
1366 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1370 sage: set_random_seed()
1371 sage: n = ZZ.random_element(1,5)
1372 sage: B = RealSymmetricEJA._denormalized_basis(n,QQ)
1373 sage: all( M.is_symmetric() for M in B)
1377 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1381 for j
in range(i
+1):
1382 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
1386 Sij
= Eij
+ Eij
.transpose()
1392 def _max_random_instance_size():
1393 return 4 # Dimension 10
1396 def __init__(self
, n
, field
=AA
, **kwargs
):
1397 basis
= self
._denormalized
_basis
(n
, field
)
1398 super(RealSymmetricEJA
, self
).__init
__(field
,
1402 self
.rank
.set_cache(n
)
1405 class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1409 Embed the n-by-n complex matrix ``M`` into the space of real
1410 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1411 bi` to the block matrix ``[[a,b],[-b,a]]``.
1415 sage: from mjo.eja.eja_algebra import \
1416 ....: ComplexMatrixEuclideanJordanAlgebra
1420 sage: F = QuadraticField(-1, 'I')
1421 sage: x1 = F(4 - 2*i)
1422 sage: x2 = F(1 + 2*i)
1425 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1426 sage: ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1435 Embedding is a homomorphism (isomorphism, in fact)::
1437 sage: set_random_seed()
1438 sage: n_max = ComplexMatrixEuclideanJordanAlgebra._max_random_instance_size()
1439 sage: n = ZZ.random_element(n_max)
1440 sage: F = QuadraticField(-1, 'I')
1441 sage: X = random_matrix(F, n)
1442 sage: Y = random_matrix(F, n)
1443 sage: Xe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X)
1444 sage: Ye = ComplexMatrixEuclideanJordanAlgebra.real_embed(Y)
1445 sage: XYe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1452 raise ValueError("the matrix 'M' must be square")
1454 # We don't need any adjoined elements...
1455 field
= M
.base_ring().base_ring()
1459 a
= z
.list()[0] # real part, I guess
1460 b
= z
.list()[1] # imag part, I guess
1461 blocks
.append(matrix(field
, 2, [[a
,b
],[-b
,a
]]))
1463 return matrix
.block(field
, n
, blocks
)
1467 def real_unembed(M
):
1469 The inverse of _embed_complex_matrix().
1473 sage: from mjo.eja.eja_algebra import \
1474 ....: ComplexMatrixEuclideanJordanAlgebra
1478 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1479 ....: [-2, 1, -4, 3],
1480 ....: [ 9, 10, 11, 12],
1481 ....: [-10, 9, -12, 11] ])
1482 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(A)
1484 [ 10*I + 9 12*I + 11]
1488 Unembedding is the inverse of embedding::
1490 sage: set_random_seed()
1491 sage: F = QuadraticField(-1, 'I')
1492 sage: M = random_matrix(F, 3)
1493 sage: Me = ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1494 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1500 raise ValueError("the matrix 'M' must be square")
1501 if not n
.mod(2).is_zero():
1502 raise ValueError("the matrix 'M' must be a complex embedding")
1504 # If "M" was normalized, its base ring might have roots
1505 # adjoined and they can stick around after unembedding.
1506 field
= M
.base_ring()
1507 R
= PolynomialRing(field
, 'z')
1510 # Sage doesn't know how to embed AA into QQbar, i.e. how
1511 # to adjoin sqrt(-1) to AA.
1514 F
= field
.extension(z
**2 + 1, 'I', embedding
=CLF(-1).sqrt())
1517 # Go top-left to bottom-right (reading order), converting every
1518 # 2-by-2 block we see to a single complex element.
1520 for k
in range(n
/2):
1521 for j
in range(n
/2):
1522 submat
= M
[2*k
:2*k
+2,2*j
:2*j
+2]
1523 if submat
[0,0] != submat
[1,1]:
1524 raise ValueError('bad on-diagonal submatrix')
1525 if submat
[0,1] != -submat
[1,0]:
1526 raise ValueError('bad off-diagonal submatrix')
1527 z
= submat
[0,0] + submat
[0,1]*i
1530 return matrix(F
, n
/2, elements
)
1534 def natural_inner_product(cls
,X
,Y
):
1536 Compute a natural inner product in this algebra directly from
1541 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1545 This gives the same answer as the slow, default method implemented
1546 in :class:`MatrixEuclideanJordanAlgebra`::
1548 sage: set_random_seed()
1549 sage: J = ComplexHermitianEJA.random_instance()
1550 sage: x,y = J.random_elements(2)
1551 sage: Xe = x.natural_representation()
1552 sage: Ye = y.natural_representation()
1553 sage: X = ComplexHermitianEJA.real_unembed(Xe)
1554 sage: Y = ComplexHermitianEJA.real_unembed(Ye)
1555 sage: expected = (X*Y).trace().real()
1556 sage: actual = ComplexHermitianEJA.natural_inner_product(Xe,Ye)
1557 sage: actual == expected
1561 return RealMatrixEuclideanJordanAlgebra
.natural_inner_product(X
,Y
)/2
1564 class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra
):
1566 The rank-n simple EJA consisting of complex Hermitian n-by-n
1567 matrices over the real numbers, the usual symmetric Jordan product,
1568 and the real-part-of-trace inner product. It has dimension `n^2` over
1573 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1577 In theory, our "field" can be any subfield of the reals::
1579 sage: ComplexHermitianEJA(2, RDF)
1580 Euclidean Jordan algebra of dimension 4 over Real Double Field
1581 sage: ComplexHermitianEJA(2, RR)
1582 Euclidean Jordan algebra of dimension 4 over Real Field with
1583 53 bits of precision
1587 The dimension of this algebra is `n^2`::
1589 sage: set_random_seed()
1590 sage: n_max = ComplexHermitianEJA._max_random_instance_size()
1591 sage: n = ZZ.random_element(1, n_max)
1592 sage: J = ComplexHermitianEJA(n)
1593 sage: J.dimension() == n^2
1596 The Jordan multiplication is what we think it is::
1598 sage: set_random_seed()
1599 sage: J = ComplexHermitianEJA.random_instance()
1600 sage: x,y = J.random_elements(2)
1601 sage: actual = (x*y).natural_representation()
1602 sage: X = x.natural_representation()
1603 sage: Y = y.natural_representation()
1604 sage: expected = (X*Y + Y*X)/2
1605 sage: actual == expected
1607 sage: J(expected) == x*y
1610 We can change the generator prefix::
1612 sage: ComplexHermitianEJA(2, prefix='z').gens()
1615 Our natural basis is normalized with respect to the natural inner
1616 product unless we specify otherwise::
1618 sage: set_random_seed()
1619 sage: J = ComplexHermitianEJA.random_instance()
1620 sage: all( b.norm() == 1 for b in J.gens() )
1623 Since our natural basis is normalized with respect to the natural
1624 inner product, and since we know that this algebra is an EJA, any
1625 left-multiplication operator's matrix will be symmetric because
1626 natural->EJA basis representation is an isometry and within the EJA
1627 the operator is self-adjoint by the Jordan axiom::
1629 sage: set_random_seed()
1630 sage: x = ComplexHermitianEJA.random_instance().random_element()
1631 sage: x.operator().matrix().is_symmetric()
1634 We can construct the (trivial) algebra of rank zero::
1636 sage: ComplexHermitianEJA(0)
1637 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1642 def _denormalized_basis(cls
, n
, field
):
1644 Returns a basis for the space of complex Hermitian n-by-n matrices.
1646 Why do we embed these? Basically, because all of numerical linear
1647 algebra assumes that you're working with vectors consisting of `n`
1648 entries from a field and scalars from the same field. There's no way
1649 to tell SageMath that (for example) the vectors contain complex
1650 numbers, while the scalar field is real.
1654 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1658 sage: set_random_seed()
1659 sage: n = ZZ.random_element(1,5)
1660 sage: field = QuadraticField(2, 'sqrt2')
1661 sage: B = ComplexHermitianEJA._denormalized_basis(n, field)
1662 sage: all( M.is_symmetric() for M in B)
1666 R
= PolynomialRing(field
, 'z')
1668 F
= field
.extension(z
**2 + 1, 'I')
1671 # This is like the symmetric case, but we need to be careful:
1673 # * We want conjugate-symmetry, not just symmetry.
1674 # * The diagonal will (as a result) be real.
1678 for j
in range(i
+1):
1679 Eij
= matrix(F
, n
, lambda k
,l
: k
==i
and l
==j
)
1681 Sij
= cls
.real_embed(Eij
)
1684 # The second one has a minus because it's conjugated.
1685 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
1687 Sij_imag
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
1690 # Since we embedded these, we can drop back to the "field" that we
1691 # started with instead of the complex extension "F".
1692 return ( s
.change_ring(field
) for s
in S
)
1695 def __init__(self
, n
, field
=AA
, **kwargs
):
1696 basis
= self
._denormalized
_basis
(n
,field
)
1697 super(ComplexHermitianEJA
,self
).__init
__(field
,
1701 self
.rank
.set_cache(n
)
1704 class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1708 Embed the n-by-n quaternion matrix ``M`` into the space of real
1709 matrices of size 4n-by-4n by first sending each quaternion entry `z
1710 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
1711 c+di],[-c + di, a-bi]]`, and then embedding those into a real
1716 sage: from mjo.eja.eja_algebra import \
1717 ....: QuaternionMatrixEuclideanJordanAlgebra
1721 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1722 sage: i,j,k = Q.gens()
1723 sage: x = 1 + 2*i + 3*j + 4*k
1724 sage: M = matrix(Q, 1, [[x]])
1725 sage: QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1731 Embedding is a homomorphism (isomorphism, in fact)::
1733 sage: set_random_seed()
1734 sage: n_max = QuaternionMatrixEuclideanJordanAlgebra._max_random_instance_size()
1735 sage: n = ZZ.random_element(n_max)
1736 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1737 sage: X = random_matrix(Q, n)
1738 sage: Y = random_matrix(Q, n)
1739 sage: Xe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X)
1740 sage: Ye = QuaternionMatrixEuclideanJordanAlgebra.real_embed(Y)
1741 sage: XYe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1746 quaternions
= M
.base_ring()
1749 raise ValueError("the matrix 'M' must be square")
1751 F
= QuadraticField(-1, 'I')
1756 t
= z
.coefficient_tuple()
1761 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
1762 [-c
+ d
*i
, a
- b
*i
]])
1763 realM
= ComplexMatrixEuclideanJordanAlgebra
.real_embed(cplxM
)
1764 blocks
.append(realM
)
1766 # We should have real entries by now, so use the realest field
1767 # we've got for the return value.
1768 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
1773 def real_unembed(M
):
1775 The inverse of _embed_quaternion_matrix().
1779 sage: from mjo.eja.eja_algebra import \
1780 ....: QuaternionMatrixEuclideanJordanAlgebra
1784 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
1785 ....: [-2, 1, -4, 3],
1786 ....: [-3, 4, 1, -2],
1787 ....: [-4, -3, 2, 1]])
1788 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(M)
1789 [1 + 2*i + 3*j + 4*k]
1793 Unembedding is the inverse of embedding::
1795 sage: set_random_seed()
1796 sage: Q = QuaternionAlgebra(QQ, -1, -1)
1797 sage: M = random_matrix(Q, 3)
1798 sage: Me = QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1799 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1805 raise ValueError("the matrix 'M' must be square")
1806 if not n
.mod(4).is_zero():
1807 raise ValueError("the matrix 'M' must be a quaternion embedding")
1809 # Use the base ring of the matrix to ensure that its entries can be
1810 # multiplied by elements of the quaternion algebra.
1811 field
= M
.base_ring()
1812 Q
= QuaternionAlgebra(field
,-1,-1)
1815 # Go top-left to bottom-right (reading order), converting every
1816 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
1819 for l
in range(n
/4):
1820 for m
in range(n
/4):
1821 submat
= ComplexMatrixEuclideanJordanAlgebra
.real_unembed(
1822 M
[4*l
:4*l
+4,4*m
:4*m
+4] )
1823 if submat
[0,0] != submat
[1,1].conjugate():
1824 raise ValueError('bad on-diagonal submatrix')
1825 if submat
[0,1] != -submat
[1,0].conjugate():
1826 raise ValueError('bad off-diagonal submatrix')
1827 z
= submat
[0,0].real()
1828 z
+= submat
[0,0].imag()*i
1829 z
+= submat
[0,1].real()*j
1830 z
+= submat
[0,1].imag()*k
1833 return matrix(Q
, n
/4, elements
)
1837 def natural_inner_product(cls
,X
,Y
):
1839 Compute a natural inner product in this algebra directly from
1844 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1848 This gives the same answer as the slow, default method implemented
1849 in :class:`MatrixEuclideanJordanAlgebra`::
1851 sage: set_random_seed()
1852 sage: J = QuaternionHermitianEJA.random_instance()
1853 sage: x,y = J.random_elements(2)
1854 sage: Xe = x.natural_representation()
1855 sage: Ye = y.natural_representation()
1856 sage: X = QuaternionHermitianEJA.real_unembed(Xe)
1857 sage: Y = QuaternionHermitianEJA.real_unembed(Ye)
1858 sage: expected = (X*Y).trace().coefficient_tuple()[0]
1859 sage: actual = QuaternionHermitianEJA.natural_inner_product(Xe,Ye)
1860 sage: actual == expected
1864 return RealMatrixEuclideanJordanAlgebra
.natural_inner_product(X
,Y
)/4
1867 class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra
):
1869 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
1870 matrices, the usual symmetric Jordan product, and the
1871 real-part-of-trace inner product. It has dimension `2n^2 - n` over
1876 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1880 In theory, our "field" can be any subfield of the reals::
1882 sage: QuaternionHermitianEJA(2, RDF)
1883 Euclidean Jordan algebra of dimension 6 over Real Double Field
1884 sage: QuaternionHermitianEJA(2, RR)
1885 Euclidean Jordan algebra of dimension 6 over Real Field with
1886 53 bits of precision
1890 The dimension of this algebra is `2*n^2 - n`::
1892 sage: set_random_seed()
1893 sage: n_max = QuaternionHermitianEJA._max_random_instance_size()
1894 sage: n = ZZ.random_element(1, n_max)
1895 sage: J = QuaternionHermitianEJA(n)
1896 sage: J.dimension() == 2*(n^2) - n
1899 The Jordan multiplication is what we think it is::
1901 sage: set_random_seed()
1902 sage: J = QuaternionHermitianEJA.random_instance()
1903 sage: x,y = J.random_elements(2)
1904 sage: actual = (x*y).natural_representation()
1905 sage: X = x.natural_representation()
1906 sage: Y = y.natural_representation()
1907 sage: expected = (X*Y + Y*X)/2
1908 sage: actual == expected
1910 sage: J(expected) == x*y
1913 We can change the generator prefix::
1915 sage: QuaternionHermitianEJA(2, prefix='a').gens()
1916 (a0, a1, a2, a3, a4, a5)
1918 Our natural basis is normalized with respect to the natural inner
1919 product unless we specify otherwise::
1921 sage: set_random_seed()
1922 sage: J = QuaternionHermitianEJA.random_instance()
1923 sage: all( b.norm() == 1 for b in J.gens() )
1926 Since our natural basis is normalized with respect to the natural
1927 inner product, and since we know that this algebra is an EJA, any
1928 left-multiplication operator's matrix will be symmetric because
1929 natural->EJA basis representation is an isometry and within the EJA
1930 the operator is self-adjoint by the Jordan axiom::
1932 sage: set_random_seed()
1933 sage: x = QuaternionHermitianEJA.random_instance().random_element()
1934 sage: x.operator().matrix().is_symmetric()
1937 We can construct the (trivial) algebra of rank zero::
1939 sage: QuaternionHermitianEJA(0)
1940 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1944 def _denormalized_basis(cls
, n
, field
):
1946 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
1948 Why do we embed these? Basically, because all of numerical
1949 linear algebra assumes that you're working with vectors consisting
1950 of `n` entries from a field and scalars from the same field. There's
1951 no way to tell SageMath that (for example) the vectors contain
1952 complex numbers, while the scalar field is real.
1956 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1960 sage: set_random_seed()
1961 sage: n = ZZ.random_element(1,5)
1962 sage: B = QuaternionHermitianEJA._denormalized_basis(n,QQ)
1963 sage: all( M.is_symmetric() for M in B )
1967 Q
= QuaternionAlgebra(QQ
,-1,-1)
1970 # This is like the symmetric case, but we need to be careful:
1972 # * We want conjugate-symmetry, not just symmetry.
1973 # * The diagonal will (as a result) be real.
1977 for j
in range(i
+1):
1978 Eij
= matrix(Q
, n
, lambda k
,l
: k
==i
and l
==j
)
1980 Sij
= cls
.real_embed(Eij
)
1983 # The second, third, and fourth ones have a minus
1984 # because they're conjugated.
1985 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
1987 Sij_I
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
1989 Sij_J
= cls
.real_embed(J
*Eij
- J
*Eij
.transpose())
1991 Sij_K
= cls
.real_embed(K
*Eij
- K
*Eij
.transpose())
1994 # Since we embedded these, we can drop back to the "field" that we
1995 # started with instead of the quaternion algebra "Q".
1996 return ( s
.change_ring(field
) for s
in S
)
1999 def __init__(self
, n
, field
=AA
, **kwargs
):
2000 basis
= self
._denormalized
_basis
(n
,field
)
2001 super(QuaternionHermitianEJA
,self
).__init
__(field
,
2005 self
.rank
.set_cache(n
)
2008 class HadamardEJA(RationalBasisEuclideanJordanAlgebra
):
2010 Return the Euclidean Jordan Algebra corresponding to the set
2011 `R^n` under the Hadamard product.
2013 Note: this is nothing more than the Cartesian product of ``n``
2014 copies of the spin algebra. Once Cartesian product algebras
2015 are implemented, this can go.
2019 sage: from mjo.eja.eja_algebra import HadamardEJA
2023 This multiplication table can be verified by hand::
2025 sage: J = HadamardEJA(3)
2026 sage: e0,e1,e2 = J.gens()
2042 We can change the generator prefix::
2044 sage: HadamardEJA(3, prefix='r').gens()
2048 def __init__(self
, n
, field
=AA
, **kwargs
):
2049 V
= VectorSpace(field
, n
)
2050 mult_table
= [ [ V
.gen(i
)*(i
== j
) for j
in range(n
) ]
2053 super(HadamardEJA
, self
).__init
__(field
,
2057 self
.rank
.set_cache(n
)
2060 def _max_random_instance_size():
2063 def inner_product(self
, x
, y
):
2065 Faster to reimplement than to use natural representations.
2069 sage: from mjo.eja.eja_algebra import HadamardEJA
2073 Ensure that this is the usual inner product for the algebras
2076 sage: set_random_seed()
2077 sage: J = HadamardEJA.random_instance()
2078 sage: x,y = J.random_elements(2)
2079 sage: X = x.natural_representation()
2080 sage: Y = y.natural_representation()
2081 sage: x.inner_product(y) == J.natural_inner_product(X,Y)
2085 return x
.to_vector().inner_product(y
.to_vector())
2088 class BilinearFormEJA(RationalBasisEuclideanJordanAlgebra
):
2090 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2091 with the half-trace inner product and jordan product ``x*y =
2092 (x0*y0 + <B*x_bar,y_bar>, x0*y_bar + y0*x_bar)`` where ``B`` is a
2093 symmetric positive-definite "bilinear form" matrix. It has
2094 dimension `n` over the reals, and reduces to the ``JordanSpinEJA``
2095 when ``B`` is the identity matrix of order ``n-1``.
2099 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
2100 ....: JordanSpinEJA)
2104 When no bilinear form is specified, the identity matrix is used,
2105 and the resulting algebra is the Jordan spin algebra::
2107 sage: J0 = BilinearFormEJA(3)
2108 sage: J1 = JordanSpinEJA(3)
2109 sage: J0.multiplication_table() == J0.multiplication_table()
2114 We can create a zero-dimensional algebra::
2116 sage: J = BilinearFormEJA(0)
2120 We can check the multiplication condition given in the Jordan, von
2121 Neumann, and Wigner paper (and also discussed on my "On the
2122 symmetry..." paper). Note that this relies heavily on the standard
2123 choice of basis, as does anything utilizing the bilinear form matrix::
2125 sage: set_random_seed()
2126 sage: n = ZZ.random_element(5)
2127 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2128 sage: B = M.transpose()*M
2129 sage: J = BilinearFormEJA(n, B=B)
2130 sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
2131 sage: V = J.vector_space()
2132 sage: sis = [ J.from_vector(V([0] + (M.inverse()*ei).list()))
2133 ....: for ei in eis ]
2134 sage: actual = [ sis[i]*sis[j]
2135 ....: for i in range(n-1)
2136 ....: for j in range(n-1) ]
2137 sage: expected = [ J.one() if i == j else J.zero()
2138 ....: for i in range(n-1)
2139 ....: for j in range(n-1) ]
2140 sage: actual == expected
2143 def __init__(self
, n
, field
=AA
, B
=None, **kwargs
):
2145 self
._B
= matrix
.identity(field
, max(0,n
-1))
2149 V
= VectorSpace(field
, n
)
2150 mult_table
= [[V
.zero() for j
in range(n
)] for i
in range(n
)]
2159 z0
= x0
*y0
+ (self
._B
*xbar
).inner_product(ybar
)
2160 zbar
= y0
*xbar
+ x0
*ybar
2161 z
= V([z0
] + zbar
.list())
2162 mult_table
[i
][j
] = z
2164 # The rank of this algebra is two, unless we're in a
2165 # one-dimensional ambient space (because the rank is bounded
2166 # by the ambient dimension).
2167 super(BilinearFormEJA
, self
).__init
__(field
,
2171 self
.rank
.set_cache(min(n
,2))
2174 def _max_random_instance_size():
2177 def inner_product(self
, x
, y
):
2179 Half of the trace inner product.
2181 This is defined so that the special case of the Jordan spin
2182 algebra gets the usual inner product.
2186 sage: from mjo.eja.eja_algebra import BilinearFormEJA
2190 Ensure that this is one-half of the trace inner-product when
2191 the algebra isn't just the reals (when ``n`` isn't one). This
2192 is in Faraut and Koranyi, and also my "On the symmetry..."
2195 sage: set_random_seed()
2196 sage: n = ZZ.random_element(2,5)
2197 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2198 sage: B = M.transpose()*M
2199 sage: J = BilinearFormEJA(n, B=B)
2200 sage: x = J.random_element()
2201 sage: y = J.random_element()
2202 sage: x.inner_product(y) == (x*y).trace()/2
2206 xvec
= x
.to_vector()
2208 yvec
= y
.to_vector()
2210 return x
[0]*y
[0] + (self
._B
*xbar
).inner_product(ybar
)
2213 class JordanSpinEJA(BilinearFormEJA
):
2215 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2216 with the usual inner product and jordan product ``x*y =
2217 (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2222 sage: from mjo.eja.eja_algebra import JordanSpinEJA
2226 This multiplication table can be verified by hand::
2228 sage: J = JordanSpinEJA(4)
2229 sage: e0,e1,e2,e3 = J.gens()
2245 We can change the generator prefix::
2247 sage: JordanSpinEJA(2, prefix='B').gens()
2252 Ensure that we have the usual inner product on `R^n`::
2254 sage: set_random_seed()
2255 sage: J = JordanSpinEJA.random_instance()
2256 sage: x,y = J.random_elements(2)
2257 sage: X = x.natural_representation()
2258 sage: Y = y.natural_representation()
2259 sage: x.inner_product(y) == J.natural_inner_product(X,Y)
2263 def __init__(self
, n
, field
=AA
, **kwargs
):
2264 # This is a special case of the BilinearFormEJA with the identity
2265 # matrix as its bilinear form.
2266 super(JordanSpinEJA
, self
).__init
__(n
, field
, **kwargs
)
2269 class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra
):
2271 The trivial Euclidean Jordan algebra consisting of only a zero element.
2275 sage: from mjo.eja.eja_algebra import TrivialEJA
2279 sage: J = TrivialEJA()
2286 sage: 7*J.one()*12*J.one()
2288 sage: J.one().inner_product(J.one())
2290 sage: J.one().norm()
2292 sage: J.one().subalgebra_generated_by()
2293 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2298 def __init__(self
, field
=AA
, **kwargs
):
2300 super(TrivialEJA
, self
).__init
__(field
,
2304 # The rank is zero using my definition, namely the dimension of the
2305 # largest subalgebra generated by any element.
2306 self
.rank
.set_cache(0)
2309 class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra
):
2311 The external (orthogonal) direct sum of two other Euclidean Jordan
2312 algebras. Essentially the Cartesian product of its two factors.
2313 Every Euclidean Jordan algebra decomposes into an orthogonal
2314 direct sum of simple Euclidean Jordan algebras, so no generality
2315 is lost by providing only this construction.
2319 sage: from mjo.eja.eja_algebra import (HadamardEJA,
2320 ....: RealSymmetricEJA,
2325 sage: J1 = HadamardEJA(2)
2326 sage: J2 = RealSymmetricEJA(3)
2327 sage: J = DirectSumEJA(J1,J2)
2334 def __init__(self
, J1
, J2
, field
=AA
, **kwargs
):
2335 self
._factors
= (J1
, J2
)
2339 V
= VectorSpace(field
, n
)
2340 mult_table
= [ [ V
.zero() for j
in range(n
) ]
2344 p
= (J1
.monomial(i
)*J1
.monomial(j
)).to_vector()
2345 mult_table
[i
][j
] = V(p
.list() + [field
.zero()]*n2
)
2349 p
= (J2
.monomial(i
)*J2
.monomial(j
)).to_vector()
2350 mult_table
[n1
+i
][n1
+j
] = V([field
.zero()]*n1
+ p
.list())
2352 super(DirectSumEJA
, self
).__init
__(field
,
2356 self
.rank
.set_cache(J1
.rank() + J2
.rank())
2361 Return the pair of this algebra's factors.
2365 sage: from mjo.eja.eja_algebra import (HadamardEJA,
2366 ....: JordanSpinEJA,
2371 sage: J1 = HadamardEJA(2,QQ)
2372 sage: J2 = JordanSpinEJA(3,QQ)
2373 sage: J = DirectSumEJA(J1,J2)
2375 (Euclidean Jordan algebra of dimension 2 over Rational Field,
2376 Euclidean Jordan algebra of dimension 3 over Rational Field)
2379 return self
._factors
2381 def projections(self
):
2383 Return a pair of projections onto this algebra's factors.
2387 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
2388 ....: ComplexHermitianEJA,
2393 sage: J1 = JordanSpinEJA(2)
2394 sage: J2 = ComplexHermitianEJA(2)
2395 sage: J = DirectSumEJA(J1,J2)
2396 sage: (pi_left, pi_right) = J.projections()
2397 sage: J.one().to_vector()
2399 sage: pi_left(J.one()).to_vector()
2401 sage: pi_right(J.one()).to_vector()
2405 (J1
,J2
) = self
.factors()
2407 pi_left
= lambda x
: J1
.from_vector(x
.to_vector()[:n
])
2408 pi_right
= lambda x
: J2
.from_vector(x
.to_vector()[n
:])
2409 return (pi_left
, pi_right
)
2411 def inclusions(self
):
2413 Return the pair of inclusion maps from our factors into us.
2417 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
2418 ....: RealSymmetricEJA,
2423 sage: J1 = JordanSpinEJA(3)
2424 sage: J2 = RealSymmetricEJA(2)
2425 sage: J = DirectSumEJA(J1,J2)
2426 sage: (iota_left, iota_right) = J.inclusions()
2427 sage: iota_left(J1.zero()) == J.zero()
2429 sage: iota_right(J2.zero()) == J.zero()
2431 sage: J1.one().to_vector()
2433 sage: iota_left(J1.one()).to_vector()
2435 sage: J2.one().to_vector()
2437 sage: iota_right(J2.one()).to_vector()
2439 sage: J.one().to_vector()
2443 (J1
,J2
) = self
.factors()
2445 V_basis
= self
.vector_space().basis()
2446 I1
= matrix
.column(self
.base_ring(), V_basis
[:n
])
2447 I2
= matrix
.column(self
.base_ring(), V_basis
[n
:])
2448 iota_left
= lambda x
: self
.from_vector(I1
*x
.to_vector())
2449 iota_right
= lambda x
: self
.from_vector(I2
*+x
.to_vector())
2450 return (iota_left
, iota_right
)
2452 def inner_product(self
, x
, y
):
2454 The standard Cartesian inner-product.
2456 We project ``x`` and ``y`` onto our factors, and add up the
2457 inner-products from the subalgebras.
2462 sage: from mjo.eja.eja_algebra import (HadamardEJA,
2463 ....: QuaternionHermitianEJA,
2468 sage: J1 = HadamardEJA(3)
2469 sage: J2 = QuaternionHermitianEJA(2,QQ,normalize_basis=False)
2470 sage: J = DirectSumEJA(J1,J2)
2475 sage: x1.inner_product(x2)
2477 sage: y1.inner_product(y2)
2479 sage: J.one().inner_product(J.one())
2483 (pi_left
, pi_right
) = self
.projections()
2489 return (x1
.inner_product(y1
) + x2
.inner_product(y2
))