2 Euclidean Jordan Algebras. These are formally-real Jordan Algebras;
3 specifically those where u^2 + v^2 = 0 implies that u = v = 0. They
4 are used in optimization, and have some additional nice methods beyond
5 what can be supported in a general Jordan Algebra.
10 sage: from mjo.eja.eja_algebra import random_eja
15 Euclidean Jordan algebra of dimension...
19 from itertools
import repeat
21 from sage
.algebras
.quatalg
.quaternion_algebra
import QuaternionAlgebra
22 from sage
.categories
.magmatic_algebras
import MagmaticAlgebras
23 from sage
.combinat
.free_module
import CombinatorialFreeModule
24 from sage
.matrix
.constructor
import matrix
25 from sage
.matrix
.matrix_space
import MatrixSpace
26 from sage
.misc
.cachefunc
import cached_method
27 from sage
.misc
.table
import table
28 from sage
.modules
.free_module
import FreeModule
, VectorSpace
29 from sage
.rings
.all
import (ZZ
, QQ
, AA
, QQbar
, RR
, RLF
, CLF
,
32 from mjo
.eja
.eja_element
import FiniteDimensionalEuclideanJordanAlgebraElement
33 from mjo
.eja
.eja_operator
import FiniteDimensionalEuclideanJordanAlgebraOperator
34 from mjo
.eja
.eja_utils
import _mat2vec
36 class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule
):
38 The lowest-level class for representing a Euclidean Jordan algebra.
40 def _coerce_map_from_base_ring(self
):
42 Disable the map from the base ring into the algebra.
44 Performing a nonsense conversion like this automatically
45 is counterpedagogical. The fallback is to try the usual
46 element constructor, which should also fail.
50 sage: from mjo.eja.eja_algebra import random_eja
54 sage: set_random_seed()
55 sage: J = random_eja()
57 Traceback (most recent call last):
59 ValueError: not an element of this algebra
76 * field -- the scalar field for this algebra (must be real)
78 * multiplication_table -- the multiplication table for this
79 algebra's implicit basis. Only the lower-triangular portion
80 of the table is used, since the multiplication is assumed
85 sage: from mjo.eja.eja_algebra import (
86 ....: FiniteDimensionalEuclideanJordanAlgebra,
92 By definition, Jordan multiplication commutes::
94 sage: set_random_seed()
95 sage: J = random_eja()
96 sage: x,y = J.random_elements(2)
100 An error is raised if the Jordan product is not commutative::
102 sage: JP = ((1,2),(0,0))
103 sage: IP = ((1,0),(0,1))
104 sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,JP,IP)
105 Traceback (most recent call last):
107 ValueError: Jordan product is not commutative
109 An error is raised if the inner-product is not commutative::
111 sage: JP = ((1,0),(0,1))
112 sage: IP = ((1,2),(0,0))
113 sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,JP,IP)
114 Traceback (most recent call last):
116 ValueError: inner-product is not commutative
120 The ``field`` we're given must be real with ``check_field=True``::
122 sage: JordanSpinEJA(2,QQbar)
123 Traceback (most recent call last):
125 ValueError: scalar field is not real
127 The multiplication table must be square with ``check_axioms=True``::
129 sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,((),()),((1,),))
130 Traceback (most recent call last):
132 ValueError: multiplication table is not square
134 The multiplication and inner-product tables must be the same
135 size (and in particular, the inner-product table must also be
136 square) with ``check_axioms=True``::
138 sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,((1,),),(()))
139 Traceback (most recent call last):
141 ValueError: multiplication and inner-product tables are
143 sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,((1,),),((1,2),))
144 Traceback (most recent call last):
146 ValueError: multiplication and inner-product tables are
151 if not field
.is_subring(RR
):
152 # Note: this does return true for the real algebraic
153 # field, the rationals, and any quadratic field where
154 # we've specified a real embedding.
155 raise ValueError("scalar field is not real")
158 # The multiplication and inner-product tables should be square
159 # if the user wants us to verify them. And we verify them as
160 # soon as possible, because we want to exploit their symmetry.
161 n
= len(multiplication_table
)
163 if not all( len(l
) == n
for l
in multiplication_table
):
164 raise ValueError("multiplication table is not square")
166 # If the multiplication table is square, we can check if
167 # the inner-product table is square by comparing it to the
168 # multiplication table's dimensions.
169 msg
= "multiplication and inner-product tables are different sizes"
170 if not len(inner_product_table
) == n
:
171 raise ValueError(msg
)
173 if not all( len(l
) == n
for l
in inner_product_table
):
174 raise ValueError(msg
)
176 # Check commutativity of the Jordan product (symmetry of
177 # the multiplication table) and the commutativity of the
178 # inner-product (symmetry of the inner-product table)
179 # first if we're going to check them at all.. This has to
180 # be done before we define product_on_basis(), because
181 # that method assumes that self._multiplication_table is
182 # symmetric. And it has to be done before we build
183 # self._inner_product_matrix, because the process used to
184 # construct it assumes symmetry as well.
185 if not all( multiplication_table
[j
][i
]
186 == multiplication_table
[i
][j
]
188 for j
in range(i
+1) ):
189 raise ValueError("Jordan product is not commutative")
191 if not all( inner_product_table
[j
][i
]
192 == inner_product_table
[i
][j
]
194 for j
in range(i
+1) ):
195 raise ValueError("inner-product is not commutative")
197 self
._matrix
_basis
= matrix_basis
200 category
= MagmaticAlgebras(field
).FiniteDimensional()
201 category
= category
.WithBasis().Unital()
203 fda
= super(FiniteDimensionalEuclideanJordanAlgebra
, self
)
208 self
.print_options(bracket
='')
210 # The multiplication table we're given is necessarily in terms
211 # of vectors, because we don't have an algebra yet for
212 # anything to be an element of. However, it's faster in the
213 # long run to have the multiplication table be in terms of
214 # algebra elements. We do this after calling the superclass
215 # constructor so that from_vector() knows what to do.
217 # Note: we take advantage of symmetry here, and only store
218 # the lower-triangular portion of the table.
219 self
._multiplication
_table
= [ [ self
.vector_space().zero()
220 for j
in range(i
+1) ]
225 elt
= self
.from_vector(multiplication_table
[i
][j
])
226 self
._multiplication
_table
[i
][j
] = elt
228 self
._multiplication
_table
= tuple(map(tuple, self
._multiplication
_table
))
230 # Save our inner product as a matrix, since the efficiency of
231 # matrix multiplication will usually outweigh the fact that we
232 # have to store a redundant upper- or lower-triangular part.
233 # Pre-cache the fact that these are Hermitian (real symmetric,
234 # in fact) in case some e.g. matrix multiplication routine can
235 # take advantage of it.
236 ip_matrix_constructor
= lambda i
,j
: inner_product_table
[i
][j
] if j
<= i
else inner_product_table
[j
][i
]
237 self
._inner
_product
_matrix
= matrix(field
, n
, ip_matrix_constructor
)
238 self
._inner
_product
_matrix
._cache
= {'hermitian': True}
239 self
._inner
_product
_matrix
.set_immutable()
242 if not self
._is
_jordanian
():
243 raise ValueError("Jordan identity does not hold")
244 if not self
._inner
_product
_is
_associative
():
245 raise ValueError("inner product is not associative")
247 def _element_constructor_(self
, elt
):
249 Construct an element of this algebra from its vector or matrix
252 This gets called only after the parent element _call_ method
253 fails to find a coercion for the argument.
257 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
259 ....: RealSymmetricEJA)
263 The identity in `S^n` is converted to the identity in the EJA::
265 sage: J = RealSymmetricEJA(3)
266 sage: I = matrix.identity(QQ,3)
267 sage: J(I) == J.one()
270 This skew-symmetric matrix can't be represented in the EJA::
272 sage: J = RealSymmetricEJA(3)
273 sage: A = matrix(QQ,3, lambda i,j: i-j)
275 Traceback (most recent call last):
277 ValueError: not an element of this algebra
281 Ensure that we can convert any element of the two non-matrix
282 simple algebras (whose matrix representations are columns)
283 back and forth faithfully::
285 sage: set_random_seed()
286 sage: J = HadamardEJA.random_instance()
287 sage: x = J.random_element()
288 sage: J(x.to_vector().column()) == x
290 sage: J = JordanSpinEJA.random_instance()
291 sage: x = J.random_element()
292 sage: J(x.to_vector().column()) == x
295 msg
= "not an element of this algebra"
297 # The superclass implementation of random_element()
298 # needs to be able to coerce "0" into the algebra.
300 elif elt
in self
.base_ring():
301 # Ensure that no base ring -> algebra coercion is performed
302 # by this method. There's some stupidity in sage that would
303 # otherwise propagate to this method; for example, sage thinks
304 # that the integer 3 belongs to the space of 2-by-2 matrices.
305 raise ValueError(msg
)
307 if elt
not in self
.matrix_space():
308 raise ValueError(msg
)
310 # Thanks for nothing! Matrix spaces aren't vector spaces in
311 # Sage, so we have to figure out its matrix-basis coordinates
312 # ourselves. We use the basis space's ring instead of the
313 # element's ring because the basis space might be an algebraic
314 # closure whereas the base ring of the 3-by-3 identity matrix
315 # could be QQ instead of QQbar.
316 V
= VectorSpace(self
.base_ring(), elt
.nrows()*elt
.ncols())
317 W
= V
.span_of_basis( _mat2vec(s
) for s
in self
.matrix_basis() )
320 coords
= W
.coordinate_vector(_mat2vec(elt
))
321 except ArithmeticError: # vector is not in free module
322 raise ValueError(msg
)
324 return self
.from_vector(coords
)
328 Return a string representation of ``self``.
332 sage: from mjo.eja.eja_algebra import JordanSpinEJA
336 Ensure that it says what we think it says::
338 sage: JordanSpinEJA(2, field=AA)
339 Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
340 sage: JordanSpinEJA(3, field=RDF)
341 Euclidean Jordan algebra of dimension 3 over Real Double Field
344 fmt
= "Euclidean Jordan algebra of dimension {} over {}"
345 return fmt
.format(self
.dimension(), self
.base_ring())
347 def product_on_basis(self
, i
, j
):
348 # We only stored the lower-triangular portion of the
349 # multiplication table.
351 return self
._multiplication
_table
[i
][j
]
353 return self
._multiplication
_table
[j
][i
]
355 def _is_commutative(self
):
357 Whether or not this algebra's multiplication table is commutative.
359 This method should of course always return ``True``, unless
360 this algebra was constructed with ``check_axioms=False`` and
361 passed an invalid multiplication table.
363 return all( self
.product_on_basis(i
,j
) == self
.product_on_basis(i
,j
)
364 for i
in range(self
.dimension())
365 for j
in range(self
.dimension()) )
367 def _is_jordanian(self
):
369 Whether or not this algebra's multiplication table respects the
370 Jordan identity `(x^{2})(xy) = x(x^{2}y)`.
372 We only check one arrangement of `x` and `y`, so for a
373 ``True`` result to be truly true, you should also check
374 :meth:`_is_commutative`. This method should of course always
375 return ``True``, unless this algebra was constructed with
376 ``check_axioms=False`` and passed an invalid multiplication table.
378 return all( (self
.monomial(i
)**2)*(self
.monomial(i
)*self
.monomial(j
))
380 (self
.monomial(i
))*((self
.monomial(i
)**2)*self
.monomial(j
))
381 for i
in range(self
.dimension())
382 for j
in range(self
.dimension()) )
384 def _inner_product_is_associative(self
):
386 Return whether or not this algebra's inner product `B` is
387 associative; that is, whether or not `B(xy,z) = B(x,yz)`.
389 This method should of course always return ``True``, unless
390 this algebra was constructed with ``check_axioms=False`` and
391 passed an invalid multiplication table.
394 # Used to check whether or not something is zero in an inexact
395 # ring. This number is sufficient to allow the construction of
396 # QuaternionHermitianEJA(2, RDF) with check_axioms=True.
399 for i
in range(self
.dimension()):
400 for j
in range(self
.dimension()):
401 for k
in range(self
.dimension()):
405 diff
= (x
*y
).inner_product(z
) - x
.inner_product(y
*z
)
407 if self
.base_ring().is_exact():
411 if diff
.abs() > epsilon
:
417 def characteristic_polynomial_of(self
):
419 Return the algebra's "characteristic polynomial of" function,
420 which is itself a multivariate polynomial that, when evaluated
421 at the coordinates of some algebra element, returns that
422 element's characteristic polynomial.
424 The resulting polynomial has `n+1` variables, where `n` is the
425 dimension of this algebra. The first `n` variables correspond to
426 the coordinates of an algebra element: when evaluated at the
427 coordinates of an algebra element with respect to a certain
428 basis, the result is a univariate polynomial (in the one
429 remaining variable ``t``), namely the characteristic polynomial
434 sage: from mjo.eja.eja_algebra import JordanSpinEJA, TrivialEJA
438 The characteristic polynomial in the spin algebra is given in
439 Alizadeh, Example 11.11::
441 sage: J = JordanSpinEJA(3)
442 sage: p = J.characteristic_polynomial_of(); p
443 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
444 sage: xvec = J.one().to_vector()
448 By definition, the characteristic polynomial is a monic
449 degree-zero polynomial in a rank-zero algebra. Note that
450 Cayley-Hamilton is indeed satisfied since the polynomial
451 ``1`` evaluates to the identity element of the algebra on
454 sage: J = TrivialEJA()
455 sage: J.characteristic_polynomial_of()
462 # The list of coefficient polynomials a_0, a_1, a_2, ..., a_(r-1).
463 a
= self
._charpoly
_coefficients
()
465 # We go to a bit of trouble here to reorder the
466 # indeterminates, so that it's easier to evaluate the
467 # characteristic polynomial at x's coordinates and get back
468 # something in terms of t, which is what we want.
469 S
= PolynomialRing(self
.base_ring(),'t')
473 S
= PolynomialRing(S
, R
.variable_names())
476 return (t
**r
+ sum( a
[k
]*(t
**k
) for k
in range(r
) ))
478 def coordinate_polynomial_ring(self
):
480 The multivariate polynomial ring in which this algebra's
481 :meth:`characteristic_polynomial_of` lives.
485 sage: from mjo.eja.eja_algebra import (HadamardEJA,
486 ....: RealSymmetricEJA)
490 sage: J = HadamardEJA(2)
491 sage: J.coordinate_polynomial_ring()
492 Multivariate Polynomial Ring in X1, X2...
493 sage: J = RealSymmetricEJA(3,QQ,orthonormalize=False)
494 sage: J.coordinate_polynomial_ring()
495 Multivariate Polynomial Ring in X1, X2, X3, X4, X5, X6...
498 var_names
= tuple( "X%d" % z
for z
in range(1, self
.dimension()+1) )
499 return PolynomialRing(self
.base_ring(), var_names
)
501 def inner_product(self
, x
, y
):
503 The inner product associated with this Euclidean Jordan algebra.
505 Defaults to the trace inner product, but can be overridden by
506 subclasses if they are sure that the necessary properties are
511 sage: from mjo.eja.eja_algebra import (random_eja,
513 ....: BilinearFormEJA)
517 Our inner product is "associative," which means the following for
518 a symmetric bilinear form::
520 sage: set_random_seed()
521 sage: J = random_eja()
522 sage: x,y,z = J.random_elements(3)
523 sage: (x*y).inner_product(z) == y.inner_product(x*z)
528 Ensure that this is the usual inner product for the algebras
531 sage: set_random_seed()
532 sage: J = HadamardEJA.random_instance()
533 sage: x,y = J.random_elements(2)
534 sage: actual = x.inner_product(y)
535 sage: expected = x.to_vector().inner_product(y.to_vector())
536 sage: actual == expected
539 Ensure that this is one-half of the trace inner-product in a
540 BilinearFormEJA that isn't just the reals (when ``n`` isn't
541 one). This is in Faraut and Koranyi, and also my "On the
544 sage: set_random_seed()
545 sage: J = BilinearFormEJA.random_instance()
546 sage: n = J.dimension()
547 sage: x = J.random_element()
548 sage: y = J.random_element()
549 sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
552 B
= self
._inner
_product
_matrix
553 return (B
*x
.to_vector()).inner_product(y
.to_vector())
556 def is_trivial(self
):
558 Return whether or not this algebra is trivial.
560 A trivial algebra contains only the zero element.
564 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
569 sage: J = ComplexHermitianEJA(3)
575 sage: J = TrivialEJA()
580 return self
.dimension() == 0
583 def multiplication_table(self
):
585 Return a visual representation of this algebra's multiplication
586 table (on basis elements).
590 sage: from mjo.eja.eja_algebra import JordanSpinEJA
594 sage: J = JordanSpinEJA(4)
595 sage: J.multiplication_table()
596 +----++----+----+----+----+
597 | * || e0 | e1 | e2 | e3 |
598 +====++====+====+====+====+
599 | e0 || e0 | e1 | e2 | e3 |
600 +----++----+----+----+----+
601 | e1 || e1 | e0 | 0 | 0 |
602 +----++----+----+----+----+
603 | e2 || e2 | 0 | e0 | 0 |
604 +----++----+----+----+----+
605 | e3 || e3 | 0 | 0 | e0 |
606 +----++----+----+----+----+
610 M
= [ [ self
.zero() for j
in range(n
) ]
614 M
[i
][j
] = self
._multiplication
_table
[i
][j
]
618 # Prepend the left "header" column entry Can't do this in
619 # the loop because it messes up the symmetry.
620 M
[i
] = [self
.monomial(i
)] + M
[i
]
622 # Prepend the header row.
623 M
= [["*"] + list(self
.gens())] + M
624 return table(M
, header_row
=True, header_column
=True, frame
=True)
627 def matrix_basis(self
):
629 Return an (often more natural) representation of this algebras
630 basis as an ordered tuple of matrices.
632 Every finite-dimensional Euclidean Jordan Algebra is a, up to
633 Jordan isomorphism, a direct sum of five simple
634 algebras---four of which comprise Hermitian matrices. And the
635 last type of algebra can of course be thought of as `n`-by-`1`
636 column matrices (ambiguusly called column vectors) to avoid
637 special cases. As a result, matrices (and column vectors) are
638 a natural representation format for Euclidean Jordan algebra
641 But, when we construct an algebra from a basis of matrices,
642 those matrix representations are lost in favor of coordinate
643 vectors *with respect to* that basis. We could eventually
644 convert back if we tried hard enough, but having the original
645 representations handy is valuable enough that we simply store
646 them and return them from this method.
648 Why implement this for non-matrix algebras? Avoiding special
649 cases for the :class:`BilinearFormEJA` pays with simplicity in
650 its own right. But mainly, we would like to be able to assume
651 that elements of a :class:`DirectSumEJA` can be displayed
652 nicely, without having to have special classes for direct sums
653 one of whose components was a matrix algebra.
657 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
658 ....: RealSymmetricEJA)
662 sage: J = RealSymmetricEJA(2)
664 Finite family {0: e0, 1: e1, 2: e2}
665 sage: J.matrix_basis()
667 [1 0] [ 0 0.7071067811865475?] [0 0]
668 [0 0], [0.7071067811865475? 0], [0 1]
673 sage: J = JordanSpinEJA(2)
675 Finite family {0: e0, 1: e1}
676 sage: J.matrix_basis()
682 if self
._matrix
_basis
is None:
683 M
= self
.matrix_space()
684 return tuple( M(b
.to_vector()) for b
in self
.basis() )
686 return self
._matrix
_basis
689 def matrix_space(self
):
691 Return the matrix space in which this algebra's elements live, if
692 we think of them as matrices (including column vectors of the
695 Generally this will be an `n`-by-`1` column-vector space,
696 except when the algebra is trivial. There it's `n`-by-`n`
697 (where `n` is zero), to ensure that two elements of the matrix
698 space (empty matrices) can be multiplied.
700 Matrix algebras override this with something more useful.
702 if self
.is_trivial():
703 return MatrixSpace(self
.base_ring(), 0)
704 elif self
._matrix
_basis
is None or len(self
._matrix
_basis
) == 0:
705 return MatrixSpace(self
.base_ring(), self
.dimension(), 1)
707 return self
._matrix
_basis
[0].matrix_space()
713 Return the unit element of this algebra.
717 sage: from mjo.eja.eja_algebra import (HadamardEJA,
722 sage: J = HadamardEJA(5)
724 e0 + e1 + e2 + e3 + e4
728 The identity element acts like the identity::
730 sage: set_random_seed()
731 sage: J = random_eja()
732 sage: x = J.random_element()
733 sage: J.one()*x == x and x*J.one() == x
736 The matrix of the unit element's operator is the identity::
738 sage: set_random_seed()
739 sage: J = random_eja()
740 sage: actual = J.one().operator().matrix()
741 sage: expected = matrix.identity(J.base_ring(), J.dimension())
742 sage: actual == expected
745 Ensure that the cached unit element (often precomputed by
746 hand) agrees with the computed one::
748 sage: set_random_seed()
749 sage: J = random_eja()
750 sage: cached = J.one()
751 sage: J.one.clear_cache()
752 sage: J.one() == cached
756 # We can brute-force compute the matrices of the operators
757 # that correspond to the basis elements of this algebra.
758 # If some linear combination of those basis elements is the
759 # algebra identity, then the same linear combination of
760 # their matrices has to be the identity matrix.
762 # Of course, matrices aren't vectors in sage, so we have to
763 # appeal to the "long vectors" isometry.
764 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
766 # Now we use basic linear algebra to find the coefficients,
767 # of the matrices-as-vectors-linear-combination, which should
768 # work for the original algebra basis too.
769 A
= matrix(self
.base_ring(), oper_vecs
)
771 # We used the isometry on the left-hand side already, but we
772 # still need to do it for the right-hand side. Recall that we
773 # wanted something that summed to the identity matrix.
774 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
776 # Now if there's an identity element in the algebra, this
777 # should work. We solve on the left to avoid having to
778 # transpose the matrix "A".
779 return self
.from_vector(A
.solve_left(b
))
782 def peirce_decomposition(self
, c
):
784 The Peirce decomposition of this algebra relative to the
787 In the future, this can be extended to a complete system of
788 orthogonal idempotents.
792 - ``c`` -- an idempotent of this algebra.
796 A triple (J0, J5, J1) containing two subalgebras and one subspace
799 - ``J0`` -- the algebra on the eigenspace of ``c.operator()``
800 corresponding to the eigenvalue zero.
802 - ``J5`` -- the eigenspace (NOT a subalgebra) of ``c.operator()``
803 corresponding to the eigenvalue one-half.
805 - ``J1`` -- the algebra on the eigenspace of ``c.operator()``
806 corresponding to the eigenvalue one.
808 These are the only possible eigenspaces for that operator, and this
809 algebra is a direct sum of them. The spaces ``J0`` and ``J1`` are
810 orthogonal, and are subalgebras of this algebra with the appropriate
815 sage: from mjo.eja.eja_algebra import random_eja, RealSymmetricEJA
819 The canonical example comes from the symmetric matrices, which
820 decompose into diagonal and off-diagonal parts::
822 sage: J = RealSymmetricEJA(3)
823 sage: C = matrix(QQ, [ [1,0,0],
827 sage: J0,J5,J1 = J.peirce_decomposition(c)
829 Euclidean Jordan algebra of dimension 1...
831 Vector space of degree 6 and dimension 2...
833 Euclidean Jordan algebra of dimension 3...
834 sage: J0.one().to_matrix()
838 sage: orig_df = AA.options.display_format
839 sage: AA.options.display_format = 'radical'
840 sage: J.from_vector(J5.basis()[0]).to_matrix()
844 sage: J.from_vector(J5.basis()[1]).to_matrix()
848 sage: AA.options.display_format = orig_df
849 sage: J1.one().to_matrix()
856 Every algebra decomposes trivially with respect to its identity
859 sage: set_random_seed()
860 sage: J = random_eja()
861 sage: J0,J5,J1 = J.peirce_decomposition(J.one())
862 sage: J0.dimension() == 0 and J5.dimension() == 0
864 sage: J1.superalgebra() == J and J1.dimension() == J.dimension()
867 The decomposition is into eigenspaces, and its components are
868 therefore necessarily orthogonal. Moreover, the identity
869 elements in the two subalgebras are the projections onto their
870 respective subspaces of the superalgebra's identity element::
872 sage: set_random_seed()
873 sage: J = random_eja()
874 sage: x = J.random_element()
875 sage: if not J.is_trivial():
876 ....: while x.is_nilpotent():
877 ....: x = J.random_element()
878 sage: c = x.subalgebra_idempotent()
879 sage: J0,J5,J1 = J.peirce_decomposition(c)
881 sage: for (w,y,z) in zip(J0.basis(), J5.basis(), J1.basis()):
882 ....: w = w.superalgebra_element()
883 ....: y = J.from_vector(y)
884 ....: z = z.superalgebra_element()
885 ....: ipsum += w.inner_product(y).abs()
886 ....: ipsum += w.inner_product(z).abs()
887 ....: ipsum += y.inner_product(z).abs()
890 sage: J1(c) == J1.one()
892 sage: J0(J.one() - c) == J0.one()
896 if not c
.is_idempotent():
897 raise ValueError("element is not idempotent: %s" % c
)
899 from mjo
.eja
.eja_subalgebra
import FiniteDimensionalEuclideanJordanSubalgebra
901 # Default these to what they should be if they turn out to be
902 # trivial, because eigenspaces_left() won't return eigenvalues
903 # corresponding to trivial spaces (e.g. it returns only the
904 # eigenspace corresponding to lambda=1 if you take the
905 # decomposition relative to the identity element).
906 trivial
= FiniteDimensionalEuclideanJordanSubalgebra(self
, ())
907 J0
= trivial
# eigenvalue zero
908 J5
= VectorSpace(self
.base_ring(), 0) # eigenvalue one-half
909 J1
= trivial
# eigenvalue one
911 for (eigval
, eigspace
) in c
.operator().matrix().right_eigenspaces():
912 if eigval
== ~
(self
.base_ring()(2)):
915 gens
= tuple( self
.from_vector(b
) for b
in eigspace
.basis() )
916 subalg
= FiniteDimensionalEuclideanJordanSubalgebra(self
,
924 raise ValueError("unexpected eigenvalue: %s" % eigval
)
929 def random_element(self
, thorough
=False):
931 Return a random element of this algebra.
933 Our algebra superclass method only returns a linear
934 combination of at most two basis elements. We instead
935 want the vector space "random element" method that
936 returns a more diverse selection.
940 - ``thorough`` -- (boolean; default False) whether or not we
941 should generate irrational coefficients for the random
942 element when our base ring is irrational; this slows the
943 algebra operations to a crawl, but any truly random method
947 # For a general base ring... maybe we can trust this to do the
948 # right thing? Unlikely, but.
949 V
= self
.vector_space()
950 v
= V
.random_element()
952 if self
.base_ring() is AA
:
953 # The "random element" method of the algebraic reals is
954 # stupid at the moment, and only returns integers between
955 # -2 and 2, inclusive:
957 # https://trac.sagemath.org/ticket/30875
959 # Instead, we implement our own "random vector" method,
960 # and then coerce that into the algebra. We use the vector
961 # space degree here instead of the dimension because a
962 # subalgebra could (for example) be spanned by only two
963 # vectors, each with five coordinates. We need to
964 # generate all five coordinates.
966 v
*= QQbar
.random_element().real()
968 v
*= QQ
.random_element()
970 return self
.from_vector(V
.coordinate_vector(v
))
972 def random_elements(self
, count
, thorough
=False):
974 Return ``count`` random elements as a tuple.
978 - ``thorough`` -- (boolean; default False) whether or not we
979 should generate irrational coefficients for the random
980 elements when our base ring is irrational; this slows the
981 algebra operations to a crawl, but any truly random method
986 sage: from mjo.eja.eja_algebra import JordanSpinEJA
990 sage: J = JordanSpinEJA(3)
991 sage: x,y,z = J.random_elements(3)
992 sage: all( [ x in J, y in J, z in J ])
994 sage: len( J.random_elements(10) ) == 10
998 return tuple( self
.random_element(thorough
)
999 for idx
in range(count
) )
1003 def _charpoly_coefficients(self
):
1005 The `r` polynomial coefficients of the "characteristic polynomial
1008 n
= self
.dimension()
1009 R
= self
.coordinate_polynomial_ring()
1011 F
= R
.fraction_field()
1014 # From a result in my book, these are the entries of the
1015 # basis representation of L_x.
1016 return sum( vars[k
]*self
.monomial(k
).operator().matrix()[i
,j
]
1019 L_x
= matrix(F
, n
, n
, L_x_i_j
)
1022 if self
.rank
.is_in_cache():
1024 # There's no need to pad the system with redundant
1025 # columns if we *know* they'll be redundant.
1028 # Compute an extra power in case the rank is equal to
1029 # the dimension (otherwise, we would stop at x^(r-1)).
1030 x_powers
= [ (L_x
**k
)*self
.one().to_vector()
1031 for k
in range(n
+1) ]
1032 A
= matrix
.column(F
, x_powers
[:n
])
1033 AE
= A
.extended_echelon_form()
1040 # The theory says that only the first "r" coefficients are
1041 # nonzero, and they actually live in the original polynomial
1042 # ring and not the fraction field. We negate them because
1043 # in the actual characteristic polynomial, they get moved
1044 # to the other side where x^r lives.
1045 return -A_rref
.solve_right(E
*b
).change_ring(R
)[:r
]
1050 Return the rank of this EJA.
1052 This is a cached method because we know the rank a priori for
1053 all of the algebras we can construct. Thus we can avoid the
1054 expensive ``_charpoly_coefficients()`` call unless we truly
1055 need to compute the whole characteristic polynomial.
1059 sage: from mjo.eja.eja_algebra import (HadamardEJA,
1060 ....: JordanSpinEJA,
1061 ....: RealSymmetricEJA,
1062 ....: ComplexHermitianEJA,
1063 ....: QuaternionHermitianEJA,
1068 The rank of the Jordan spin algebra is always two::
1070 sage: JordanSpinEJA(2).rank()
1072 sage: JordanSpinEJA(3).rank()
1074 sage: JordanSpinEJA(4).rank()
1077 The rank of the `n`-by-`n` Hermitian real, complex, or
1078 quaternion matrices is `n`::
1080 sage: RealSymmetricEJA(4).rank()
1082 sage: ComplexHermitianEJA(3).rank()
1084 sage: QuaternionHermitianEJA(2).rank()
1089 Ensure that every EJA that we know how to construct has a
1090 positive integer rank, unless the algebra is trivial in
1091 which case its rank will be zero::
1093 sage: set_random_seed()
1094 sage: J = random_eja()
1098 sage: r > 0 or (r == 0 and J.is_trivial())
1101 Ensure that computing the rank actually works, since the ranks
1102 of all simple algebras are known and will be cached by default::
1104 sage: set_random_seed() # long time
1105 sage: J = random_eja() # long time
1106 sage: caches = J.rank() # long time
1107 sage: J.rank.clear_cache() # long time
1108 sage: J.rank() == cached # long time
1112 return len(self
._charpoly
_coefficients
())
1115 def vector_space(self
):
1117 Return the vector space that underlies this algebra.
1121 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1125 sage: J = RealSymmetricEJA(2)
1126 sage: J.vector_space()
1127 Vector space of dimension 3 over...
1130 return self
.zero().to_vector().parent().ambient_vector_space()
1133 Element
= FiniteDimensionalEuclideanJordanAlgebraElement
1135 class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra
):
1137 New class for algebras whose supplied basis elements have all rational entries.
1141 sage: from mjo.eja.eja_algebra import BilinearFormEJA
1145 The supplied basis is orthonormalized by default::
1147 sage: B = matrix(QQ, [[1, 0, 0], [0, 25, -32], [0, -32, 41]])
1148 sage: J = BilinearFormEJA(B)
1149 sage: J.matrix_basis()
1162 orthonormalize
=True,
1169 vector_basis
= basis
1171 from sage
.structure
.element
import is_Matrix
1172 basis_is_matrices
= False
1176 if is_Matrix(basis
[0]):
1177 basis_is_matrices
= True
1178 from mjo
.eja
.eja_utils
import _vec2mat
1179 vector_basis
= tuple( map(_mat2vec
,basis
) )
1180 degree
= basis
[0].nrows()**2
1182 degree
= basis
[0].degree()
1184 V
= VectorSpace(field
, degree
)
1186 # If we were asked to orthonormalize, and if the orthonormal
1187 # basis is different from the given one, then we also want to
1188 # compute multiplication and inner-product tables for the
1189 # deorthonormalized basis. These can be used later to
1190 # construct a deorthonormalized copy of this algebra over QQ
1191 # in which several operations are much faster.
1192 self
._deortho
_multiplication
_table
= None
1193 self
._deortho
_inner
_product
_table
= None
1196 # Compute the deorthonormalized tables before we orthonormalize
1198 W
= V
.span_of_basis( vector_basis
)
1201 # If the superclass constructor is going to verify the
1202 # symmetry of this table, it has better at least be
1204 self
._deortho
_multiplication
_table
= [ [0 for j
in range(n
)]
1206 self
._deortho
_inner
_product
_table
= [ [0 for j
in range(n
)]
1209 self
._deortho
_multiplication
_table
= [ [0 for j
in range(i
+1)]
1211 self
._deortho
_inner
_product
_table
= [ [0 for j
in range(i
+1)]
1214 # Note: the Jordan and inner-products are defined in terms
1215 # of the ambient basis. It's important that their arguments
1216 # are in ambient coordinates as well.
1218 for j
in range(i
+1):
1219 # given basis w.r.t. ambient coords
1220 q_i
= vector_basis
[i
]
1221 q_j
= vector_basis
[j
]
1223 if basis_is_matrices
:
1227 elt
= jordan_product(q_i
, q_j
)
1228 ip
= inner_product(q_i
, q_j
)
1230 if basis_is_matrices
:
1231 # do another mat2vec because the multiplication
1232 # table is in terms of vectors
1235 elt
= W
.coordinate_vector(elt
)
1236 self
._deortho
_multiplication
_table
[i
][j
] = elt
1237 self
._deortho
_inner
_product
_table
[i
][j
] = ip
1239 # The tables are square if we're verifying that they
1241 self
._deortho
_multiplication
_table
[j
][i
] = elt
1242 self
._deortho
_inner
_product
_table
[j
][i
] = ip
1244 if self
._deortho
_multiplication
_table
is not None:
1245 self
._deortho
_multiplication
_table
= tuple(map(tuple, self
._deortho
_multiplication
_table
))
1246 if self
._deortho
_inner
_product
_table
is not None:
1247 self
._deortho
_inner
_product
_table
= tuple(map(tuple, self
._deortho
_inner
_product
_table
))
1249 # We overwrite the name "vector_basis" in a second, but never modify it
1250 # in place, to this effectively makes a copy of it.
1251 deortho_vector_basis
= vector_basis
1252 self
._deortho
_matrix
= None
1255 from mjo
.eja
.eja_utils
import gram_schmidt
1256 if basis_is_matrices
:
1257 vector_ip
= lambda x
,y
: inner_product(_vec2mat(x
), _vec2mat(y
))
1258 vector_basis
= gram_schmidt(vector_basis
, vector_ip
)
1260 vector_basis
= gram_schmidt(vector_basis
, inner_product
)
1262 W
= V
.span_of_basis( vector_basis
)
1264 # Normalize the "matrix" basis, too!
1265 basis
= vector_basis
1267 if basis_is_matrices
:
1268 basis
= tuple( map(_vec2mat
,basis
) )
1270 W
= V
.span_of_basis( vector_basis
)
1272 # Now "W" is the vector space of our algebra coordinates. The
1273 # variables "X1", "X2",... refer to the entries of vectors in
1274 # W. Thus to convert back and forth between the orthonormal
1275 # coordinates and the given ones, we need to stick the original
1277 U
= V
.span_of_basis( deortho_vector_basis
)
1278 self
._deortho
_matrix
= matrix( U
.coordinate_vector(q
)
1279 for q
in vector_basis
)
1281 # If the superclass constructor is going to verify the
1282 # symmetry of this table, it has better at least be
1285 mult_table
= [ [0 for j
in range(n
)] for i
in range(n
) ]
1286 ip_table
= [ [0 for j
in range(n
)] for i
in range(n
) ]
1288 mult_table
= [ [0 for j
in range(i
+1)] for i
in range(n
) ]
1289 ip_table
= [ [0 for j
in range(i
+1)] for i
in range(n
) ]
1291 # Note: the Jordan and inner-products are defined in terms
1292 # of the ambient basis. It's important that their arguments
1293 # are in ambient coordinates as well.
1295 for j
in range(i
+1):
1296 # ortho basis w.r.t. ambient coords
1297 q_i
= vector_basis
[i
]
1298 q_j
= vector_basis
[j
]
1300 if basis_is_matrices
:
1304 elt
= jordan_product(q_i
, q_j
)
1305 ip
= inner_product(q_i
, q_j
)
1307 if basis_is_matrices
:
1308 # do another mat2vec because the multiplication
1309 # table is in terms of vectors
1312 elt
= W
.coordinate_vector(elt
)
1313 mult_table
[i
][j
] = elt
1316 # The tables are square if we're verifying that they
1318 mult_table
[j
][i
] = elt
1321 if basis_is_matrices
:
1325 basis
= tuple( x
.column() for x
in basis
)
1327 super().__init
__(field
,
1332 basis
, # matrix basis
1337 def _charpoly_coefficients(self
):
1341 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
1342 ....: JordanSpinEJA)
1346 The base ring of the resulting polynomial coefficients is what
1347 it should be, and not the rationals (unless the algebra was
1348 already over the rationals)::
1350 sage: J = JordanSpinEJA(3)
1351 sage: J._charpoly_coefficients()
1352 (X1^2 - X2^2 - X3^2, -2*X1)
1353 sage: a0 = J._charpoly_coefficients()[0]
1355 Algebraic Real Field
1356 sage: a0.base_ring()
1357 Algebraic Real Field
1360 if self
.base_ring() is QQ
:
1361 # There's no need to construct *another* algebra over the
1362 # rationals if this one is already over the rationals.
1363 superclass
= super(RationalBasisEuclideanJordanAlgebra
, self
)
1364 return superclass
._charpoly
_coefficients
()
1366 # Do the computation over the rationals. The answer will be
1367 # the same, because all we've done is a change of basis.
1368 J
= FiniteDimensionalEuclideanJordanAlgebra(QQ
,
1369 self
._deortho
_multiplication
_table
,
1370 self
._deortho
_inner
_product
_table
)
1372 # Change back from QQ to our real base ring
1373 a
= ( a_i
.change_ring(self
.base_ring())
1374 for a_i
in J
._charpoly
_coefficients
() )
1376 # Now convert the coordinate variables back to the
1377 # deorthonormalized ones.
1378 R
= self
.coordinate_polynomial_ring()
1379 from sage
.modules
.free_module_element
import vector
1380 X
= vector(R
, R
.gens())
1381 BX
= self
._deortho
_matrix
*X
1383 subs_dict
= { X[i]: BX[i] for i in range(len(X)) }
1384 return tuple( a_i
.subs(subs_dict
) for a_i
in a
)
1386 class ConcreteEuclideanJordanAlgebra(RationalBasisEuclideanJordanAlgebra
):
1388 A class for the Euclidean Jordan algebras that we know by name.
1390 These are the Jordan algebras whose basis, multiplication table,
1391 rank, and so on are known a priori. More to the point, they are
1392 the Euclidean Jordan algebras for which we are able to conjure up
1393 a "random instance."
1397 sage: from mjo.eja.eja_algebra import ConcreteEuclideanJordanAlgebra
1401 Our basis is normalized with respect to the algebra's inner
1402 product, unless we specify otherwise::
1404 sage: set_random_seed()
1405 sage: J = ConcreteEuclideanJordanAlgebra.random_instance()
1406 sage: all( b.norm() == 1 for b in J.gens() )
1409 Since our basis is orthonormal with respect to the algebra's inner
1410 product, and since we know that this algebra is an EJA, any
1411 left-multiplication operator's matrix will be symmetric because
1412 natural->EJA basis representation is an isometry and within the
1413 EJA the operator is self-adjoint by the Jordan axiom::
1415 sage: set_random_seed()
1416 sage: J = ConcreteEuclideanJordanAlgebra.random_instance()
1417 sage: x = J.random_element()
1418 sage: x.operator().is_self_adjoint()
1423 def _max_random_instance_size():
1425 Return an integer "size" that is an upper bound on the size of
1426 this algebra when it is used in a random test
1427 case. Unfortunately, the term "size" is ambiguous -- when
1428 dealing with `R^n` under either the Hadamard or Jordan spin
1429 product, the "size" refers to the dimension `n`. When dealing
1430 with a matrix algebra (real symmetric or complex/quaternion
1431 Hermitian), it refers to the size of the matrix, which is far
1432 less than the dimension of the underlying vector space.
1434 This method must be implemented in each subclass.
1436 raise NotImplementedError
1439 def random_instance(cls
, *args
, **kwargs
):
1441 Return a random instance of this type of algebra.
1443 This method should be implemented in each subclass.
1445 from sage
.misc
.prandom
import choice
1446 eja_class
= choice(cls
.__subclasses
__())
1448 # These all bubble up to the RationalBasisEuclideanJordanAlgebra
1449 # superclass constructor, so any (kw)args valid there are also
1451 return eja_class
.random_instance(*args
, **kwargs
)
1454 class MatrixEuclideanJordanAlgebra
:
1458 Embed the matrix ``M`` into a space of real matrices.
1460 The matrix ``M`` can have entries in any field at the moment:
1461 the real numbers, complex numbers, or quaternions. And although
1462 they are not a field, we can probably support octonions at some
1463 point, too. This function returns a real matrix that "acts like"
1464 the original with respect to matrix multiplication; i.e.
1466 real_embed(M*N) = real_embed(M)*real_embed(N)
1469 raise NotImplementedError
1473 def real_unembed(M
):
1475 The inverse of :meth:`real_embed`.
1477 raise NotImplementedError
1480 def jordan_product(X
,Y
):
1481 return (X
*Y
+ Y
*X
)/2
1484 def trace_inner_product(cls
,X
,Y
):
1485 Xu
= cls
.real_unembed(X
)
1486 Yu
= cls
.real_unembed(Y
)
1487 tr
= (Xu
*Yu
).trace()
1490 # Works in QQ, AA, RDF, et cetera.
1492 except AttributeError:
1493 # A quaternion doesn't have a real() method, but does
1494 # have coefficient_tuple() method that returns the
1495 # coefficients of 1, i, j, and k -- in that order.
1496 return tr
.coefficient_tuple()[0]
1499 class RealMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1503 The identity function, for embedding real matrices into real
1509 def real_unembed(M
):
1511 The identity function, for unembedding real matrices from real
1517 class RealSymmetricEJA(ConcreteEuclideanJordanAlgebra
,
1518 RealMatrixEuclideanJordanAlgebra
):
1520 The rank-n simple EJA consisting of real symmetric n-by-n
1521 matrices, the usual symmetric Jordan product, and the trace inner
1522 product. It has dimension `(n^2 + n)/2` over the reals.
1526 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1530 sage: J = RealSymmetricEJA(2)
1531 sage: e0, e1, e2 = J.gens()
1539 In theory, our "field" can be any subfield of the reals::
1541 sage: RealSymmetricEJA(2, RDF)
1542 Euclidean Jordan algebra of dimension 3 over Real Double Field
1543 sage: RealSymmetricEJA(2, RR)
1544 Euclidean Jordan algebra of dimension 3 over Real Field with
1545 53 bits of precision
1549 The dimension of this algebra is `(n^2 + n) / 2`::
1551 sage: set_random_seed()
1552 sage: n_max = RealSymmetricEJA._max_random_instance_size()
1553 sage: n = ZZ.random_element(1, n_max)
1554 sage: J = RealSymmetricEJA(n)
1555 sage: J.dimension() == (n^2 + n)/2
1558 The Jordan multiplication is what we think it is::
1560 sage: set_random_seed()
1561 sage: J = RealSymmetricEJA.random_instance()
1562 sage: x,y = J.random_elements(2)
1563 sage: actual = (x*y).to_matrix()
1564 sage: X = x.to_matrix()
1565 sage: Y = y.to_matrix()
1566 sage: expected = (X*Y + Y*X)/2
1567 sage: actual == expected
1569 sage: J(expected) == x*y
1572 We can change the generator prefix::
1574 sage: RealSymmetricEJA(3, prefix='q').gens()
1575 (q0, q1, q2, q3, q4, q5)
1577 We can construct the (trivial) algebra of rank zero::
1579 sage: RealSymmetricEJA(0)
1580 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1584 def _denormalized_basis(cls
, n
, field
):
1586 Return a basis for the space of real symmetric n-by-n matrices.
1590 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1594 sage: set_random_seed()
1595 sage: n = ZZ.random_element(1,5)
1596 sage: B = RealSymmetricEJA._denormalized_basis(n,QQ)
1597 sage: all( M.is_symmetric() for M in B)
1601 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1605 for j
in range(i
+1):
1606 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
1610 Sij
= Eij
+ Eij
.transpose()
1616 def _max_random_instance_size():
1617 return 4 # Dimension 10
1620 def random_instance(cls
, field
=AA
, **kwargs
):
1622 Return a random instance of this type of algebra.
1624 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1625 return cls(n
, field
, **kwargs
)
1627 def __init__(self
, n
, field
=AA
, **kwargs
):
1628 basis
= self
._denormalized
_basis
(n
, field
)
1629 super(RealSymmetricEJA
, self
).__init
__(field
,
1631 self
.jordan_product
,
1632 self
.trace_inner_product
,
1634 self
.rank
.set_cache(n
)
1635 self
.one
.set_cache(self(matrix
.identity(field
,n
)))
1638 class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1642 Embed the n-by-n complex matrix ``M`` into the space of real
1643 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1644 bi` to the block matrix ``[[a,b],[-b,a]]``.
1648 sage: from mjo.eja.eja_algebra import \
1649 ....: ComplexMatrixEuclideanJordanAlgebra
1653 sage: F = QuadraticField(-1, 'I')
1654 sage: x1 = F(4 - 2*i)
1655 sage: x2 = F(1 + 2*i)
1658 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1659 sage: ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1668 Embedding is a homomorphism (isomorphism, in fact)::
1670 sage: set_random_seed()
1671 sage: n = ZZ.random_element(3)
1672 sage: F = QuadraticField(-1, 'I')
1673 sage: X = random_matrix(F, n)
1674 sage: Y = random_matrix(F, n)
1675 sage: Xe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X)
1676 sage: Ye = ComplexMatrixEuclideanJordanAlgebra.real_embed(Y)
1677 sage: XYe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1684 raise ValueError("the matrix 'M' must be square")
1686 # We don't need any adjoined elements...
1687 field
= M
.base_ring().base_ring()
1691 a
= z
.list()[0] # real part, I guess
1692 b
= z
.list()[1] # imag part, I guess
1693 blocks
.append(matrix(field
, 2, [[a
,b
],[-b
,a
]]))
1695 return matrix
.block(field
, n
, blocks
)
1699 def real_unembed(M
):
1701 The inverse of _embed_complex_matrix().
1705 sage: from mjo.eja.eja_algebra import \
1706 ....: ComplexMatrixEuclideanJordanAlgebra
1710 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1711 ....: [-2, 1, -4, 3],
1712 ....: [ 9, 10, 11, 12],
1713 ....: [-10, 9, -12, 11] ])
1714 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(A)
1716 [ 10*I + 9 12*I + 11]
1720 Unembedding is the inverse of embedding::
1722 sage: set_random_seed()
1723 sage: F = QuadraticField(-1, 'I')
1724 sage: M = random_matrix(F, 3)
1725 sage: Me = ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1726 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1732 raise ValueError("the matrix 'M' must be square")
1733 if not n
.mod(2).is_zero():
1734 raise ValueError("the matrix 'M' must be a complex embedding")
1736 # If "M" was normalized, its base ring might have roots
1737 # adjoined and they can stick around after unembedding.
1738 field
= M
.base_ring()
1739 R
= PolynomialRing(field
, 'z')
1742 # Sage doesn't know how to embed AA into QQbar, i.e. how
1743 # to adjoin sqrt(-1) to AA.
1746 F
= field
.extension(z
**2 + 1, 'I', embedding
=CLF(-1).sqrt())
1749 # Go top-left to bottom-right (reading order), converting every
1750 # 2-by-2 block we see to a single complex element.
1752 for k
in range(n
/2):
1753 for j
in range(n
/2):
1754 submat
= M
[2*k
:2*k
+2,2*j
:2*j
+2]
1755 if submat
[0,0] != submat
[1,1]:
1756 raise ValueError('bad on-diagonal submatrix')
1757 if submat
[0,1] != -submat
[1,0]:
1758 raise ValueError('bad off-diagonal submatrix')
1759 z
= submat
[0,0] + submat
[0,1]*i
1762 return matrix(F
, n
/2, elements
)
1766 def trace_inner_product(cls
,X
,Y
):
1768 Compute a matrix inner product in this algebra directly from
1773 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1777 This gives the same answer as the slow, default method implemented
1778 in :class:`MatrixEuclideanJordanAlgebra`::
1780 sage: set_random_seed()
1781 sage: J = ComplexHermitianEJA.random_instance()
1782 sage: x,y = J.random_elements(2)
1783 sage: Xe = x.to_matrix()
1784 sage: Ye = y.to_matrix()
1785 sage: X = ComplexHermitianEJA.real_unembed(Xe)
1786 sage: Y = ComplexHermitianEJA.real_unembed(Ye)
1787 sage: expected = (X*Y).trace().real()
1788 sage: actual = ComplexHermitianEJA.trace_inner_product(Xe,Ye)
1789 sage: actual == expected
1793 return RealMatrixEuclideanJordanAlgebra
.trace_inner_product(X
,Y
)/2
1796 class ComplexHermitianEJA(ConcreteEuclideanJordanAlgebra
,
1797 ComplexMatrixEuclideanJordanAlgebra
):
1799 The rank-n simple EJA consisting of complex Hermitian n-by-n
1800 matrices over the real numbers, the usual symmetric Jordan product,
1801 and the real-part-of-trace inner product. It has dimension `n^2` over
1806 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1810 In theory, our "field" can be any subfield of the reals::
1812 sage: ComplexHermitianEJA(2, RDF)
1813 Euclidean Jordan algebra of dimension 4 over Real Double Field
1814 sage: ComplexHermitianEJA(2, RR)
1815 Euclidean Jordan algebra of dimension 4 over Real Field with
1816 53 bits of precision
1820 The dimension of this algebra is `n^2`::
1822 sage: set_random_seed()
1823 sage: n_max = ComplexHermitianEJA._max_random_instance_size()
1824 sage: n = ZZ.random_element(1, n_max)
1825 sage: J = ComplexHermitianEJA(n)
1826 sage: J.dimension() == n^2
1829 The Jordan multiplication is what we think it is::
1831 sage: set_random_seed()
1832 sage: J = ComplexHermitianEJA.random_instance()
1833 sage: x,y = J.random_elements(2)
1834 sage: actual = (x*y).to_matrix()
1835 sage: X = x.to_matrix()
1836 sage: Y = y.to_matrix()
1837 sage: expected = (X*Y + Y*X)/2
1838 sage: actual == expected
1840 sage: J(expected) == x*y
1843 We can change the generator prefix::
1845 sage: ComplexHermitianEJA(2, prefix='z').gens()
1848 We can construct the (trivial) algebra of rank zero::
1850 sage: ComplexHermitianEJA(0)
1851 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1856 def _denormalized_basis(cls
, n
, field
):
1858 Returns a basis for the space of complex Hermitian n-by-n matrices.
1860 Why do we embed these? Basically, because all of numerical linear
1861 algebra assumes that you're working with vectors consisting of `n`
1862 entries from a field and scalars from the same field. There's no way
1863 to tell SageMath that (for example) the vectors contain complex
1864 numbers, while the scalar field is real.
1868 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1872 sage: set_random_seed()
1873 sage: n = ZZ.random_element(1,5)
1874 sage: field = QuadraticField(2, 'sqrt2')
1875 sage: B = ComplexHermitianEJA._denormalized_basis(n, field)
1876 sage: all( M.is_symmetric() for M in B)
1880 R
= PolynomialRing(field
, 'z')
1882 F
= field
.extension(z
**2 + 1, 'I')
1885 # This is like the symmetric case, but we need to be careful:
1887 # * We want conjugate-symmetry, not just symmetry.
1888 # * The diagonal will (as a result) be real.
1892 for j
in range(i
+1):
1893 Eij
= matrix(F
, n
, lambda k
,l
: k
==i
and l
==j
)
1895 Sij
= cls
.real_embed(Eij
)
1898 # The second one has a minus because it's conjugated.
1899 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
1901 Sij_imag
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
1904 # Since we embedded these, we can drop back to the "field" that we
1905 # started with instead of the complex extension "F".
1906 return tuple( s
.change_ring(field
) for s
in S
)
1909 def __init__(self
, n
, field
=AA
, **kwargs
):
1910 basis
= self
._denormalized
_basis
(n
,field
)
1911 super(ComplexHermitianEJA
, self
).__init
__(field
,
1913 self
.jordan_product
,
1914 self
.trace_inner_product
,
1916 self
.rank
.set_cache(n
)
1917 # TODO: pre-cache the identity!
1920 def _max_random_instance_size():
1921 return 3 # Dimension 9
1924 def random_instance(cls
, field
=AA
, **kwargs
):
1926 Return a random instance of this type of algebra.
1928 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1929 return cls(n
, field
, **kwargs
)
1931 class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1935 Embed the n-by-n quaternion matrix ``M`` into the space of real
1936 matrices of size 4n-by-4n by first sending each quaternion entry `z
1937 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
1938 c+di],[-c + di, a-bi]]`, and then embedding those into a real
1943 sage: from mjo.eja.eja_algebra import \
1944 ....: QuaternionMatrixEuclideanJordanAlgebra
1948 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1949 sage: i,j,k = Q.gens()
1950 sage: x = 1 + 2*i + 3*j + 4*k
1951 sage: M = matrix(Q, 1, [[x]])
1952 sage: QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1958 Embedding is a homomorphism (isomorphism, in fact)::
1960 sage: set_random_seed()
1961 sage: n = ZZ.random_element(2)
1962 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1963 sage: X = random_matrix(Q, n)
1964 sage: Y = random_matrix(Q, n)
1965 sage: Xe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X)
1966 sage: Ye = QuaternionMatrixEuclideanJordanAlgebra.real_embed(Y)
1967 sage: XYe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1972 quaternions
= M
.base_ring()
1975 raise ValueError("the matrix 'M' must be square")
1977 F
= QuadraticField(-1, 'I')
1982 t
= z
.coefficient_tuple()
1987 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
1988 [-c
+ d
*i
, a
- b
*i
]])
1989 realM
= ComplexMatrixEuclideanJordanAlgebra
.real_embed(cplxM
)
1990 blocks
.append(realM
)
1992 # We should have real entries by now, so use the realest field
1993 # we've got for the return value.
1994 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
1999 def real_unembed(M
):
2001 The inverse of _embed_quaternion_matrix().
2005 sage: from mjo.eja.eja_algebra import \
2006 ....: QuaternionMatrixEuclideanJordanAlgebra
2010 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
2011 ....: [-2, 1, -4, 3],
2012 ....: [-3, 4, 1, -2],
2013 ....: [-4, -3, 2, 1]])
2014 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(M)
2015 [1 + 2*i + 3*j + 4*k]
2019 Unembedding is the inverse of embedding::
2021 sage: set_random_seed()
2022 sage: Q = QuaternionAlgebra(QQ, -1, -1)
2023 sage: M = random_matrix(Q, 3)
2024 sage: Me = QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
2025 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
2031 raise ValueError("the matrix 'M' must be square")
2032 if not n
.mod(4).is_zero():
2033 raise ValueError("the matrix 'M' must be a quaternion embedding")
2035 # Use the base ring of the matrix to ensure that its entries can be
2036 # multiplied by elements of the quaternion algebra.
2037 field
= M
.base_ring()
2038 Q
= QuaternionAlgebra(field
,-1,-1)
2041 # Go top-left to bottom-right (reading order), converting every
2042 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
2045 for l
in range(n
/4):
2046 for m
in range(n
/4):
2047 submat
= ComplexMatrixEuclideanJordanAlgebra
.real_unembed(
2048 M
[4*l
:4*l
+4,4*m
:4*m
+4] )
2049 if submat
[0,0] != submat
[1,1].conjugate():
2050 raise ValueError('bad on-diagonal submatrix')
2051 if submat
[0,1] != -submat
[1,0].conjugate():
2052 raise ValueError('bad off-diagonal submatrix')
2053 z
= submat
[0,0].real()
2054 z
+= submat
[0,0].imag()*i
2055 z
+= submat
[0,1].real()*j
2056 z
+= submat
[0,1].imag()*k
2059 return matrix(Q
, n
/4, elements
)
2063 def trace_inner_product(cls
,X
,Y
):
2065 Compute a matrix inner product in this algebra directly from
2070 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2074 This gives the same answer as the slow, default method implemented
2075 in :class:`MatrixEuclideanJordanAlgebra`::
2077 sage: set_random_seed()
2078 sage: J = QuaternionHermitianEJA.random_instance()
2079 sage: x,y = J.random_elements(2)
2080 sage: Xe = x.to_matrix()
2081 sage: Ye = y.to_matrix()
2082 sage: X = QuaternionHermitianEJA.real_unembed(Xe)
2083 sage: Y = QuaternionHermitianEJA.real_unembed(Ye)
2084 sage: expected = (X*Y).trace().coefficient_tuple()[0]
2085 sage: actual = QuaternionHermitianEJA.trace_inner_product(Xe,Ye)
2086 sage: actual == expected
2090 return RealMatrixEuclideanJordanAlgebra
.trace_inner_product(X
,Y
)/4
2093 class QuaternionHermitianEJA(ConcreteEuclideanJordanAlgebra
,
2094 QuaternionMatrixEuclideanJordanAlgebra
):
2096 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
2097 matrices, the usual symmetric Jordan product, and the
2098 real-part-of-trace inner product. It has dimension `2n^2 - n` over
2103 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2107 In theory, our "field" can be any subfield of the reals::
2109 sage: QuaternionHermitianEJA(2, RDF)
2110 Euclidean Jordan algebra of dimension 6 over Real Double Field
2111 sage: QuaternionHermitianEJA(2, RR)
2112 Euclidean Jordan algebra of dimension 6 over Real Field with
2113 53 bits of precision
2117 The dimension of this algebra is `2*n^2 - n`::
2119 sage: set_random_seed()
2120 sage: n_max = QuaternionHermitianEJA._max_random_instance_size()
2121 sage: n = ZZ.random_element(1, n_max)
2122 sage: J = QuaternionHermitianEJA(n)
2123 sage: J.dimension() == 2*(n^2) - n
2126 The Jordan multiplication is what we think it is::
2128 sage: set_random_seed()
2129 sage: J = QuaternionHermitianEJA.random_instance()
2130 sage: x,y = J.random_elements(2)
2131 sage: actual = (x*y).to_matrix()
2132 sage: X = x.to_matrix()
2133 sage: Y = y.to_matrix()
2134 sage: expected = (X*Y + Y*X)/2
2135 sage: actual == expected
2137 sage: J(expected) == x*y
2140 We can change the generator prefix::
2142 sage: QuaternionHermitianEJA(2, prefix='a').gens()
2143 (a0, a1, a2, a3, a4, a5)
2145 We can construct the (trivial) algebra of rank zero::
2147 sage: QuaternionHermitianEJA(0)
2148 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2152 def _denormalized_basis(cls
, n
, field
):
2154 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
2156 Why do we embed these? Basically, because all of numerical
2157 linear algebra assumes that you're working with vectors consisting
2158 of `n` entries from a field and scalars from the same field. There's
2159 no way to tell SageMath that (for example) the vectors contain
2160 complex numbers, while the scalar field is real.
2164 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2168 sage: set_random_seed()
2169 sage: n = ZZ.random_element(1,5)
2170 sage: B = QuaternionHermitianEJA._denormalized_basis(n,QQ)
2171 sage: all( M.is_symmetric() for M in B )
2175 Q
= QuaternionAlgebra(QQ
,-1,-1)
2178 # This is like the symmetric case, but we need to be careful:
2180 # * We want conjugate-symmetry, not just symmetry.
2181 # * The diagonal will (as a result) be real.
2185 for j
in range(i
+1):
2186 Eij
= matrix(Q
, n
, lambda k
,l
: k
==i
and l
==j
)
2188 Sij
= cls
.real_embed(Eij
)
2191 # The second, third, and fourth ones have a minus
2192 # because they're conjugated.
2193 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
2195 Sij_I
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
2197 Sij_J
= cls
.real_embed(J
*Eij
- J
*Eij
.transpose())
2199 Sij_K
= cls
.real_embed(K
*Eij
- K
*Eij
.transpose())
2202 # Since we embedded these, we can drop back to the "field" that we
2203 # started with instead of the quaternion algebra "Q".
2204 return tuple( s
.change_ring(field
) for s
in S
)
2207 def __init__(self
, n
, field
=AA
, **kwargs
):
2208 basis
= self
._denormalized
_basis
(n
,field
)
2209 super(QuaternionHermitianEJA
, self
).__init
__(field
,
2211 self
.jordan_product
,
2212 self
.trace_inner_product
,
2214 self
.rank
.set_cache(n
)
2215 # TODO: cache one()!
2218 def _max_random_instance_size():
2220 The maximum rank of a random QuaternionHermitianEJA.
2222 return 2 # Dimension 6
2225 def random_instance(cls
, field
=AA
, **kwargs
):
2227 Return a random instance of this type of algebra.
2229 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2230 return cls(n
, field
, **kwargs
)
2233 class HadamardEJA(ConcreteEuclideanJordanAlgebra
):
2235 Return the Euclidean Jordan Algebra corresponding to the set
2236 `R^n` under the Hadamard product.
2238 Note: this is nothing more than the Cartesian product of ``n``
2239 copies of the spin algebra. Once Cartesian product algebras
2240 are implemented, this can go.
2244 sage: from mjo.eja.eja_algebra import HadamardEJA
2248 This multiplication table can be verified by hand::
2250 sage: J = HadamardEJA(3)
2251 sage: e0,e1,e2 = J.gens()
2267 We can change the generator prefix::
2269 sage: HadamardEJA(3, prefix='r').gens()
2273 def __init__(self
, n
, field
=AA
, **kwargs
):
2274 V
= VectorSpace(field
, n
)
2277 def jordan_product(x
,y
):
2278 return V([ xi
*yi
for (xi
,yi
) in zip(x
,y
) ])
2279 def inner_product(x
,y
):
2280 return x
.inner_product(y
)
2282 super(HadamardEJA
, self
).__init
__(field
,
2287 self
.rank
.set_cache(n
)
2290 self
.one
.set_cache( self
.zero() )
2292 self
.one
.set_cache( sum(self
.gens()) )
2295 def _max_random_instance_size():
2297 The maximum dimension of a random HadamardEJA.
2302 def random_instance(cls
, field
=AA
, **kwargs
):
2304 Return a random instance of this type of algebra.
2306 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2307 return cls(n
, field
, **kwargs
)
2310 class BilinearFormEJA(ConcreteEuclideanJordanAlgebra
):
2312 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2313 with the half-trace inner product and jordan product ``x*y =
2314 (<Bx,y>,y_bar>, x0*y_bar + y0*x_bar)`` where `B = 1 \times B22` is
2315 a symmetric positive-definite "bilinear form" matrix. Its
2316 dimension is the size of `B`, and it has rank two in dimensions
2317 larger than two. It reduces to the ``JordanSpinEJA`` when `B` is
2318 the identity matrix of order ``n``.
2320 We insist that the one-by-one upper-left identity block of `B` be
2321 passed in as well so that we can be passed a matrix of size zero
2322 to construct a trivial algebra.
2326 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
2327 ....: JordanSpinEJA)
2331 When no bilinear form is specified, the identity matrix is used,
2332 and the resulting algebra is the Jordan spin algebra::
2334 sage: B = matrix.identity(AA,3)
2335 sage: J0 = BilinearFormEJA(B)
2336 sage: J1 = JordanSpinEJA(3)
2337 sage: J0.multiplication_table() == J0.multiplication_table()
2340 An error is raised if the matrix `B` does not correspond to a
2341 positive-definite bilinear form::
2343 sage: B = matrix.random(QQ,2,3)
2344 sage: J = BilinearFormEJA(B)
2345 Traceback (most recent call last):
2347 ValueError: bilinear form is not positive-definite
2348 sage: B = matrix.zero(QQ,3)
2349 sage: J = BilinearFormEJA(B)
2350 Traceback (most recent call last):
2352 ValueError: bilinear form is not positive-definite
2356 We can create a zero-dimensional algebra::
2358 sage: B = matrix.identity(AA,0)
2359 sage: J = BilinearFormEJA(B)
2363 We can check the multiplication condition given in the Jordan, von
2364 Neumann, and Wigner paper (and also discussed on my "On the
2365 symmetry..." paper). Note that this relies heavily on the standard
2366 choice of basis, as does anything utilizing the bilinear form
2367 matrix. We opt not to orthonormalize the basis, because if we
2368 did, we would have to normalize the `s_{i}` in a similar manner::
2370 sage: set_random_seed()
2371 sage: n = ZZ.random_element(5)
2372 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2373 sage: B11 = matrix.identity(QQ,1)
2374 sage: B22 = M.transpose()*M
2375 sage: B = block_matrix(2,2,[ [B11,0 ],
2377 sage: J = BilinearFormEJA(B, orthonormalize=False)
2378 sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
2379 sage: V = J.vector_space()
2380 sage: sis = [ J( V([0] + (M.inverse()*ei).list()).column() )
2381 ....: for ei in eis ]
2382 sage: actual = [ sis[i]*sis[j]
2383 ....: for i in range(n-1)
2384 ....: for j in range(n-1) ]
2385 sage: expected = [ J.one() if i == j else J.zero()
2386 ....: for i in range(n-1)
2387 ....: for j in range(n-1) ]
2388 sage: actual == expected
2391 def __init__(self
, B
, field
=AA
, **kwargs
):
2392 if not B
.is_positive_definite():
2393 raise ValueError("bilinear form is not positive-definite")
2396 V
= VectorSpace(field
, n
)
2398 def inner_product(x
,y
):
2399 return (B
*x
).inner_product(y
)
2401 def jordan_product(x
,y
):
2406 z0
= inner_product(x
,y
)
2407 zbar
= y0
*xbar
+ x0
*ybar
2408 return V([z0
] + zbar
.list())
2410 super(BilinearFormEJA
, self
).__init
__(field
,
2416 # The rank of this algebra is two, unless we're in a
2417 # one-dimensional ambient space (because the rank is bounded
2418 # by the ambient dimension).
2419 self
.rank
.set_cache(min(n
,2))
2422 self
.one
.set_cache( self
.zero() )
2424 self
.one
.set_cache( self
.monomial(0) )
2427 def _max_random_instance_size():
2429 The maximum dimension of a random BilinearFormEJA.
2434 def random_instance(cls
, field
=AA
, **kwargs
):
2436 Return a random instance of this algebra.
2438 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2440 B
= matrix
.identity(field
, n
)
2441 return cls(B
, field
, **kwargs
)
2443 B11
= matrix
.identity(field
,1)
2444 M
= matrix
.random(field
, n
-1)
2445 I
= matrix
.identity(field
, n
-1)
2446 alpha
= field
.zero()
2447 while alpha
.is_zero():
2448 alpha
= field
.random_element().abs()
2449 B22
= M
.transpose()*M
+ alpha
*I
2451 from sage
.matrix
.special
import block_matrix
2452 B
= block_matrix(2,2, [ [B11
, ZZ(0) ],
2455 return cls(B
, field
, **kwargs
)
2458 class JordanSpinEJA(BilinearFormEJA
):
2460 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2461 with the usual inner product and jordan product ``x*y =
2462 (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2467 sage: from mjo.eja.eja_algebra import JordanSpinEJA
2471 This multiplication table can be verified by hand::
2473 sage: J = JordanSpinEJA(4)
2474 sage: e0,e1,e2,e3 = J.gens()
2490 We can change the generator prefix::
2492 sage: JordanSpinEJA(2, prefix='B').gens()
2497 Ensure that we have the usual inner product on `R^n`::
2499 sage: set_random_seed()
2500 sage: J = JordanSpinEJA.random_instance()
2501 sage: x,y = J.random_elements(2)
2502 sage: actual = x.inner_product(y)
2503 sage: expected = x.to_vector().inner_product(y.to_vector())
2504 sage: actual == expected
2508 def __init__(self
, n
, field
=AA
, **kwargs
):
2509 # This is a special case of the BilinearFormEJA with the identity
2510 # matrix as its bilinear form.
2511 B
= matrix
.identity(field
, n
)
2512 super(JordanSpinEJA
, self
).__init
__(B
, field
, **kwargs
)
2515 def _max_random_instance_size():
2517 The maximum dimension of a random JordanSpinEJA.
2522 def random_instance(cls
, field
=AA
, **kwargs
):
2524 Return a random instance of this type of algebra.
2526 Needed here to override the implementation for ``BilinearFormEJA``.
2528 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2529 return cls(n
, field
, **kwargs
)
2532 class TrivialEJA(ConcreteEuclideanJordanAlgebra
):
2534 The trivial Euclidean Jordan algebra consisting of only a zero element.
2538 sage: from mjo.eja.eja_algebra import TrivialEJA
2542 sage: J = TrivialEJA()
2549 sage: 7*J.one()*12*J.one()
2551 sage: J.one().inner_product(J.one())
2553 sage: J.one().norm()
2555 sage: J.one().subalgebra_generated_by()
2556 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2561 def __init__(self
, field
=AA
, **kwargs
):
2562 jordan_product
= lambda x
,y
: x
2563 inner_product
= lambda x
,y
: field(0)
2565 super(TrivialEJA
, self
).__init
__(field
,
2570 # The rank is zero using my definition, namely the dimension of the
2571 # largest subalgebra generated by any element.
2572 self
.rank
.set_cache(0)
2573 self
.one
.set_cache( self
.zero() )
2576 def random_instance(cls
, field
=AA
, **kwargs
):
2577 # We don't take a "size" argument so the superclass method is
2578 # inappropriate for us.
2579 return cls(field
, **kwargs
)
2581 class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra
):
2583 The external (orthogonal) direct sum of two other Euclidean Jordan
2584 algebras. Essentially the Cartesian product of its two factors.
2585 Every Euclidean Jordan algebra decomposes into an orthogonal
2586 direct sum of simple Euclidean Jordan algebras, so no generality
2587 is lost by providing only this construction.
2591 sage: from mjo.eja.eja_algebra import (random_eja,
2593 ....: RealSymmetricEJA,
2598 sage: J1 = HadamardEJA(2)
2599 sage: J2 = RealSymmetricEJA(3)
2600 sage: J = DirectSumEJA(J1,J2)
2608 The external direct sum construction is only valid when the two factors
2609 have the same base ring; an error is raised otherwise::
2611 sage: set_random_seed()
2612 sage: J1 = random_eja(AA)
2613 sage: J2 = random_eja(QQ,orthonormalize=False)
2614 sage: J = DirectSumEJA(J1,J2)
2615 Traceback (most recent call last):
2617 ValueError: algebras must share the same base field
2620 def __init__(self
, J1
, J2
, **kwargs
):
2621 if J1
.base_ring() != J2
.base_ring():
2622 raise ValueError("algebras must share the same base field")
2623 field
= J1
.base_ring()
2625 self
._factors
= (J1
, J2
)
2629 V
= VectorSpace(field
, n
)
2630 mult_table
= [ [ V
.zero() for j
in range(i
+1) ]
2633 for j
in range(i
+1):
2634 p
= (J1
.monomial(i
)*J1
.monomial(j
)).to_vector()
2635 mult_table
[i
][j
] = V(p
.list() + [field
.zero()]*n2
)
2638 for j
in range(i
+1):
2639 p
= (J2
.monomial(i
)*J2
.monomial(j
)).to_vector()
2640 mult_table
[n1
+i
][n1
+j
] = V([field
.zero()]*n1
+ p
.list())
2642 # TODO: build the IP table here from the two constituent IP
2643 # matrices (it'll be block diagonal, I think).
2644 ip_table
= [ [ field
.zero() for j
in range(i
+1) ]
2646 super(DirectSumEJA
, self
).__init
__(field
,
2651 self
.rank
.set_cache(J1
.rank() + J2
.rank())
2656 Return the pair of this algebra's factors.
2660 sage: from mjo.eja.eja_algebra import (HadamardEJA,
2661 ....: JordanSpinEJA,
2666 sage: J1 = HadamardEJA(2,QQ)
2667 sage: J2 = JordanSpinEJA(3,QQ)
2668 sage: J = DirectSumEJA(J1,J2)
2670 (Euclidean Jordan algebra of dimension 2 over Rational Field,
2671 Euclidean Jordan algebra of dimension 3 over Rational Field)
2674 return self
._factors
2676 def projections(self
):
2678 Return a pair of projections onto this algebra's factors.
2682 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
2683 ....: ComplexHermitianEJA,
2688 sage: J1 = JordanSpinEJA(2)
2689 sage: J2 = ComplexHermitianEJA(2)
2690 sage: J = DirectSumEJA(J1,J2)
2691 sage: (pi_left, pi_right) = J.projections()
2692 sage: J.one().to_vector()
2694 sage: pi_left(J.one()).to_vector()
2696 sage: pi_right(J.one()).to_vector()
2700 (J1
,J2
) = self
.factors()
2703 V_basis
= self
.vector_space().basis()
2704 # Need to specify the dimensions explicitly so that we don't
2705 # wind up with a zero-by-zero matrix when we want e.g. a
2706 # zero-by-two matrix (important for composing things).
2707 P1
= matrix(self
.base_ring(), m
, m
+n
, V_basis
[:m
])
2708 P2
= matrix(self
.base_ring(), n
, m
+n
, V_basis
[m
:])
2709 pi_left
= FiniteDimensionalEuclideanJordanAlgebraOperator(self
,J1
,P1
)
2710 pi_right
= FiniteDimensionalEuclideanJordanAlgebraOperator(self
,J2
,P2
)
2711 return (pi_left
, pi_right
)
2713 def inclusions(self
):
2715 Return the pair of inclusion maps from our factors into us.
2719 sage: from mjo.eja.eja_algebra import (random_eja,
2720 ....: JordanSpinEJA,
2721 ....: RealSymmetricEJA,
2726 sage: J1 = JordanSpinEJA(3)
2727 sage: J2 = RealSymmetricEJA(2)
2728 sage: J = DirectSumEJA(J1,J2)
2729 sage: (iota_left, iota_right) = J.inclusions()
2730 sage: iota_left(J1.zero()) == J.zero()
2732 sage: iota_right(J2.zero()) == J.zero()
2734 sage: J1.one().to_vector()
2736 sage: iota_left(J1.one()).to_vector()
2738 sage: J2.one().to_vector()
2740 sage: iota_right(J2.one()).to_vector()
2742 sage: J.one().to_vector()
2747 Composing a projection with the corresponding inclusion should
2748 produce the identity map, and mismatching them should produce
2751 sage: set_random_seed()
2752 sage: J1 = random_eja()
2753 sage: J2 = random_eja()
2754 sage: J = DirectSumEJA(J1,J2)
2755 sage: (iota_left, iota_right) = J.inclusions()
2756 sage: (pi_left, pi_right) = J.projections()
2757 sage: pi_left*iota_left == J1.one().operator()
2759 sage: pi_right*iota_right == J2.one().operator()
2761 sage: (pi_left*iota_right).is_zero()
2763 sage: (pi_right*iota_left).is_zero()
2767 (J1
,J2
) = self
.factors()
2770 V_basis
= self
.vector_space().basis()
2771 # Need to specify the dimensions explicitly so that we don't
2772 # wind up with a zero-by-zero matrix when we want e.g. a
2773 # two-by-zero matrix (important for composing things).
2774 I1
= matrix
.column(self
.base_ring(), m
, m
+n
, V_basis
[:m
])
2775 I2
= matrix
.column(self
.base_ring(), n
, m
+n
, V_basis
[m
:])
2776 iota_left
= FiniteDimensionalEuclideanJordanAlgebraOperator(J1
,self
,I1
)
2777 iota_right
= FiniteDimensionalEuclideanJordanAlgebraOperator(J2
,self
,I2
)
2778 return (iota_left
, iota_right
)
2780 def inner_product(self
, x
, y
):
2782 The standard Cartesian inner-product.
2784 We project ``x`` and ``y`` onto our factors, and add up the
2785 inner-products from the subalgebras.
2790 sage: from mjo.eja.eja_algebra import (HadamardEJA,
2791 ....: QuaternionHermitianEJA,
2796 sage: J1 = HadamardEJA(3,QQ)
2797 sage: J2 = QuaternionHermitianEJA(2,QQ,orthonormalize=False)
2798 sage: J = DirectSumEJA(J1,J2)
2803 sage: x1.inner_product(x2)
2805 sage: y1.inner_product(y2)
2807 sage: J.one().inner_product(J.one())
2811 (pi_left
, pi_right
) = self
.projections()
2817 return (x1
.inner_product(y1
) + x2
.inner_product(y2
))
2821 random_eja
= ConcreteEuclideanJordanAlgebra
.random_instance