2 Euclidean Jordan Algebras. These are formally-real Jordan Algebras;
3 specifically those where u^2 + v^2 = 0 implies that u = v = 0. They
4 are used in optimization, and have some additional nice methods beyond
5 what can be supported in a general Jordan Algebra.
10 sage: from mjo.eja.eja_algebra import random_eja
15 Euclidean Jordan algebra of dimension...
19 from itertools
import repeat
21 from sage
.algebras
.quatalg
.quaternion_algebra
import QuaternionAlgebra
22 from sage
.categories
.magmatic_algebras
import MagmaticAlgebras
23 from sage
.combinat
.free_module
import CombinatorialFreeModule
24 from sage
.matrix
.constructor
import matrix
25 from sage
.matrix
.matrix_space
import MatrixSpace
26 from sage
.misc
.cachefunc
import cached_method
27 from sage
.misc
.table
import table
28 from sage
.modules
.free_module
import FreeModule
, VectorSpace
29 from sage
.rings
.all
import (ZZ
, QQ
, AA
, QQbar
, RR
, RLF
, CLF
,
32 from mjo
.eja
.eja_element
import FiniteDimensionalEuclideanJordanAlgebraElement
33 from mjo
.eja
.eja_operator
import FiniteDimensionalEuclideanJordanAlgebraOperator
34 from mjo
.eja
.eja_utils
import _mat2vec
36 class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule
):
38 The lowest-level class for representing a Euclidean Jordan algebra.
40 def _coerce_map_from_base_ring(self
):
42 Disable the map from the base ring into the algebra.
44 Performing a nonsense conversion like this automatically
45 is counterpedagogical. The fallback is to try the usual
46 element constructor, which should also fail.
50 sage: from mjo.eja.eja_algebra import random_eja
54 sage: set_random_seed()
55 sage: J = random_eja()
57 Traceback (most recent call last):
59 ValueError: not an element of this algebra
76 * field -- the scalar field for this algebra (must be real)
78 * multiplication_table -- the multiplication table for this
79 algebra's implicit basis. Only the lower-triangular portion
80 of the table is used, since the multiplication is assumed
85 sage: from mjo.eja.eja_algebra import (
86 ....: FiniteDimensionalEuclideanJordanAlgebra,
92 By definition, Jordan multiplication commutes::
94 sage: set_random_seed()
95 sage: J = random_eja()
96 sage: x,y = J.random_elements(2)
100 An error is raised if the Jordan product is not commutative::
102 sage: JP = ((1,2),(0,0))
103 sage: IP = ((1,0),(0,1))
104 sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,JP,IP)
105 Traceback (most recent call last):
107 ValueError: Jordan product is not commutative
109 An error is raised if the inner-product is not commutative::
111 sage: JP = ((1,0),(0,1))
112 sage: IP = ((1,2),(0,0))
113 sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,JP,IP)
114 Traceback (most recent call last):
116 ValueError: inner-product is not commutative
120 The ``field`` we're given must be real with ``check_field=True``::
122 sage: JordanSpinEJA(2, field=QQbar)
123 Traceback (most recent call last):
125 ValueError: scalar field is not real
126 sage: JordanSpinEJA(2, field=QQbar, check_field=False)
127 Euclidean Jordan algebra of dimension 2 over Algebraic Field
129 The multiplication table must be square with ``check_axioms=True``::
131 sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,((),()),((1,),))
132 Traceback (most recent call last):
134 ValueError: multiplication table is not square
136 The multiplication and inner-product tables must be the same
137 size (and in particular, the inner-product table must also be
138 square) with ``check_axioms=True``::
140 sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,((1,),),(()))
141 Traceback (most recent call last):
143 ValueError: multiplication and inner-product tables are
145 sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,((1,),),((1,2),))
146 Traceback (most recent call last):
148 ValueError: multiplication and inner-product tables are
153 if not field
.is_subring(RR
):
154 # Note: this does return true for the real algebraic
155 # field, the rationals, and any quadratic field where
156 # we've specified a real embedding.
157 raise ValueError("scalar field is not real")
160 # The multiplication and inner-product tables should be square
161 # if the user wants us to verify them. And we verify them as
162 # soon as possible, because we want to exploit their symmetry.
163 n
= len(multiplication_table
)
165 if not all( len(l
) == n
for l
in multiplication_table
):
166 raise ValueError("multiplication table is not square")
168 # If the multiplication table is square, we can check if
169 # the inner-product table is square by comparing it to the
170 # multiplication table's dimensions.
171 msg
= "multiplication and inner-product tables are different sizes"
172 if not len(inner_product_table
) == n
:
173 raise ValueError(msg
)
175 if not all( len(l
) == n
for l
in inner_product_table
):
176 raise ValueError(msg
)
178 # Check commutativity of the Jordan product (symmetry of
179 # the multiplication table) and the commutativity of the
180 # inner-product (symmetry of the inner-product table)
181 # first if we're going to check them at all.. This has to
182 # be done before we define product_on_basis(), because
183 # that method assumes that self._multiplication_table is
184 # symmetric. And it has to be done before we build
185 # self._inner_product_matrix, because the process used to
186 # construct it assumes symmetry as well.
187 if not all( multiplication_table
[j
][i
]
188 == multiplication_table
[i
][j
]
190 for j
in range(i
+1) ):
191 raise ValueError("Jordan product is not commutative")
193 if not all( inner_product_table
[j
][i
]
194 == inner_product_table
[i
][j
]
196 for j
in range(i
+1) ):
197 raise ValueError("inner-product is not commutative")
199 self
._matrix
_basis
= matrix_basis
202 category
= MagmaticAlgebras(field
).FiniteDimensional()
203 category
= category
.WithBasis().Unital()
205 fda
= super(FiniteDimensionalEuclideanJordanAlgebra
, self
)
210 self
.print_options(bracket
='')
212 # The multiplication table we're given is necessarily in terms
213 # of vectors, because we don't have an algebra yet for
214 # anything to be an element of. However, it's faster in the
215 # long run to have the multiplication table be in terms of
216 # algebra elements. We do this after calling the superclass
217 # constructor so that from_vector() knows what to do.
219 # Note: we take advantage of symmetry here, and only store
220 # the lower-triangular portion of the table.
221 self
._multiplication
_table
= [ [ self
.vector_space().zero()
222 for j
in range(i
+1) ]
227 elt
= self
.from_vector(multiplication_table
[i
][j
])
228 self
._multiplication
_table
[i
][j
] = elt
230 self
._multiplication
_table
= tuple(map(tuple, self
._multiplication
_table
))
232 # Save our inner product as a matrix, since the efficiency of
233 # matrix multiplication will usually outweigh the fact that we
234 # have to store a redundant upper- or lower-triangular part.
235 # Pre-cache the fact that these are Hermitian (real symmetric,
236 # in fact) in case some e.g. matrix multiplication routine can
237 # take advantage of it.
238 ip_matrix_constructor
= lambda i
,j
: inner_product_table
[i
][j
] if j
<= i
else inner_product_table
[j
][i
]
239 self
._inner
_product
_matrix
= matrix(field
, n
, ip_matrix_constructor
)
240 self
._inner
_product
_matrix
._cache
= {'hermitian': True}
241 self
._inner
_product
_matrix
.set_immutable()
244 if not self
._is
_jordanian
():
245 raise ValueError("Jordan identity does not hold")
246 if not self
._inner
_product
_is
_associative
():
247 raise ValueError("inner product is not associative")
249 def _element_constructor_(self
, elt
):
251 Construct an element of this algebra from its vector or matrix
254 This gets called only after the parent element _call_ method
255 fails to find a coercion for the argument.
259 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
261 ....: RealSymmetricEJA)
265 The identity in `S^n` is converted to the identity in the EJA::
267 sage: J = RealSymmetricEJA(3)
268 sage: I = matrix.identity(QQ,3)
269 sage: J(I) == J.one()
272 This skew-symmetric matrix can't be represented in the EJA::
274 sage: J = RealSymmetricEJA(3)
275 sage: A = matrix(QQ,3, lambda i,j: i-j)
277 Traceback (most recent call last):
279 ValueError: not an element of this algebra
283 Ensure that we can convert any element of the two non-matrix
284 simple algebras (whose matrix representations are columns)
285 back and forth faithfully::
287 sage: set_random_seed()
288 sage: J = HadamardEJA.random_instance()
289 sage: x = J.random_element()
290 sage: J(x.to_vector().column()) == x
292 sage: J = JordanSpinEJA.random_instance()
293 sage: x = J.random_element()
294 sage: J(x.to_vector().column()) == x
298 msg
= "not an element of this algebra"
300 # The superclass implementation of random_element()
301 # needs to be able to coerce "0" into the algebra.
303 elif elt
in self
.base_ring():
304 # Ensure that no base ring -> algebra coercion is performed
305 # by this method. There's some stupidity in sage that would
306 # otherwise propagate to this method; for example, sage thinks
307 # that the integer 3 belongs to the space of 2-by-2 matrices.
308 raise ValueError(msg
)
310 if elt
not in self
.matrix_space():
311 raise ValueError(msg
)
313 # Thanks for nothing! Matrix spaces aren't vector spaces in
314 # Sage, so we have to figure out its matrix-basis coordinates
315 # ourselves. We use the basis space's ring instead of the
316 # element's ring because the basis space might be an algebraic
317 # closure whereas the base ring of the 3-by-3 identity matrix
318 # could be QQ instead of QQbar.
320 # We pass check=False because the matrix basis is "guaranteed"
321 # to be linearly independent... right? Ha ha.
322 V
= VectorSpace(self
.base_ring(), elt
.nrows()*elt
.ncols())
323 W
= V
.span_of_basis( (_mat2vec(s
) for s
in self
.matrix_basis()),
327 coords
= W
.coordinate_vector(_mat2vec(elt
))
328 except ArithmeticError: # vector is not in free module
329 raise ValueError(msg
)
331 return self
.from_vector(coords
)
335 Return a string representation of ``self``.
339 sage: from mjo.eja.eja_algebra import JordanSpinEJA
343 Ensure that it says what we think it says::
345 sage: JordanSpinEJA(2, field=AA)
346 Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
347 sage: JordanSpinEJA(3, field=RDF)
348 Euclidean Jordan algebra of dimension 3 over Real Double Field
351 fmt
= "Euclidean Jordan algebra of dimension {} over {}"
352 return fmt
.format(self
.dimension(), self
.base_ring())
354 def product_on_basis(self
, i
, j
):
355 # We only stored the lower-triangular portion of the
356 # multiplication table.
358 return self
._multiplication
_table
[i
][j
]
360 return self
._multiplication
_table
[j
][i
]
362 def _is_commutative(self
):
364 Whether or not this algebra's multiplication table is commutative.
366 This method should of course always return ``True``, unless
367 this algebra was constructed with ``check_axioms=False`` and
368 passed an invalid multiplication table.
370 return all( self
.product_on_basis(i
,j
) == self
.product_on_basis(i
,j
)
371 for i
in range(self
.dimension())
372 for j
in range(self
.dimension()) )
374 def _is_jordanian(self
):
376 Whether or not this algebra's multiplication table respects the
377 Jordan identity `(x^{2})(xy) = x(x^{2}y)`.
379 We only check one arrangement of `x` and `y`, so for a
380 ``True`` result to be truly true, you should also check
381 :meth:`_is_commutative`. This method should of course always
382 return ``True``, unless this algebra was constructed with
383 ``check_axioms=False`` and passed an invalid multiplication table.
385 return all( (self
.monomial(i
)**2)*(self
.monomial(i
)*self
.monomial(j
))
387 (self
.monomial(i
))*((self
.monomial(i
)**2)*self
.monomial(j
))
388 for i
in range(self
.dimension())
389 for j
in range(self
.dimension()) )
391 def _inner_product_is_associative(self
):
393 Return whether or not this algebra's inner product `B` is
394 associative; that is, whether or not `B(xy,z) = B(x,yz)`.
396 This method should of course always return ``True``, unless
397 this algebra was constructed with ``check_axioms=False`` and
398 passed an invalid multiplication table.
401 # Used to check whether or not something is zero in an inexact
402 # ring. This number is sufficient to allow the construction of
403 # QuaternionHermitianEJA(2, field=RDF) with check_axioms=True.
406 for i
in range(self
.dimension()):
407 for j
in range(self
.dimension()):
408 for k
in range(self
.dimension()):
412 diff
= (x
*y
).inner_product(z
) - x
.inner_product(y
*z
)
414 if self
.base_ring().is_exact():
418 if diff
.abs() > epsilon
:
424 def characteristic_polynomial_of(self
):
426 Return the algebra's "characteristic polynomial of" function,
427 which is itself a multivariate polynomial that, when evaluated
428 at the coordinates of some algebra element, returns that
429 element's characteristic polynomial.
431 The resulting polynomial has `n+1` variables, where `n` is the
432 dimension of this algebra. The first `n` variables correspond to
433 the coordinates of an algebra element: when evaluated at the
434 coordinates of an algebra element with respect to a certain
435 basis, the result is a univariate polynomial (in the one
436 remaining variable ``t``), namely the characteristic polynomial
441 sage: from mjo.eja.eja_algebra import JordanSpinEJA, TrivialEJA
445 The characteristic polynomial in the spin algebra is given in
446 Alizadeh, Example 11.11::
448 sage: J = JordanSpinEJA(3)
449 sage: p = J.characteristic_polynomial_of(); p
450 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
451 sage: xvec = J.one().to_vector()
455 By definition, the characteristic polynomial is a monic
456 degree-zero polynomial in a rank-zero algebra. Note that
457 Cayley-Hamilton is indeed satisfied since the polynomial
458 ``1`` evaluates to the identity element of the algebra on
461 sage: J = TrivialEJA()
462 sage: J.characteristic_polynomial_of()
469 # The list of coefficient polynomials a_0, a_1, a_2, ..., a_(r-1).
470 a
= self
._charpoly
_coefficients
()
472 # We go to a bit of trouble here to reorder the
473 # indeterminates, so that it's easier to evaluate the
474 # characteristic polynomial at x's coordinates and get back
475 # something in terms of t, which is what we want.
476 S
= PolynomialRing(self
.base_ring(),'t')
480 S
= PolynomialRing(S
, R
.variable_names())
483 return (t
**r
+ sum( a
[k
]*(t
**k
) for k
in range(r
) ))
485 def coordinate_polynomial_ring(self
):
487 The multivariate polynomial ring in which this algebra's
488 :meth:`characteristic_polynomial_of` lives.
492 sage: from mjo.eja.eja_algebra import (HadamardEJA,
493 ....: RealSymmetricEJA)
497 sage: J = HadamardEJA(2)
498 sage: J.coordinate_polynomial_ring()
499 Multivariate Polynomial Ring in X1, X2...
500 sage: J = RealSymmetricEJA(3,field=QQ,orthonormalize=False)
501 sage: J.coordinate_polynomial_ring()
502 Multivariate Polynomial Ring in X1, X2, X3, X4, X5, X6...
505 var_names
= tuple( "X%d" % z
for z
in range(1, self
.dimension()+1) )
506 return PolynomialRing(self
.base_ring(), var_names
)
508 def inner_product(self
, x
, y
):
510 The inner product associated with this Euclidean Jordan algebra.
512 Defaults to the trace inner product, but can be overridden by
513 subclasses if they are sure that the necessary properties are
518 sage: from mjo.eja.eja_algebra import (random_eja,
520 ....: BilinearFormEJA)
524 Our inner product is "associative," which means the following for
525 a symmetric bilinear form::
527 sage: set_random_seed()
528 sage: J = random_eja()
529 sage: x,y,z = J.random_elements(3)
530 sage: (x*y).inner_product(z) == y.inner_product(x*z)
535 Ensure that this is the usual inner product for the algebras
538 sage: set_random_seed()
539 sage: J = HadamardEJA.random_instance()
540 sage: x,y = J.random_elements(2)
541 sage: actual = x.inner_product(y)
542 sage: expected = x.to_vector().inner_product(y.to_vector())
543 sage: actual == expected
546 Ensure that this is one-half of the trace inner-product in a
547 BilinearFormEJA that isn't just the reals (when ``n`` isn't
548 one). This is in Faraut and Koranyi, and also my "On the
551 sage: set_random_seed()
552 sage: J = BilinearFormEJA.random_instance()
553 sage: n = J.dimension()
554 sage: x = J.random_element()
555 sage: y = J.random_element()
556 sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
559 B
= self
._inner
_product
_matrix
560 return (B
*x
.to_vector()).inner_product(y
.to_vector())
563 def is_trivial(self
):
565 Return whether or not this algebra is trivial.
567 A trivial algebra contains only the zero element.
571 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
576 sage: J = ComplexHermitianEJA(3)
582 sage: J = TrivialEJA()
587 return self
.dimension() == 0
590 def multiplication_table(self
):
592 Return a visual representation of this algebra's multiplication
593 table (on basis elements).
597 sage: from mjo.eja.eja_algebra import JordanSpinEJA
601 sage: J = JordanSpinEJA(4)
602 sage: J.multiplication_table()
603 +----++----+----+----+----+
604 | * || e0 | e1 | e2 | e3 |
605 +====++====+====+====+====+
606 | e0 || e0 | e1 | e2 | e3 |
607 +----++----+----+----+----+
608 | e1 || e1 | e0 | 0 | 0 |
609 +----++----+----+----+----+
610 | e2 || e2 | 0 | e0 | 0 |
611 +----++----+----+----+----+
612 | e3 || e3 | 0 | 0 | e0 |
613 +----++----+----+----+----+
617 M
= [ [ self
.zero() for j
in range(n
) ]
621 M
[i
][j
] = self
._multiplication
_table
[i
][j
]
625 # Prepend the left "header" column entry Can't do this in
626 # the loop because it messes up the symmetry.
627 M
[i
] = [self
.monomial(i
)] + M
[i
]
629 # Prepend the header row.
630 M
= [["*"] + list(self
.gens())] + M
631 return table(M
, header_row
=True, header_column
=True, frame
=True)
634 def matrix_basis(self
):
636 Return an (often more natural) representation of this algebras
637 basis as an ordered tuple of matrices.
639 Every finite-dimensional Euclidean Jordan Algebra is a, up to
640 Jordan isomorphism, a direct sum of five simple
641 algebras---four of which comprise Hermitian matrices. And the
642 last type of algebra can of course be thought of as `n`-by-`1`
643 column matrices (ambiguusly called column vectors) to avoid
644 special cases. As a result, matrices (and column vectors) are
645 a natural representation format for Euclidean Jordan algebra
648 But, when we construct an algebra from a basis of matrices,
649 those matrix representations are lost in favor of coordinate
650 vectors *with respect to* that basis. We could eventually
651 convert back if we tried hard enough, but having the original
652 representations handy is valuable enough that we simply store
653 them and return them from this method.
655 Why implement this for non-matrix algebras? Avoiding special
656 cases for the :class:`BilinearFormEJA` pays with simplicity in
657 its own right. But mainly, we would like to be able to assume
658 that elements of a :class:`DirectSumEJA` can be displayed
659 nicely, without having to have special classes for direct sums
660 one of whose components was a matrix algebra.
664 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
665 ....: RealSymmetricEJA)
669 sage: J = RealSymmetricEJA(2)
671 Finite family {0: e0, 1: e1, 2: e2}
672 sage: J.matrix_basis()
674 [1 0] [ 0 0.7071067811865475?] [0 0]
675 [0 0], [0.7071067811865475? 0], [0 1]
680 sage: J = JordanSpinEJA(2)
682 Finite family {0: e0, 1: e1}
683 sage: J.matrix_basis()
689 if self
._matrix
_basis
is None:
690 M
= self
.matrix_space()
691 return tuple( M(b
.to_vector()) for b
in self
.basis() )
693 return self
._matrix
_basis
696 def matrix_space(self
):
698 Return the matrix space in which this algebra's elements live, if
699 we think of them as matrices (including column vectors of the
702 Generally this will be an `n`-by-`1` column-vector space,
703 except when the algebra is trivial. There it's `n`-by-`n`
704 (where `n` is zero), to ensure that two elements of the matrix
705 space (empty matrices) can be multiplied.
707 Matrix algebras override this with something more useful.
709 if self
.is_trivial():
710 return MatrixSpace(self
.base_ring(), 0)
711 elif self
._matrix
_basis
is None or len(self
._matrix
_basis
) == 0:
712 return MatrixSpace(self
.base_ring(), self
.dimension(), 1)
714 return self
._matrix
_basis
[0].matrix_space()
720 Return the unit element of this algebra.
724 sage: from mjo.eja.eja_algebra import (HadamardEJA,
729 sage: J = HadamardEJA(5)
731 e0 + e1 + e2 + e3 + e4
735 The identity element acts like the identity::
737 sage: set_random_seed()
738 sage: J = random_eja()
739 sage: x = J.random_element()
740 sage: J.one()*x == x and x*J.one() == x
743 The matrix of the unit element's operator is the identity::
745 sage: set_random_seed()
746 sage: J = random_eja()
747 sage: actual = J.one().operator().matrix()
748 sage: expected = matrix.identity(J.base_ring(), J.dimension())
749 sage: actual == expected
752 Ensure that the cached unit element (often precomputed by
753 hand) agrees with the computed one::
755 sage: set_random_seed()
756 sage: J = random_eja()
757 sage: cached = J.one()
758 sage: J.one.clear_cache()
759 sage: J.one() == cached
763 # We can brute-force compute the matrices of the operators
764 # that correspond to the basis elements of this algebra.
765 # If some linear combination of those basis elements is the
766 # algebra identity, then the same linear combination of
767 # their matrices has to be the identity matrix.
769 # Of course, matrices aren't vectors in sage, so we have to
770 # appeal to the "long vectors" isometry.
771 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
773 # Now we use basic linear algebra to find the coefficients,
774 # of the matrices-as-vectors-linear-combination, which should
775 # work for the original algebra basis too.
776 A
= matrix(self
.base_ring(), oper_vecs
)
778 # We used the isometry on the left-hand side already, but we
779 # still need to do it for the right-hand side. Recall that we
780 # wanted something that summed to the identity matrix.
781 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
783 # Now if there's an identity element in the algebra, this
784 # should work. We solve on the left to avoid having to
785 # transpose the matrix "A".
786 return self
.from_vector(A
.solve_left(b
))
789 def peirce_decomposition(self
, c
):
791 The Peirce decomposition of this algebra relative to the
794 In the future, this can be extended to a complete system of
795 orthogonal idempotents.
799 - ``c`` -- an idempotent of this algebra.
803 A triple (J0, J5, J1) containing two subalgebras and one subspace
806 - ``J0`` -- the algebra on the eigenspace of ``c.operator()``
807 corresponding to the eigenvalue zero.
809 - ``J5`` -- the eigenspace (NOT a subalgebra) of ``c.operator()``
810 corresponding to the eigenvalue one-half.
812 - ``J1`` -- the algebra on the eigenspace of ``c.operator()``
813 corresponding to the eigenvalue one.
815 These are the only possible eigenspaces for that operator, and this
816 algebra is a direct sum of them. The spaces ``J0`` and ``J1`` are
817 orthogonal, and are subalgebras of this algebra with the appropriate
822 sage: from mjo.eja.eja_algebra import random_eja, RealSymmetricEJA
826 The canonical example comes from the symmetric matrices, which
827 decompose into diagonal and off-diagonal parts::
829 sage: J = RealSymmetricEJA(3)
830 sage: C = matrix(QQ, [ [1,0,0],
834 sage: J0,J5,J1 = J.peirce_decomposition(c)
836 Euclidean Jordan algebra of dimension 1...
838 Vector space of degree 6 and dimension 2...
840 Euclidean Jordan algebra of dimension 3...
841 sage: J0.one().to_matrix()
845 sage: orig_df = AA.options.display_format
846 sage: AA.options.display_format = 'radical'
847 sage: J.from_vector(J5.basis()[0]).to_matrix()
851 sage: J.from_vector(J5.basis()[1]).to_matrix()
855 sage: AA.options.display_format = orig_df
856 sage: J1.one().to_matrix()
863 Every algebra decomposes trivially with respect to its identity
866 sage: set_random_seed()
867 sage: J = random_eja()
868 sage: J0,J5,J1 = J.peirce_decomposition(J.one())
869 sage: J0.dimension() == 0 and J5.dimension() == 0
871 sage: J1.superalgebra() == J and J1.dimension() == J.dimension()
874 The decomposition is into eigenspaces, and its components are
875 therefore necessarily orthogonal. Moreover, the identity
876 elements in the two subalgebras are the projections onto their
877 respective subspaces of the superalgebra's identity element::
879 sage: set_random_seed()
880 sage: J = random_eja()
881 sage: x = J.random_element()
882 sage: if not J.is_trivial():
883 ....: while x.is_nilpotent():
884 ....: x = J.random_element()
885 sage: c = x.subalgebra_idempotent()
886 sage: J0,J5,J1 = J.peirce_decomposition(c)
888 sage: for (w,y,z) in zip(J0.basis(), J5.basis(), J1.basis()):
889 ....: w = w.superalgebra_element()
890 ....: y = J.from_vector(y)
891 ....: z = z.superalgebra_element()
892 ....: ipsum += w.inner_product(y).abs()
893 ....: ipsum += w.inner_product(z).abs()
894 ....: ipsum += y.inner_product(z).abs()
897 sage: J1(c) == J1.one()
899 sage: J0(J.one() - c) == J0.one()
903 if not c
.is_idempotent():
904 raise ValueError("element is not idempotent: %s" % c
)
906 from mjo
.eja
.eja_subalgebra
import FiniteDimensionalEuclideanJordanSubalgebra
908 # Default these to what they should be if they turn out to be
909 # trivial, because eigenspaces_left() won't return eigenvalues
910 # corresponding to trivial spaces (e.g. it returns only the
911 # eigenspace corresponding to lambda=1 if you take the
912 # decomposition relative to the identity element).
913 trivial
= FiniteDimensionalEuclideanJordanSubalgebra(self
, ())
914 J0
= trivial
# eigenvalue zero
915 J5
= VectorSpace(self
.base_ring(), 0) # eigenvalue one-half
916 J1
= trivial
# eigenvalue one
918 for (eigval
, eigspace
) in c
.operator().matrix().right_eigenspaces():
919 if eigval
== ~
(self
.base_ring()(2)):
922 gens
= tuple( self
.from_vector(b
) for b
in eigspace
.basis() )
923 subalg
= FiniteDimensionalEuclideanJordanSubalgebra(self
,
931 raise ValueError("unexpected eigenvalue: %s" % eigval
)
936 def random_element(self
, thorough
=False):
938 Return a random element of this algebra.
940 Our algebra superclass method only returns a linear
941 combination of at most two basis elements. We instead
942 want the vector space "random element" method that
943 returns a more diverse selection.
947 - ``thorough`` -- (boolean; default False) whether or not we
948 should generate irrational coefficients for the random
949 element when our base ring is irrational; this slows the
950 algebra operations to a crawl, but any truly random method
954 # For a general base ring... maybe we can trust this to do the
955 # right thing? Unlikely, but.
956 V
= self
.vector_space()
957 v
= V
.random_element()
959 if self
.base_ring() is AA
:
960 # The "random element" method of the algebraic reals is
961 # stupid at the moment, and only returns integers between
962 # -2 and 2, inclusive:
964 # https://trac.sagemath.org/ticket/30875
966 # Instead, we implement our own "random vector" method,
967 # and then coerce that into the algebra. We use the vector
968 # space degree here instead of the dimension because a
969 # subalgebra could (for example) be spanned by only two
970 # vectors, each with five coordinates. We need to
971 # generate all five coordinates.
973 v
*= QQbar
.random_element().real()
975 v
*= QQ
.random_element()
977 return self
.from_vector(V
.coordinate_vector(v
))
979 def random_elements(self
, count
, thorough
=False):
981 Return ``count`` random elements as a tuple.
985 - ``thorough`` -- (boolean; default False) whether or not we
986 should generate irrational coefficients for the random
987 elements when our base ring is irrational; this slows the
988 algebra operations to a crawl, but any truly random method
993 sage: from mjo.eja.eja_algebra import JordanSpinEJA
997 sage: J = JordanSpinEJA(3)
998 sage: x,y,z = J.random_elements(3)
999 sage: all( [ x in J, y in J, z in J ])
1001 sage: len( J.random_elements(10) ) == 10
1005 return tuple( self
.random_element(thorough
)
1006 for idx
in range(count
) )
1010 def _charpoly_coefficients(self
):
1012 The `r` polynomial coefficients of the "characteristic polynomial
1015 n
= self
.dimension()
1016 R
= self
.coordinate_polynomial_ring()
1018 F
= R
.fraction_field()
1021 # From a result in my book, these are the entries of the
1022 # basis representation of L_x.
1023 return sum( vars[k
]*self
.monomial(k
).operator().matrix()[i
,j
]
1026 L_x
= matrix(F
, n
, n
, L_x_i_j
)
1029 if self
.rank
.is_in_cache():
1031 # There's no need to pad the system with redundant
1032 # columns if we *know* they'll be redundant.
1035 # Compute an extra power in case the rank is equal to
1036 # the dimension (otherwise, we would stop at x^(r-1)).
1037 x_powers
= [ (L_x
**k
)*self
.one().to_vector()
1038 for k
in range(n
+1) ]
1039 A
= matrix
.column(F
, x_powers
[:n
])
1040 AE
= A
.extended_echelon_form()
1047 # The theory says that only the first "r" coefficients are
1048 # nonzero, and they actually live in the original polynomial
1049 # ring and not the fraction field. We negate them because
1050 # in the actual characteristic polynomial, they get moved
1051 # to the other side where x^r lives.
1052 return -A_rref
.solve_right(E
*b
).change_ring(R
)[:r
]
1057 Return the rank of this EJA.
1059 This is a cached method because we know the rank a priori for
1060 all of the algebras we can construct. Thus we can avoid the
1061 expensive ``_charpoly_coefficients()`` call unless we truly
1062 need to compute the whole characteristic polynomial.
1066 sage: from mjo.eja.eja_algebra import (HadamardEJA,
1067 ....: JordanSpinEJA,
1068 ....: RealSymmetricEJA,
1069 ....: ComplexHermitianEJA,
1070 ....: QuaternionHermitianEJA,
1075 The rank of the Jordan spin algebra is always two::
1077 sage: JordanSpinEJA(2).rank()
1079 sage: JordanSpinEJA(3).rank()
1081 sage: JordanSpinEJA(4).rank()
1084 The rank of the `n`-by-`n` Hermitian real, complex, or
1085 quaternion matrices is `n`::
1087 sage: RealSymmetricEJA(4).rank()
1089 sage: ComplexHermitianEJA(3).rank()
1091 sage: QuaternionHermitianEJA(2).rank()
1096 Ensure that every EJA that we know how to construct has a
1097 positive integer rank, unless the algebra is trivial in
1098 which case its rank will be zero::
1100 sage: set_random_seed()
1101 sage: J = random_eja()
1105 sage: r > 0 or (r == 0 and J.is_trivial())
1108 Ensure that computing the rank actually works, since the ranks
1109 of all simple algebras are known and will be cached by default::
1111 sage: set_random_seed() # long time
1112 sage: J = random_eja() # long time
1113 sage: caches = J.rank() # long time
1114 sage: J.rank.clear_cache() # long time
1115 sage: J.rank() == cached # long time
1119 return len(self
._charpoly
_coefficients
())
1122 def vector_space(self
):
1124 Return the vector space that underlies this algebra.
1128 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1132 sage: J = RealSymmetricEJA(2)
1133 sage: J.vector_space()
1134 Vector space of dimension 3 over...
1137 return self
.zero().to_vector().parent().ambient_vector_space()
1140 Element
= FiniteDimensionalEuclideanJordanAlgebraElement
1142 class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra
):
1144 New class for algebras whose supplied basis elements have all rational entries.
1148 sage: from mjo.eja.eja_algebra import BilinearFormEJA
1152 The supplied basis is orthonormalized by default::
1154 sage: B = matrix(QQ, [[1, 0, 0], [0, 25, -32], [0, -32, 41]])
1155 sage: J = BilinearFormEJA(B)
1156 sage: J.matrix_basis()
1169 orthonormalize
=True,
1176 # Abuse the check_field parameter to check that the entries of
1177 # out basis (in ambient coordinates) are in the field QQ.
1178 if not all( all(b_i
in QQ
for b_i
in b
.list()) for b
in basis
):
1179 raise TypeError("basis not rational")
1181 # Temporary(?) hack to ensure that the matrix and vector bases
1182 # are over the same ring.
1183 basis
= tuple( b
.change_ring(field
) for b
in basis
)
1186 vector_basis
= basis
1188 from sage
.structure
.element
import is_Matrix
1189 basis_is_matrices
= False
1193 if is_Matrix(basis
[0]):
1194 basis_is_matrices
= True
1195 from mjo
.eja
.eja_utils
import _vec2mat
1196 vector_basis
= tuple( map(_mat2vec
,basis
) )
1197 degree
= basis
[0].nrows()**2
1199 degree
= basis
[0].degree()
1201 V
= VectorSpace(field
, degree
)
1203 # If we were asked to orthonormalize, and if the orthonormal
1204 # basis is different from the given one, then we also want to
1205 # compute multiplication and inner-product tables for the
1206 # deorthonormalized basis. These can be used later to
1207 # construct a deorthonormalized copy of this algebra over QQ
1208 # in which several operations are much faster.
1209 self
._rational
_algebra
= None
1212 if self
.base_ring() is not QQ
:
1213 # There's no point in constructing the extra algebra if this
1214 # one is already rational. If the original basis is rational
1215 # but normalization would make it irrational, then this whole
1216 # constructor will just fail anyway as it tries to stick an
1217 # irrational number into a rational algebra.
1219 # Note: the same Jordan and inner-products work here,
1220 # because they are necessarily defined with respect to
1221 # ambient coordinates and not any particular basis.
1222 self
._rational
_algebra
= RationalBasisEuclideanJordanAlgebra(
1227 orthonormalize
=False,
1233 # Compute the deorthonormalized tables before we orthonormalize
1234 # the given basis. The "check" parameter here guarantees that
1235 # the basis is linearly-independent.
1236 W
= V
.span_of_basis( vector_basis
, check
=check_axioms
)
1238 # Note: the Jordan and inner-products are defined in terms
1239 # of the ambient basis. It's important that their arguments
1240 # are in ambient coordinates as well.
1242 for j
in range(i
+1):
1243 # given basis w.r.t. ambient coords
1244 q_i
= vector_basis
[i
]
1245 q_j
= vector_basis
[j
]
1247 if basis_is_matrices
:
1251 elt
= jordan_product(q_i
, q_j
)
1252 ip
= inner_product(q_i
, q_j
)
1254 if basis_is_matrices
:
1255 # do another mat2vec because the multiplication
1256 # table is in terms of vectors
1259 # We overwrite the name "vector_basis" in a second, but never modify it
1260 # in place, to this effectively makes a copy of it.
1261 deortho_vector_basis
= vector_basis
1262 self
._deortho
_matrix
= None
1265 from mjo
.eja
.eja_utils
import gram_schmidt
1266 if basis_is_matrices
:
1267 vector_ip
= lambda x
,y
: inner_product(_vec2mat(x
), _vec2mat(y
))
1268 vector_basis
= gram_schmidt(vector_basis
, vector_ip
)
1270 vector_basis
= gram_schmidt(vector_basis
, inner_product
)
1272 # Normalize the "matrix" basis, too!
1273 basis
= vector_basis
1275 if basis_is_matrices
:
1276 basis
= tuple( map(_vec2mat
,basis
) )
1278 W
= V
.span_of_basis( vector_basis
, check
=check_axioms
)
1280 # Now "W" is the vector space of our algebra coordinates. The
1281 # variables "X1", "X2",... refer to the entries of vectors in
1282 # W. Thus to convert back and forth between the orthonormal
1283 # coordinates and the given ones, we need to stick the original
1285 U
= V
.span_of_basis( deortho_vector_basis
, check
=check_axioms
)
1286 self
._deortho
_matrix
= matrix( U
.coordinate_vector(q
)
1287 for q
in vector_basis
)
1289 # If the superclass constructor is going to verify the
1290 # symmetry of this table, it has better at least be
1293 mult_table
= [ [0 for j
in range(n
)] for i
in range(n
) ]
1294 ip_table
= [ [0 for j
in range(n
)] for i
in range(n
) ]
1296 mult_table
= [ [0 for j
in range(i
+1)] for i
in range(n
) ]
1297 ip_table
= [ [0 for j
in range(i
+1)] for i
in range(n
) ]
1299 # Note: the Jordan and inner-products are defined in terms
1300 # of the ambient basis. It's important that their arguments
1301 # are in ambient coordinates as well.
1303 for j
in range(i
+1):
1304 # ortho basis w.r.t. ambient coords
1305 q_i
= vector_basis
[i
]
1306 q_j
= vector_basis
[j
]
1308 if basis_is_matrices
:
1312 elt
= jordan_product(q_i
, q_j
)
1313 ip
= inner_product(q_i
, q_j
)
1315 if basis_is_matrices
:
1316 # do another mat2vec because the multiplication
1317 # table is in terms of vectors
1320 elt
= W
.coordinate_vector(elt
)
1321 mult_table
[i
][j
] = elt
1324 # The tables are square if we're verifying that they
1326 mult_table
[j
][i
] = elt
1329 if basis_is_matrices
:
1333 basis
= tuple( x
.column() for x
in basis
)
1335 super().__init
__(field
,
1340 basis
, # matrix basis
1345 def _charpoly_coefficients(self
):
1349 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
1350 ....: JordanSpinEJA)
1354 The base ring of the resulting polynomial coefficients is what
1355 it should be, and not the rationals (unless the algebra was
1356 already over the rationals)::
1358 sage: J = JordanSpinEJA(3)
1359 sage: J._charpoly_coefficients()
1360 (X1^2 - X2^2 - X3^2, -2*X1)
1361 sage: a0 = J._charpoly_coefficients()[0]
1363 Algebraic Real Field
1364 sage: a0.base_ring()
1365 Algebraic Real Field
1368 if self
.base_ring() is QQ
or self
._rational
_algebra
is None:
1369 # There's no need to construct *another* algebra over the
1370 # rationals if this one is already over the
1371 # rationals. Likewise, if we never orthonormalized our
1372 # basis, we might as well just use the given one.
1373 superclass
= super(RationalBasisEuclideanJordanAlgebra
, self
)
1374 return superclass
._charpoly
_coefficients
()
1376 # Do the computation over the rationals. The answer will be
1377 # the same, because all we've done is a change of basis.
1378 # Then, change back from QQ to our real base ring
1379 a
= ( a_i
.change_ring(self
.base_ring())
1380 for a_i
in self
._rational
_algebra
._charpoly
_coefficients
() )
1382 # Now convert the coordinate variables back to the
1383 # deorthonormalized ones.
1384 R
= self
.coordinate_polynomial_ring()
1385 from sage
.modules
.free_module_element
import vector
1386 X
= vector(R
, R
.gens())
1387 BX
= self
._deortho
_matrix
*X
1389 subs_dict
= { X[i]: BX[i] for i in range(len(X)) }
1390 return tuple( a_i
.subs(subs_dict
) for a_i
in a
)
1392 class ConcreteEuclideanJordanAlgebra(RationalBasisEuclideanJordanAlgebra
):
1394 A class for the Euclidean Jordan algebras that we know by name.
1396 These are the Jordan algebras whose basis, multiplication table,
1397 rank, and so on are known a priori. More to the point, they are
1398 the Euclidean Jordan algebras for which we are able to conjure up
1399 a "random instance."
1403 sage: from mjo.eja.eja_algebra import ConcreteEuclideanJordanAlgebra
1407 Our basis is normalized with respect to the algebra's inner
1408 product, unless we specify otherwise::
1410 sage: set_random_seed()
1411 sage: J = ConcreteEuclideanJordanAlgebra.random_instance()
1412 sage: all( b.norm() == 1 for b in J.gens() )
1415 Since our basis is orthonormal with respect to the algebra's inner
1416 product, and since we know that this algebra is an EJA, any
1417 left-multiplication operator's matrix will be symmetric because
1418 natural->EJA basis representation is an isometry and within the
1419 EJA the operator is self-adjoint by the Jordan axiom::
1421 sage: set_random_seed()
1422 sage: J = ConcreteEuclideanJordanAlgebra.random_instance()
1423 sage: x = J.random_element()
1424 sage: x.operator().is_self_adjoint()
1429 def _max_random_instance_size():
1431 Return an integer "size" that is an upper bound on the size of
1432 this algebra when it is used in a random test
1433 case. Unfortunately, the term "size" is ambiguous -- when
1434 dealing with `R^n` under either the Hadamard or Jordan spin
1435 product, the "size" refers to the dimension `n`. When dealing
1436 with a matrix algebra (real symmetric or complex/quaternion
1437 Hermitian), it refers to the size of the matrix, which is far
1438 less than the dimension of the underlying vector space.
1440 This method must be implemented in each subclass.
1442 raise NotImplementedError
1445 def random_instance(cls
, *args
, **kwargs
):
1447 Return a random instance of this type of algebra.
1449 This method should be implemented in each subclass.
1451 from sage
.misc
.prandom
import choice
1452 eja_class
= choice(cls
.__subclasses
__())
1454 # These all bubble up to the RationalBasisEuclideanJordanAlgebra
1455 # superclass constructor, so any (kw)args valid there are also
1457 return eja_class
.random_instance(*args
, **kwargs
)
1460 class MatrixEuclideanJordanAlgebra
:
1462 def dimension_over_reals():
1464 The dimension of this matrix's base ring over the reals.
1466 The reals are dimension one over themselves, obviously; that's
1467 just `\mathbb{R}^{1}`. Likewise, the complex numbers `a + bi`
1468 have dimension two. Finally, the quaternions have dimension
1469 four over the reals.
1471 This is used to determine the size of the matrix returned from
1472 :meth:`real_embed`, among other things.
1474 raise NotImplementedError
1477 def real_embed(cls
,M
):
1479 Embed the matrix ``M`` into a space of real matrices.
1481 The matrix ``M`` can have entries in any field at the moment:
1482 the real numbers, complex numbers, or quaternions. And although
1483 they are not a field, we can probably support octonions at some
1484 point, too. This function returns a real matrix that "acts like"
1485 the original with respect to matrix multiplication; i.e.
1487 real_embed(M*N) = real_embed(M)*real_embed(N)
1490 if M
.ncols() != M
.nrows():
1491 raise ValueError("the matrix 'M' must be square")
1496 def real_unembed(cls
,M
):
1498 The inverse of :meth:`real_embed`.
1500 if M
.ncols() != M
.nrows():
1501 raise ValueError("the matrix 'M' must be square")
1502 if not ZZ(M
.nrows()).mod(cls
.dimension_over_reals()).is_zero():
1503 raise ValueError("the matrix 'M' must be a real embedding")
1507 def jordan_product(X
,Y
):
1508 return (X
*Y
+ Y
*X
)/2
1511 def trace_inner_product(cls
,X
,Y
):
1513 Compute the trace inner-product of two real-embeddings.
1517 sage: from mjo.eja.eja_algebra import (RealSymmetricEJA,
1518 ....: ComplexHermitianEJA,
1519 ....: QuaternionHermitianEJA)
1523 This gives the same answer as it would if we computed the trace
1524 from the unembedded (original) matrices::
1526 sage: set_random_seed()
1527 sage: J = RealSymmetricEJA.random_instance()
1528 sage: x,y = J.random_elements(2)
1529 sage: Xe = x.to_matrix()
1530 sage: Ye = y.to_matrix()
1531 sage: X = J.real_unembed(Xe)
1532 sage: Y = J.real_unembed(Ye)
1533 sage: expected = (X*Y).trace()
1534 sage: actual = J.trace_inner_product(Xe,Ye)
1535 sage: actual == expected
1540 sage: set_random_seed()
1541 sage: J = ComplexHermitianEJA.random_instance()
1542 sage: x,y = J.random_elements(2)
1543 sage: Xe = x.to_matrix()
1544 sage: Ye = y.to_matrix()
1545 sage: X = J.real_unembed(Xe)
1546 sage: Y = J.real_unembed(Ye)
1547 sage: expected = (X*Y).trace().real()
1548 sage: actual = J.trace_inner_product(Xe,Ye)
1549 sage: actual == expected
1554 sage: set_random_seed()
1555 sage: J = QuaternionHermitianEJA.random_instance()
1556 sage: x,y = J.random_elements(2)
1557 sage: Xe = x.to_matrix()
1558 sage: Ye = y.to_matrix()
1559 sage: X = J.real_unembed(Xe)
1560 sage: Y = J.real_unembed(Ye)
1561 sage: expected = (X*Y).trace().coefficient_tuple()[0]
1562 sage: actual = J.trace_inner_product(Xe,Ye)
1563 sage: actual == expected
1567 Xu
= cls
.real_unembed(X
)
1568 Yu
= cls
.real_unembed(Y
)
1569 tr
= (Xu
*Yu
).trace()
1572 # Works in QQ, AA, RDF, et cetera.
1574 except AttributeError:
1575 # A quaternion doesn't have a real() method, but does
1576 # have coefficient_tuple() method that returns the
1577 # coefficients of 1, i, j, and k -- in that order.
1578 return tr
.coefficient_tuple()[0]
1581 class RealMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1583 def dimension_over_reals():
1587 class RealSymmetricEJA(ConcreteEuclideanJordanAlgebra
,
1588 RealMatrixEuclideanJordanAlgebra
):
1590 The rank-n simple EJA consisting of real symmetric n-by-n
1591 matrices, the usual symmetric Jordan product, and the trace inner
1592 product. It has dimension `(n^2 + n)/2` over the reals.
1596 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1600 sage: J = RealSymmetricEJA(2)
1601 sage: e0, e1, e2 = J.gens()
1609 In theory, our "field" can be any subfield of the reals::
1611 sage: RealSymmetricEJA(2, field=RDF)
1612 Euclidean Jordan algebra of dimension 3 over Real Double Field
1613 sage: RealSymmetricEJA(2, field=RR)
1614 Euclidean Jordan algebra of dimension 3 over Real Field with
1615 53 bits of precision
1619 The dimension of this algebra is `(n^2 + n) / 2`::
1621 sage: set_random_seed()
1622 sage: n_max = RealSymmetricEJA._max_random_instance_size()
1623 sage: n = ZZ.random_element(1, n_max)
1624 sage: J = RealSymmetricEJA(n)
1625 sage: J.dimension() == (n^2 + n)/2
1628 The Jordan multiplication is what we think it is::
1630 sage: set_random_seed()
1631 sage: J = RealSymmetricEJA.random_instance()
1632 sage: x,y = J.random_elements(2)
1633 sage: actual = (x*y).to_matrix()
1634 sage: X = x.to_matrix()
1635 sage: Y = y.to_matrix()
1636 sage: expected = (X*Y + Y*X)/2
1637 sage: actual == expected
1639 sage: J(expected) == x*y
1642 We can change the generator prefix::
1644 sage: RealSymmetricEJA(3, prefix='q').gens()
1645 (q0, q1, q2, q3, q4, q5)
1647 We can construct the (trivial) algebra of rank zero::
1649 sage: RealSymmetricEJA(0)
1650 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1654 def _denormalized_basis(cls
, n
):
1656 Return a basis for the space of real symmetric n-by-n matrices.
1660 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1664 sage: set_random_seed()
1665 sage: n = ZZ.random_element(1,5)
1666 sage: B = RealSymmetricEJA._denormalized_basis(n)
1667 sage: all( M.is_symmetric() for M in B)
1671 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1675 for j
in range(i
+1):
1676 Eij
= matrix(ZZ
, n
, lambda k
,l
: k
==i
and l
==j
)
1680 Sij
= Eij
+ Eij
.transpose()
1686 def _max_random_instance_size():
1687 return 4 # Dimension 10
1690 def random_instance(cls
, **kwargs
):
1692 Return a random instance of this type of algebra.
1694 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1695 return cls(n
, **kwargs
)
1697 def __init__(self
, n
, **kwargs
):
1698 # We know this is a valid EJA, but will double-check
1699 # if the user passes check_axioms=True.
1700 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
1702 super(RealSymmetricEJA
, self
).__init
__(self
._denormalized
_basis
(n
),
1703 self
.jordan_product
,
1704 self
.trace_inner_product
,
1707 # TODO: this could be factored out somehow, but is left here
1708 # because the MatrixEuclideanJordanAlgebra is not presently
1709 # a subclass of the FDEJA class that defines rank() and one().
1710 self
.rank
.set_cache(n
)
1711 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
1712 self
.one
.set_cache(self(idV
))
1716 class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1718 def dimension_over_reals():
1722 def real_embed(cls
,M
):
1724 Embed the n-by-n complex matrix ``M`` into the space of real
1725 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1726 bi` to the block matrix ``[[a,b],[-b,a]]``.
1730 sage: from mjo.eja.eja_algebra import \
1731 ....: ComplexMatrixEuclideanJordanAlgebra
1735 sage: F = QuadraticField(-1, 'I')
1736 sage: x1 = F(4 - 2*i)
1737 sage: x2 = F(1 + 2*i)
1740 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1741 sage: ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1750 Embedding is a homomorphism (isomorphism, in fact)::
1752 sage: set_random_seed()
1753 sage: n = ZZ.random_element(3)
1754 sage: F = QuadraticField(-1, 'I')
1755 sage: X = random_matrix(F, n)
1756 sage: Y = random_matrix(F, n)
1757 sage: Xe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X)
1758 sage: Ye = ComplexMatrixEuclideanJordanAlgebra.real_embed(Y)
1759 sage: XYe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1764 super(ComplexMatrixEuclideanJordanAlgebra
,cls
).real_embed(M
)
1767 # We don't need any adjoined elements...
1768 field
= M
.base_ring().base_ring()
1772 a
= z
.list()[0] # real part, I guess
1773 b
= z
.list()[1] # imag part, I guess
1774 blocks
.append(matrix(field
, 2, [[a
,b
],[-b
,a
]]))
1776 return matrix
.block(field
, n
, blocks
)
1780 def real_unembed(cls
,M
):
1782 The inverse of _embed_complex_matrix().
1786 sage: from mjo.eja.eja_algebra import \
1787 ....: ComplexMatrixEuclideanJordanAlgebra
1791 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1792 ....: [-2, 1, -4, 3],
1793 ....: [ 9, 10, 11, 12],
1794 ....: [-10, 9, -12, 11] ])
1795 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(A)
1797 [ 10*I + 9 12*I + 11]
1801 Unembedding is the inverse of embedding::
1803 sage: set_random_seed()
1804 sage: F = QuadraticField(-1, 'I')
1805 sage: M = random_matrix(F, 3)
1806 sage: Me = ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1807 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1811 super(ComplexMatrixEuclideanJordanAlgebra
,cls
).real_unembed(M
)
1813 d
= cls
.dimension_over_reals()
1815 # If "M" was normalized, its base ring might have roots
1816 # adjoined and they can stick around after unembedding.
1817 field
= M
.base_ring()
1818 R
= PolynomialRing(field
, 'z')
1821 # Sage doesn't know how to adjoin the complex "i" (the root of
1822 # x^2 + 1) to a field in a general way. Here, we just enumerate
1823 # all of the cases that I have cared to support so far.
1825 # Sage doesn't know how to embed AA into QQbar, i.e. how
1826 # to adjoin sqrt(-1) to AA.
1828 elif not field
.is_exact():
1830 F
= field
.complex_field()
1832 # Works for QQ and... maybe some other fields.
1833 F
= field
.extension(z
**2 + 1, 'I', embedding
=CLF(-1).sqrt())
1836 # Go top-left to bottom-right (reading order), converting every
1837 # 2-by-2 block we see to a single complex element.
1839 for k
in range(n
/d
):
1840 for j
in range(n
/d
):
1841 submat
= M
[d
*k
:d
*k
+d
,d
*j
:d
*j
+d
]
1842 if submat
[0,0] != submat
[1,1]:
1843 raise ValueError('bad on-diagonal submatrix')
1844 if submat
[0,1] != -submat
[1,0]:
1845 raise ValueError('bad off-diagonal submatrix')
1846 z
= submat
[0,0] + submat
[0,1]*i
1849 return matrix(F
, n
/d
, elements
)
1852 class ComplexHermitianEJA(ConcreteEuclideanJordanAlgebra
,
1853 ComplexMatrixEuclideanJordanAlgebra
):
1855 The rank-n simple EJA consisting of complex Hermitian n-by-n
1856 matrices over the real numbers, the usual symmetric Jordan product,
1857 and the real-part-of-trace inner product. It has dimension `n^2` over
1862 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1866 In theory, our "field" can be any subfield of the reals::
1868 sage: ComplexHermitianEJA(2, field=RDF)
1869 Euclidean Jordan algebra of dimension 4 over Real Double Field
1870 sage: ComplexHermitianEJA(2, field=RR)
1871 Euclidean Jordan algebra of dimension 4 over Real Field with
1872 53 bits of precision
1876 The dimension of this algebra is `n^2`::
1878 sage: set_random_seed()
1879 sage: n_max = ComplexHermitianEJA._max_random_instance_size()
1880 sage: n = ZZ.random_element(1, n_max)
1881 sage: J = ComplexHermitianEJA(n)
1882 sage: J.dimension() == n^2
1885 The Jordan multiplication is what we think it is::
1887 sage: set_random_seed()
1888 sage: J = ComplexHermitianEJA.random_instance()
1889 sage: x,y = J.random_elements(2)
1890 sage: actual = (x*y).to_matrix()
1891 sage: X = x.to_matrix()
1892 sage: Y = y.to_matrix()
1893 sage: expected = (X*Y + Y*X)/2
1894 sage: actual == expected
1896 sage: J(expected) == x*y
1899 We can change the generator prefix::
1901 sage: ComplexHermitianEJA(2, prefix='z').gens()
1904 We can construct the (trivial) algebra of rank zero::
1906 sage: ComplexHermitianEJA(0)
1907 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1912 def _denormalized_basis(cls
, n
):
1914 Returns a basis for the space of complex Hermitian n-by-n matrices.
1916 Why do we embed these? Basically, because all of numerical linear
1917 algebra assumes that you're working with vectors consisting of `n`
1918 entries from a field and scalars from the same field. There's no way
1919 to tell SageMath that (for example) the vectors contain complex
1920 numbers, while the scalar field is real.
1924 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1928 sage: set_random_seed()
1929 sage: n = ZZ.random_element(1,5)
1930 sage: field = QuadraticField(2, 'sqrt2')
1931 sage: B = ComplexHermitianEJA._denormalized_basis(n)
1932 sage: all( M.is_symmetric() for M in B)
1937 R
= PolynomialRing(field
, 'z')
1939 F
= field
.extension(z
**2 + 1, 'I')
1942 # This is like the symmetric case, but we need to be careful:
1944 # * We want conjugate-symmetry, not just symmetry.
1945 # * The diagonal will (as a result) be real.
1949 for j
in range(i
+1):
1950 Eij
= matrix(F
, n
, lambda k
,l
: k
==i
and l
==j
)
1952 Sij
= cls
.real_embed(Eij
)
1955 # The second one has a minus because it's conjugated.
1956 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
1958 Sij_imag
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
1961 # Since we embedded these, we can drop back to the "field" that we
1962 # started with instead of the complex extension "F".
1963 return tuple( s
.change_ring(field
) for s
in S
)
1966 def __init__(self
, n
, **kwargs
):
1967 # We know this is a valid EJA, but will double-check
1968 # if the user passes check_axioms=True.
1969 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
1971 super(ComplexHermitianEJA
, self
).__init
__(self
._denormalized
_basis
(n
),
1972 self
.jordan_product
,
1973 self
.trace_inner_product
,
1975 # TODO: this could be factored out somehow, but is left here
1976 # because the MatrixEuclideanJordanAlgebra is not presently
1977 # a subclass of the FDEJA class that defines rank() and one().
1978 self
.rank
.set_cache(n
)
1979 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
1980 self
.one
.set_cache(self(idV
))
1983 def _max_random_instance_size():
1984 return 3 # Dimension 9
1987 def random_instance(cls
, **kwargs
):
1989 Return a random instance of this type of algebra.
1991 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1992 return cls(n
, **kwargs
)
1994 class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1996 def dimension_over_reals():
2000 def real_embed(cls
,M
):
2002 Embed the n-by-n quaternion matrix ``M`` into the space of real
2003 matrices of size 4n-by-4n by first sending each quaternion entry `z
2004 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
2005 c+di],[-c + di, a-bi]]`, and then embedding those into a real
2010 sage: from mjo.eja.eja_algebra import \
2011 ....: QuaternionMatrixEuclideanJordanAlgebra
2015 sage: Q = QuaternionAlgebra(QQ,-1,-1)
2016 sage: i,j,k = Q.gens()
2017 sage: x = 1 + 2*i + 3*j + 4*k
2018 sage: M = matrix(Q, 1, [[x]])
2019 sage: QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
2025 Embedding is a homomorphism (isomorphism, in fact)::
2027 sage: set_random_seed()
2028 sage: n = ZZ.random_element(2)
2029 sage: Q = QuaternionAlgebra(QQ,-1,-1)
2030 sage: X = random_matrix(Q, n)
2031 sage: Y = random_matrix(Q, n)
2032 sage: Xe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X)
2033 sage: Ye = QuaternionMatrixEuclideanJordanAlgebra.real_embed(Y)
2034 sage: XYe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X*Y)
2039 super(QuaternionMatrixEuclideanJordanAlgebra
,cls
).real_embed(M
)
2040 quaternions
= M
.base_ring()
2043 F
= QuadraticField(-1, 'I')
2048 t
= z
.coefficient_tuple()
2053 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
2054 [-c
+ d
*i
, a
- b
*i
]])
2055 realM
= ComplexMatrixEuclideanJordanAlgebra
.real_embed(cplxM
)
2056 blocks
.append(realM
)
2058 # We should have real entries by now, so use the realest field
2059 # we've got for the return value.
2060 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
2065 def real_unembed(cls
,M
):
2067 The inverse of _embed_quaternion_matrix().
2071 sage: from mjo.eja.eja_algebra import \
2072 ....: QuaternionMatrixEuclideanJordanAlgebra
2076 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
2077 ....: [-2, 1, -4, 3],
2078 ....: [-3, 4, 1, -2],
2079 ....: [-4, -3, 2, 1]])
2080 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(M)
2081 [1 + 2*i + 3*j + 4*k]
2085 Unembedding is the inverse of embedding::
2087 sage: set_random_seed()
2088 sage: Q = QuaternionAlgebra(QQ, -1, -1)
2089 sage: M = random_matrix(Q, 3)
2090 sage: Me = QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
2091 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
2095 super(QuaternionMatrixEuclideanJordanAlgebra
,cls
).real_unembed(M
)
2097 d
= cls
.dimension_over_reals()
2099 # Use the base ring of the matrix to ensure that its entries can be
2100 # multiplied by elements of the quaternion algebra.
2101 field
= M
.base_ring()
2102 Q
= QuaternionAlgebra(field
,-1,-1)
2105 # Go top-left to bottom-right (reading order), converting every
2106 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
2109 for l
in range(n
/d
):
2110 for m
in range(n
/d
):
2111 submat
= ComplexMatrixEuclideanJordanAlgebra
.real_unembed(
2112 M
[d
*l
:d
*l
+d
,d
*m
:d
*m
+d
] )
2113 if submat
[0,0] != submat
[1,1].conjugate():
2114 raise ValueError('bad on-diagonal submatrix')
2115 if submat
[0,1] != -submat
[1,0].conjugate():
2116 raise ValueError('bad off-diagonal submatrix')
2117 z
= submat
[0,0].real()
2118 z
+= submat
[0,0].imag()*i
2119 z
+= submat
[0,1].real()*j
2120 z
+= submat
[0,1].imag()*k
2123 return matrix(Q
, n
/d
, elements
)
2126 class QuaternionHermitianEJA(ConcreteEuclideanJordanAlgebra
,
2127 QuaternionMatrixEuclideanJordanAlgebra
):
2129 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
2130 matrices, the usual symmetric Jordan product, and the
2131 real-part-of-trace inner product. It has dimension `2n^2 - n` over
2136 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2140 In theory, our "field" can be any subfield of the reals::
2142 sage: QuaternionHermitianEJA(2, field=RDF)
2143 Euclidean Jordan algebra of dimension 6 over Real Double Field
2144 sage: QuaternionHermitianEJA(2, field=RR)
2145 Euclidean Jordan algebra of dimension 6 over Real Field with
2146 53 bits of precision
2150 The dimension of this algebra is `2*n^2 - n`::
2152 sage: set_random_seed()
2153 sage: n_max = QuaternionHermitianEJA._max_random_instance_size()
2154 sage: n = ZZ.random_element(1, n_max)
2155 sage: J = QuaternionHermitianEJA(n)
2156 sage: J.dimension() == 2*(n^2) - n
2159 The Jordan multiplication is what we think it is::
2161 sage: set_random_seed()
2162 sage: J = QuaternionHermitianEJA.random_instance()
2163 sage: x,y = J.random_elements(2)
2164 sage: actual = (x*y).to_matrix()
2165 sage: X = x.to_matrix()
2166 sage: Y = y.to_matrix()
2167 sage: expected = (X*Y + Y*X)/2
2168 sage: actual == expected
2170 sage: J(expected) == x*y
2173 We can change the generator prefix::
2175 sage: QuaternionHermitianEJA(2, prefix='a').gens()
2176 (a0, a1, a2, a3, a4, a5)
2178 We can construct the (trivial) algebra of rank zero::
2180 sage: QuaternionHermitianEJA(0)
2181 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2185 def _denormalized_basis(cls
, n
):
2187 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
2189 Why do we embed these? Basically, because all of numerical
2190 linear algebra assumes that you're working with vectors consisting
2191 of `n` entries from a field and scalars from the same field. There's
2192 no way to tell SageMath that (for example) the vectors contain
2193 complex numbers, while the scalar field is real.
2197 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2201 sage: set_random_seed()
2202 sage: n = ZZ.random_element(1,5)
2203 sage: B = QuaternionHermitianEJA._denormalized_basis(n)
2204 sage: all( M.is_symmetric() for M in B )
2209 Q
= QuaternionAlgebra(QQ
,-1,-1)
2212 # This is like the symmetric case, but we need to be careful:
2214 # * We want conjugate-symmetry, not just symmetry.
2215 # * The diagonal will (as a result) be real.
2219 for j
in range(i
+1):
2220 Eij
= matrix(Q
, n
, lambda k
,l
: k
==i
and l
==j
)
2222 Sij
= cls
.real_embed(Eij
)
2225 # The second, third, and fourth ones have a minus
2226 # because they're conjugated.
2227 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
2229 Sij_I
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
2231 Sij_J
= cls
.real_embed(J
*Eij
- J
*Eij
.transpose())
2233 Sij_K
= cls
.real_embed(K
*Eij
- K
*Eij
.transpose())
2236 # Since we embedded these, we can drop back to the "field" that we
2237 # started with instead of the quaternion algebra "Q".
2238 return tuple( s
.change_ring(field
) for s
in S
)
2241 def __init__(self
, n
, **kwargs
):
2242 # We know this is a valid EJA, but will double-check
2243 # if the user passes check_axioms=True.
2244 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2246 super(QuaternionHermitianEJA
, self
).__init
__(self
._denormalized
_basis
(n
),
2247 self
.jordan_product
,
2248 self
.trace_inner_product
,
2250 # TODO: this could be factored out somehow, but is left here
2251 # because the MatrixEuclideanJordanAlgebra is not presently
2252 # a subclass of the FDEJA class that defines rank() and one().
2253 self
.rank
.set_cache(n
)
2254 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
2255 self
.one
.set_cache(self(idV
))
2259 def _max_random_instance_size():
2261 The maximum rank of a random QuaternionHermitianEJA.
2263 return 2 # Dimension 6
2266 def random_instance(cls
, **kwargs
):
2268 Return a random instance of this type of algebra.
2270 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2271 return cls(n
, **kwargs
)
2274 class HadamardEJA(ConcreteEuclideanJordanAlgebra
):
2276 Return the Euclidean Jordan Algebra corresponding to the set
2277 `R^n` under the Hadamard product.
2279 Note: this is nothing more than the Cartesian product of ``n``
2280 copies of the spin algebra. Once Cartesian product algebras
2281 are implemented, this can go.
2285 sage: from mjo.eja.eja_algebra import HadamardEJA
2289 This multiplication table can be verified by hand::
2291 sage: J = HadamardEJA(3)
2292 sage: e0,e1,e2 = J.gens()
2308 We can change the generator prefix::
2310 sage: HadamardEJA(3, prefix='r').gens()
2314 def __init__(self
, n
, **kwargs
):
2315 def jordan_product(x
,y
):
2317 return P(tuple( xi
*yi
for (xi
,yi
) in zip(x
,y
) ))
2318 def inner_product(x
,y
):
2319 return x
.inner_product(y
)
2321 # New defaults for keyword arguments. Don't orthonormalize
2322 # because our basis is already orthonormal with respect to our
2323 # inner-product. Don't check the axioms, because we know this
2324 # is a valid EJA... but do double-check if the user passes
2325 # check_axioms=True. Note: we DON'T override the "check_field"
2326 # default here, because the user can pass in a field!
2327 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2328 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2331 standard_basis
= FreeModule(ZZ
, n
).basis()
2332 super(HadamardEJA
, self
).__init
__(standard_basis
,
2336 self
.rank
.set_cache(n
)
2339 self
.one
.set_cache( self
.zero() )
2341 self
.one
.set_cache( sum(self
.gens()) )
2344 def _max_random_instance_size():
2346 The maximum dimension of a random HadamardEJA.
2351 def random_instance(cls
, **kwargs
):
2353 Return a random instance of this type of algebra.
2355 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2356 return cls(n
, **kwargs
)
2359 class BilinearFormEJA(ConcreteEuclideanJordanAlgebra
):
2361 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2362 with the half-trace inner product and jordan product ``x*y =
2363 (<Bx,y>,y_bar>, x0*y_bar + y0*x_bar)`` where `B = 1 \times B22` is
2364 a symmetric positive-definite "bilinear form" matrix. Its
2365 dimension is the size of `B`, and it has rank two in dimensions
2366 larger than two. It reduces to the ``JordanSpinEJA`` when `B` is
2367 the identity matrix of order ``n``.
2369 We insist that the one-by-one upper-left identity block of `B` be
2370 passed in as well so that we can be passed a matrix of size zero
2371 to construct a trivial algebra.
2375 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
2376 ....: JordanSpinEJA)
2380 When no bilinear form is specified, the identity matrix is used,
2381 and the resulting algebra is the Jordan spin algebra::
2383 sage: B = matrix.identity(AA,3)
2384 sage: J0 = BilinearFormEJA(B)
2385 sage: J1 = JordanSpinEJA(3)
2386 sage: J0.multiplication_table() == J0.multiplication_table()
2389 An error is raised if the matrix `B` does not correspond to a
2390 positive-definite bilinear form::
2392 sage: B = matrix.random(QQ,2,3)
2393 sage: J = BilinearFormEJA(B)
2394 Traceback (most recent call last):
2396 ValueError: bilinear form is not positive-definite
2397 sage: B = matrix.zero(QQ,3)
2398 sage: J = BilinearFormEJA(B)
2399 Traceback (most recent call last):
2401 ValueError: bilinear form is not positive-definite
2405 We can create a zero-dimensional algebra::
2407 sage: B = matrix.identity(AA,0)
2408 sage: J = BilinearFormEJA(B)
2412 We can check the multiplication condition given in the Jordan, von
2413 Neumann, and Wigner paper (and also discussed on my "On the
2414 symmetry..." paper). Note that this relies heavily on the standard
2415 choice of basis, as does anything utilizing the bilinear form
2416 matrix. We opt not to orthonormalize the basis, because if we
2417 did, we would have to normalize the `s_{i}` in a similar manner::
2419 sage: set_random_seed()
2420 sage: n = ZZ.random_element(5)
2421 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2422 sage: B11 = matrix.identity(QQ,1)
2423 sage: B22 = M.transpose()*M
2424 sage: B = block_matrix(2,2,[ [B11,0 ],
2426 sage: J = BilinearFormEJA(B, orthonormalize=False)
2427 sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
2428 sage: V = J.vector_space()
2429 sage: sis = [ J( V([0] + (M.inverse()*ei).list()).column() )
2430 ....: for ei in eis ]
2431 sage: actual = [ sis[i]*sis[j]
2432 ....: for i in range(n-1)
2433 ....: for j in range(n-1) ]
2434 sage: expected = [ J.one() if i == j else J.zero()
2435 ....: for i in range(n-1)
2436 ....: for j in range(n-1) ]
2437 sage: actual == expected
2440 def __init__(self
, B
, **kwargs
):
2441 if not B
.is_positive_definite():
2442 raise ValueError("bilinear form is not positive-definite")
2444 def inner_product(x
,y
):
2445 return (B
*x
).inner_product(y
)
2447 def jordan_product(x
,y
):
2453 z0
= inner_product(x
,y
)
2454 zbar
= y0
*xbar
+ x0
*ybar
2455 return P((z0
,) + tuple(zbar
))
2457 # We know this is a valid EJA, but will double-check
2458 # if the user passes check_axioms=True.
2459 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2462 standard_basis
= FreeModule(ZZ
, n
).basis()
2463 super(BilinearFormEJA
, self
).__init
__(standard_basis
,
2468 # The rank of this algebra is two, unless we're in a
2469 # one-dimensional ambient space (because the rank is bounded
2470 # by the ambient dimension).
2471 self
.rank
.set_cache(min(n
,2))
2474 self
.one
.set_cache( self
.zero() )
2476 self
.one
.set_cache( self
.monomial(0) )
2479 def _max_random_instance_size():
2481 The maximum dimension of a random BilinearFormEJA.
2486 def random_instance(cls
, **kwargs
):
2488 Return a random instance of this algebra.
2490 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2492 B
= matrix
.identity(ZZ
, n
)
2493 return cls(B
, **kwargs
)
2495 B11
= matrix
.identity(ZZ
, 1)
2496 M
= matrix
.random(ZZ
, n
-1)
2497 I
= matrix
.identity(ZZ
, n
-1)
2499 while alpha
.is_zero():
2500 alpha
= ZZ
.random_element().abs()
2501 B22
= M
.transpose()*M
+ alpha
*I
2503 from sage
.matrix
.special
import block_matrix
2504 B
= block_matrix(2,2, [ [B11
, ZZ(0) ],
2507 return cls(B
, **kwargs
)
2510 class JordanSpinEJA(BilinearFormEJA
):
2512 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2513 with the usual inner product and jordan product ``x*y =
2514 (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2519 sage: from mjo.eja.eja_algebra import JordanSpinEJA
2523 This multiplication table can be verified by hand::
2525 sage: J = JordanSpinEJA(4)
2526 sage: e0,e1,e2,e3 = J.gens()
2542 We can change the generator prefix::
2544 sage: JordanSpinEJA(2, prefix='B').gens()
2549 Ensure that we have the usual inner product on `R^n`::
2551 sage: set_random_seed()
2552 sage: J = JordanSpinEJA.random_instance()
2553 sage: x,y = J.random_elements(2)
2554 sage: actual = x.inner_product(y)
2555 sage: expected = x.to_vector().inner_product(y.to_vector())
2556 sage: actual == expected
2560 def __init__(self
, n
, **kwargs
):
2561 # This is a special case of the BilinearFormEJA with the
2562 # identity matrix as its bilinear form.
2563 B
= matrix
.identity(ZZ
, n
)
2565 # Don't orthonormalize because our basis is already
2566 # orthonormal with respect to our inner-product.
2567 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2569 # But also don't pass check_field=False here, because the user
2570 # can pass in a field!
2571 super(JordanSpinEJA
, self
).__init
__(B
, **kwargs
)
2574 def _max_random_instance_size():
2576 The maximum dimension of a random JordanSpinEJA.
2581 def random_instance(cls
, **kwargs
):
2583 Return a random instance of this type of algebra.
2585 Needed here to override the implementation for ``BilinearFormEJA``.
2587 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2588 return cls(n
, **kwargs
)
2591 class TrivialEJA(ConcreteEuclideanJordanAlgebra
):
2593 The trivial Euclidean Jordan algebra consisting of only a zero element.
2597 sage: from mjo.eja.eja_algebra import TrivialEJA
2601 sage: J = TrivialEJA()
2608 sage: 7*J.one()*12*J.one()
2610 sage: J.one().inner_product(J.one())
2612 sage: J.one().norm()
2614 sage: J.one().subalgebra_generated_by()
2615 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2620 def __init__(self
, **kwargs
):
2621 jordan_product
= lambda x
,y
: x
2622 inner_product
= lambda x
,y
: 0
2625 # New defaults for keyword arguments
2626 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2627 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2629 super(TrivialEJA
, self
).__init
__(basis
,
2633 # The rank is zero using my definition, namely the dimension of the
2634 # largest subalgebra generated by any element.
2635 self
.rank
.set_cache(0)
2636 self
.one
.set_cache( self
.zero() )
2639 def random_instance(cls
, **kwargs
):
2640 # We don't take a "size" argument so the superclass method is
2641 # inappropriate for us.
2642 return cls(**kwargs
)
2644 class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra
):
2646 The external (orthogonal) direct sum of two other Euclidean Jordan
2647 algebras. Essentially the Cartesian product of its two factors.
2648 Every Euclidean Jordan algebra decomposes into an orthogonal
2649 direct sum of simple Euclidean Jordan algebras, so no generality
2650 is lost by providing only this construction.
2654 sage: from mjo.eja.eja_algebra import (random_eja,
2656 ....: RealSymmetricEJA,
2661 sage: J1 = HadamardEJA(2)
2662 sage: J2 = RealSymmetricEJA(3)
2663 sage: J = DirectSumEJA(J1,J2)
2671 The external direct sum construction is only valid when the two factors
2672 have the same base ring; an error is raised otherwise::
2674 sage: set_random_seed()
2675 sage: J1 = random_eja(field=AA)
2676 sage: J2 = random_eja(field=QQ,orthonormalize=False)
2677 sage: J = DirectSumEJA(J1,J2)
2678 Traceback (most recent call last):
2680 ValueError: algebras must share the same base field
2683 def __init__(self
, J1
, J2
, **kwargs
):
2684 if J1
.base_ring() != J2
.base_ring():
2685 raise ValueError("algebras must share the same base field")
2686 field
= J1
.base_ring()
2688 self
._factors
= (J1
, J2
)
2692 V
= VectorSpace(field
, n
)
2693 mult_table
= [ [ V
.zero() for j
in range(i
+1) ]
2696 for j
in range(i
+1):
2697 p
= (J1
.monomial(i
)*J1
.monomial(j
)).to_vector()
2698 mult_table
[i
][j
] = V(p
.list() + [field
.zero()]*n2
)
2701 for j
in range(i
+1):
2702 p
= (J2
.monomial(i
)*J2
.monomial(j
)).to_vector()
2703 mult_table
[n1
+i
][n1
+j
] = V([field
.zero()]*n1
+ p
.list())
2705 # TODO: build the IP table here from the two constituent IP
2706 # matrices (it'll be block diagonal, I think).
2707 ip_table
= [ [ field
.zero() for j
in range(i
+1) ]
2709 super(DirectSumEJA
, self
).__init
__(field
,
2714 self
.rank
.set_cache(J1
.rank() + J2
.rank())
2719 Return the pair of this algebra's factors.
2723 sage: from mjo.eja.eja_algebra import (HadamardEJA,
2724 ....: JordanSpinEJA,
2729 sage: J1 = HadamardEJA(2, field=QQ)
2730 sage: J2 = JordanSpinEJA(3, field=QQ)
2731 sage: J = DirectSumEJA(J1,J2)
2733 (Euclidean Jordan algebra of dimension 2 over Rational Field,
2734 Euclidean Jordan algebra of dimension 3 over Rational Field)
2737 return self
._factors
2739 def projections(self
):
2741 Return a pair of projections onto this algebra's factors.
2745 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
2746 ....: ComplexHermitianEJA,
2751 sage: J1 = JordanSpinEJA(2)
2752 sage: J2 = ComplexHermitianEJA(2)
2753 sage: J = DirectSumEJA(J1,J2)
2754 sage: (pi_left, pi_right) = J.projections()
2755 sage: J.one().to_vector()
2757 sage: pi_left(J.one()).to_vector()
2759 sage: pi_right(J.one()).to_vector()
2763 (J1
,J2
) = self
.factors()
2766 V_basis
= self
.vector_space().basis()
2767 # Need to specify the dimensions explicitly so that we don't
2768 # wind up with a zero-by-zero matrix when we want e.g. a
2769 # zero-by-two matrix (important for composing things).
2770 P1
= matrix(self
.base_ring(), m
, m
+n
, V_basis
[:m
])
2771 P2
= matrix(self
.base_ring(), n
, m
+n
, V_basis
[m
:])
2772 pi_left
= FiniteDimensionalEuclideanJordanAlgebraOperator(self
,J1
,P1
)
2773 pi_right
= FiniteDimensionalEuclideanJordanAlgebraOperator(self
,J2
,P2
)
2774 return (pi_left
, pi_right
)
2776 def inclusions(self
):
2778 Return the pair of inclusion maps from our factors into us.
2782 sage: from mjo.eja.eja_algebra import (random_eja,
2783 ....: JordanSpinEJA,
2784 ....: RealSymmetricEJA,
2789 sage: J1 = JordanSpinEJA(3)
2790 sage: J2 = RealSymmetricEJA(2)
2791 sage: J = DirectSumEJA(J1,J2)
2792 sage: (iota_left, iota_right) = J.inclusions()
2793 sage: iota_left(J1.zero()) == J.zero()
2795 sage: iota_right(J2.zero()) == J.zero()
2797 sage: J1.one().to_vector()
2799 sage: iota_left(J1.one()).to_vector()
2801 sage: J2.one().to_vector()
2803 sage: iota_right(J2.one()).to_vector()
2805 sage: J.one().to_vector()
2810 Composing a projection with the corresponding inclusion should
2811 produce the identity map, and mismatching them should produce
2814 sage: set_random_seed()
2815 sage: J1 = random_eja()
2816 sage: J2 = random_eja()
2817 sage: J = DirectSumEJA(J1,J2)
2818 sage: (iota_left, iota_right) = J.inclusions()
2819 sage: (pi_left, pi_right) = J.projections()
2820 sage: pi_left*iota_left == J1.one().operator()
2822 sage: pi_right*iota_right == J2.one().operator()
2824 sage: (pi_left*iota_right).is_zero()
2826 sage: (pi_right*iota_left).is_zero()
2830 (J1
,J2
) = self
.factors()
2833 V_basis
= self
.vector_space().basis()
2834 # Need to specify the dimensions explicitly so that we don't
2835 # wind up with a zero-by-zero matrix when we want e.g. a
2836 # two-by-zero matrix (important for composing things).
2837 I1
= matrix
.column(self
.base_ring(), m
, m
+n
, V_basis
[:m
])
2838 I2
= matrix
.column(self
.base_ring(), n
, m
+n
, V_basis
[m
:])
2839 iota_left
= FiniteDimensionalEuclideanJordanAlgebraOperator(J1
,self
,I1
)
2840 iota_right
= FiniteDimensionalEuclideanJordanAlgebraOperator(J2
,self
,I2
)
2841 return (iota_left
, iota_right
)
2843 def inner_product(self
, x
, y
):
2845 The standard Cartesian inner-product.
2847 We project ``x`` and ``y`` onto our factors, and add up the
2848 inner-products from the subalgebras.
2853 sage: from mjo.eja.eja_algebra import (HadamardEJA,
2854 ....: QuaternionHermitianEJA,
2859 sage: J1 = HadamardEJA(3,field=QQ)
2860 sage: J2 = QuaternionHermitianEJA(2,field=QQ,orthonormalize=False)
2861 sage: J = DirectSumEJA(J1,J2)
2866 sage: x1.inner_product(x2)
2868 sage: y1.inner_product(y2)
2870 sage: J.one().inner_product(J.one())
2874 (pi_left
, pi_right
) = self
.projections()
2880 return (x1
.inner_product(y1
) + x2
.inner_product(y2
))
2884 random_eja
= ConcreteEuclideanJordanAlgebra
.random_instance