2 Euclidean Jordan Algebras. These are formally-real Jordan Algebras;
3 specifically those where u^2 + v^2 = 0 implies that u = v = 0. They
4 are used in optimization, and have some additional nice methods beyond
5 what can be supported in a general Jordan Algebra.
10 sage: from mjo.eja.eja_algebra import random_eja
15 Euclidean Jordan algebra of dimension...
19 from itertools
import repeat
21 from sage
.algebras
.quatalg
.quaternion_algebra
import QuaternionAlgebra
22 from sage
.categories
.magmatic_algebras
import MagmaticAlgebras
23 from sage
.combinat
.free_module
import CombinatorialFreeModule
24 from sage
.matrix
.constructor
import matrix
25 from sage
.matrix
.matrix_space
import MatrixSpace
26 from sage
.misc
.cachefunc
import cached_method
27 from sage
.misc
.table
import table
28 from sage
.modules
.free_module
import FreeModule
, VectorSpace
29 from sage
.rings
.all
import (ZZ
, QQ
, AA
, QQbar
, RR
, RLF
, CLF
,
32 from mjo
.eja
.eja_element
import FiniteDimensionalEuclideanJordanAlgebraElement
33 from mjo
.eja
.eja_operator
import FiniteDimensionalEuclideanJordanAlgebraOperator
34 from mjo
.eja
.eja_utils
import _mat2vec
36 class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule
):
38 The lowest-level class for representing a Euclidean Jordan algebra.
40 def _coerce_map_from_base_ring(self
):
42 Disable the map from the base ring into the algebra.
44 Performing a nonsense conversion like this automatically
45 is counterpedagogical. The fallback is to try the usual
46 element constructor, which should also fail.
50 sage: from mjo.eja.eja_algebra import random_eja
54 sage: set_random_seed()
55 sage: J = random_eja()
57 Traceback (most recent call last):
59 ValueError: not an element of this algebra
76 * field -- the scalar field for this algebra (must be real)
78 * multiplication_table -- the multiplication table for this
79 algebra's implicit basis. Only the lower-triangular portion
80 of the table is used, since the multiplication is assumed
85 sage: from mjo.eja.eja_algebra import (
86 ....: FiniteDimensionalEuclideanJordanAlgebra,
92 By definition, Jordan multiplication commutes::
94 sage: set_random_seed()
95 sage: J = random_eja()
96 sage: x,y = J.random_elements(2)
100 An error is raised if the Jordan product is not commutative::
102 sage: JP = ((1,2),(0,0))
103 sage: IP = ((1,0),(0,1))
104 sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,JP,IP)
105 Traceback (most recent call last):
107 ValueError: Jordan product is not commutative
109 An error is raised if the inner-product is not commutative::
111 sage: JP = ((1,0),(0,1))
112 sage: IP = ((1,2),(0,0))
113 sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,JP,IP)
114 Traceback (most recent call last):
116 ValueError: inner-product is not commutative
120 The ``field`` we're given must be real with ``check_field=True``::
122 sage: JordanSpinEJA(2,QQbar)
123 Traceback (most recent call last):
125 ValueError: scalar field is not real
127 The multiplication table must be square with ``check_axioms=True``::
129 sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,((),()),((1,),))
130 Traceback (most recent call last):
132 ValueError: multiplication table is not square
134 The multiplication and inner-product tables must be the same
135 size (and in particular, the inner-product table must also be
136 square) with ``check_axioms=True``::
138 sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,((1,),),(()))
139 Traceback (most recent call last):
141 ValueError: multiplication and inner-product tables are
143 sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,((1,),),((1,2),))
144 Traceback (most recent call last):
146 ValueError: multiplication and inner-product tables are
151 if not field
.is_subring(RR
):
152 # Note: this does return true for the real algebraic
153 # field, the rationals, and any quadratic field where
154 # we've specified a real embedding.
155 raise ValueError("scalar field is not real")
158 # The multiplication and inner-product tables should be square
159 # if the user wants us to verify them. And we verify them as
160 # soon as possible, because we want to exploit their symmetry.
161 n
= len(multiplication_table
)
163 if not all( len(l
) == n
for l
in multiplication_table
):
164 raise ValueError("multiplication table is not square")
166 # If the multiplication table is square, we can check if
167 # the inner-product table is square by comparing it to the
168 # multiplication table's dimensions.
169 msg
= "multiplication and inner-product tables are different sizes"
170 if not len(inner_product_table
) == n
:
171 raise ValueError(msg
)
173 if not all( len(l
) == n
for l
in inner_product_table
):
174 raise ValueError(msg
)
176 if not all( multiplication_table
[j
][i
]
177 == multiplication_table
[i
][j
]
179 for j
in range(i
+1) ):
180 raise ValueError("Jordan product is not commutative")
181 if not all( inner_product_table
[j
][i
]
182 == inner_product_table
[i
][j
]
184 for j
in range(i
+1) ):
185 raise ValueError("inner-product is not commutative")
186 self
._matrix
_basis
= matrix_basis
189 category
= MagmaticAlgebras(field
).FiniteDimensional()
190 category
= category
.WithBasis().Unital()
192 fda
= super(FiniteDimensionalEuclideanJordanAlgebra
, self
)
197 self
.print_options(bracket
='')
199 # The multiplication table we're given is necessarily in terms
200 # of vectors, because we don't have an algebra yet for
201 # anything to be an element of. However, it's faster in the
202 # long run to have the multiplication table be in terms of
203 # algebra elements. We do this after calling the superclass
204 # constructor so that from_vector() knows what to do.
206 # Note: we take advantage of symmetry here, and only store
207 # the lower-triangular portion of the table.
208 self
._multiplication
_table
= [ [ self
.vector_space().zero()
209 for j
in range(i
+1) ]
214 elt
= self
.from_vector(multiplication_table
[i
][j
])
215 self
._multiplication
_table
[i
][j
] = elt
217 self
._multiplication
_table
= tuple(map(tuple, self
._multiplication
_table
))
219 # Save our inner product as a matrix, since the efficiency of
220 # matrix multiplication will usually outweigh the fact that we
221 # have to store a redundant upper- or lower-triangular part.
222 # Pre-cache the fact that these are Hermitian (real symmetric,
223 # in fact) in case some e.g. matrix multiplication routine can
224 # take advantage of it.
225 self
._inner
_product
_matrix
= matrix(field
, inner_product_table
)
226 self
._inner
_product
_matrix
._cache
= {'hermitian': True}
227 self
._inner
_product
_matrix
.set_immutable()
230 if not self
._is
_jordanian
():
231 raise ValueError("Jordan identity does not hold")
232 if not self
._inner
_product
_is
_associative
():
233 raise ValueError("inner product is not associative")
235 def _element_constructor_(self
, elt
):
237 Construct an element of this algebra from its vector or matrix
240 This gets called only after the parent element _call_ method
241 fails to find a coercion for the argument.
245 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
247 ....: RealSymmetricEJA)
251 The identity in `S^n` is converted to the identity in the EJA::
253 sage: J = RealSymmetricEJA(3)
254 sage: I = matrix.identity(QQ,3)
255 sage: J(I) == J.one()
258 This skew-symmetric matrix can't be represented in the EJA::
260 sage: J = RealSymmetricEJA(3)
261 sage: A = matrix(QQ,3, lambda i,j: i-j)
263 Traceback (most recent call last):
265 ValueError: not an element of this algebra
269 Ensure that we can convert any element of the two non-matrix
270 simple algebras (whose matrix representations are columns)
271 back and forth faithfully::
273 sage: set_random_seed()
274 sage: J = HadamardEJA.random_instance()
275 sage: x = J.random_element()
276 sage: J(x.to_vector().column()) == x
278 sage: J = JordanSpinEJA.random_instance()
279 sage: x = J.random_element()
280 sage: J(x.to_vector().column()) == x
283 msg
= "not an element of this algebra"
285 # The superclass implementation of random_element()
286 # needs to be able to coerce "0" into the algebra.
288 elif elt
in self
.base_ring():
289 # Ensure that no base ring -> algebra coercion is performed
290 # by this method. There's some stupidity in sage that would
291 # otherwise propagate to this method; for example, sage thinks
292 # that the integer 3 belongs to the space of 2-by-2 matrices.
293 raise ValueError(msg
)
295 if elt
not in self
.matrix_space():
296 raise ValueError(msg
)
298 # Thanks for nothing! Matrix spaces aren't vector spaces in
299 # Sage, so we have to figure out its matrix-basis coordinates
300 # ourselves. We use the basis space's ring instead of the
301 # element's ring because the basis space might be an algebraic
302 # closure whereas the base ring of the 3-by-3 identity matrix
303 # could be QQ instead of QQbar.
304 V
= VectorSpace(self
.base_ring(), elt
.nrows()*elt
.ncols())
305 W
= V
.span_of_basis( _mat2vec(s
) for s
in self
.matrix_basis() )
308 coords
= W
.coordinate_vector(_mat2vec(elt
))
309 except ArithmeticError: # vector is not in free module
310 raise ValueError(msg
)
312 return self
.from_vector(coords
)
316 Return a string representation of ``self``.
320 sage: from mjo.eja.eja_algebra import JordanSpinEJA
324 Ensure that it says what we think it says::
326 sage: JordanSpinEJA(2, field=AA)
327 Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
328 sage: JordanSpinEJA(3, field=RDF)
329 Euclidean Jordan algebra of dimension 3 over Real Double Field
332 fmt
= "Euclidean Jordan algebra of dimension {} over {}"
333 return fmt
.format(self
.dimension(), self
.base_ring())
335 def product_on_basis(self
, i
, j
):
336 # We only stored the lower-triangular portion of the
337 # multiplication table.
339 return self
._multiplication
_table
[i
][j
]
341 return self
._multiplication
_table
[j
][i
]
343 def _is_commutative(self
):
345 Whether or not this algebra's multiplication table is commutative.
347 This method should of course always return ``True``, unless
348 this algebra was constructed with ``check_axioms=False`` and
349 passed an invalid multiplication table.
351 return all( self
.product_on_basis(i
,j
) == self
.product_on_basis(i
,j
)
352 for i
in range(self
.dimension())
353 for j
in range(self
.dimension()) )
355 def _is_jordanian(self
):
357 Whether or not this algebra's multiplication table respects the
358 Jordan identity `(x^{2})(xy) = x(x^{2}y)`.
360 We only check one arrangement of `x` and `y`, so for a
361 ``True`` result to be truly true, you should also check
362 :meth:`_is_commutative`. This method should of course always
363 return ``True``, unless this algebra was constructed with
364 ``check_axioms=False`` and passed an invalid multiplication table.
366 return all( (self
.monomial(i
)**2)*(self
.monomial(i
)*self
.monomial(j
))
368 (self
.monomial(i
))*((self
.monomial(i
)**2)*self
.monomial(j
))
369 for i
in range(self
.dimension())
370 for j
in range(self
.dimension()) )
372 def _inner_product_is_associative(self
):
374 Return whether or not this algebra's inner product `B` is
375 associative; that is, whether or not `B(xy,z) = B(x,yz)`.
377 This method should of course always return ``True``, unless
378 this algebra was constructed with ``check_axioms=False`` and
379 passed an invalid multiplication table.
382 # Used to check whether or not something is zero in an inexact
383 # ring. This number is sufficient to allow the construction of
384 # QuaternionHermitianEJA(2, RDF) with check_axioms=True.
387 for i
in range(self
.dimension()):
388 for j
in range(self
.dimension()):
389 for k
in range(self
.dimension()):
393 diff
= (x
*y
).inner_product(z
) - x
.inner_product(y
*z
)
395 if self
.base_ring().is_exact():
399 if diff
.abs() > epsilon
:
405 def characteristic_polynomial_of(self
):
407 Return the algebra's "characteristic polynomial of" function,
408 which is itself a multivariate polynomial that, when evaluated
409 at the coordinates of some algebra element, returns that
410 element's characteristic polynomial.
412 The resulting polynomial has `n+1` variables, where `n` is the
413 dimension of this algebra. The first `n` variables correspond to
414 the coordinates of an algebra element: when evaluated at the
415 coordinates of an algebra element with respect to a certain
416 basis, the result is a univariate polynomial (in the one
417 remaining variable ``t``), namely the characteristic polynomial
422 sage: from mjo.eja.eja_algebra import JordanSpinEJA, TrivialEJA
426 The characteristic polynomial in the spin algebra is given in
427 Alizadeh, Example 11.11::
429 sage: J = JordanSpinEJA(3)
430 sage: p = J.characteristic_polynomial_of(); p
431 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
432 sage: xvec = J.one().to_vector()
436 By definition, the characteristic polynomial is a monic
437 degree-zero polynomial in a rank-zero algebra. Note that
438 Cayley-Hamilton is indeed satisfied since the polynomial
439 ``1`` evaluates to the identity element of the algebra on
442 sage: J = TrivialEJA()
443 sage: J.characteristic_polynomial_of()
450 # The list of coefficient polynomials a_0, a_1, a_2, ..., a_(r-1).
451 a
= self
._charpoly
_coefficients
()
453 # We go to a bit of trouble here to reorder the
454 # indeterminates, so that it's easier to evaluate the
455 # characteristic polynomial at x's coordinates and get back
456 # something in terms of t, which is what we want.
457 S
= PolynomialRing(self
.base_ring(),'t')
461 S
= PolynomialRing(S
, R
.variable_names())
464 return (t
**r
+ sum( a
[k
]*(t
**k
) for k
in range(r
) ))
466 def coordinate_polynomial_ring(self
):
468 The multivariate polynomial ring in which this algebra's
469 :meth:`characteristic_polynomial_of` lives.
473 sage: from mjo.eja.eja_algebra import (HadamardEJA,
474 ....: RealSymmetricEJA)
478 sage: J = HadamardEJA(2)
479 sage: J.coordinate_polynomial_ring()
480 Multivariate Polynomial Ring in X1, X2...
481 sage: J = RealSymmetricEJA(3,QQ)
482 sage: J.coordinate_polynomial_ring()
483 Multivariate Polynomial Ring in X1, X2, X3, X4, X5, X6...
486 var_names
= tuple( "X%d" % z
for z
in range(1, self
.dimension()+1) )
487 return PolynomialRing(self
.base_ring(), var_names
)
489 def inner_product(self
, x
, y
):
491 The inner product associated with this Euclidean Jordan algebra.
493 Defaults to the trace inner product, but can be overridden by
494 subclasses if they are sure that the necessary properties are
499 sage: from mjo.eja.eja_algebra import (random_eja,
501 ....: BilinearFormEJA)
505 Our inner product is "associative," which means the following for
506 a symmetric bilinear form::
508 sage: set_random_seed()
509 sage: J = random_eja()
510 sage: x,y,z = J.random_elements(3)
511 sage: (x*y).inner_product(z) == y.inner_product(x*z)
516 Ensure that this is the usual inner product for the algebras
519 sage: set_random_seed()
520 sage: J = HadamardEJA.random_instance()
521 sage: x,y = J.random_elements(2)
522 sage: actual = x.inner_product(y)
523 sage: expected = x.to_vector().inner_product(y.to_vector())
524 sage: actual == expected
527 Ensure that this is one-half of the trace inner-product in a
528 BilinearFormEJA that isn't just the reals (when ``n`` isn't
529 one). This is in Faraut and Koranyi, and also my "On the
532 sage: set_random_seed()
533 sage: J = BilinearFormEJA.random_instance()
534 sage: n = J.dimension()
535 sage: x = J.random_element()
536 sage: y = J.random_element()
537 sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
540 B
= self
._inner
_product
_matrix
541 return (B
*x
.to_vector()).inner_product(y
.to_vector())
544 def is_trivial(self
):
546 Return whether or not this algebra is trivial.
548 A trivial algebra contains only the zero element.
552 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
557 sage: J = ComplexHermitianEJA(3)
563 sage: J = TrivialEJA()
568 return self
.dimension() == 0
571 def multiplication_table(self
):
573 Return a visual representation of this algebra's multiplication
574 table (on basis elements).
578 sage: from mjo.eja.eja_algebra import JordanSpinEJA
582 sage: J = JordanSpinEJA(4)
583 sage: J.multiplication_table()
584 +----++----+----+----+----+
585 | * || e0 | e1 | e2 | e3 |
586 +====++====+====+====+====+
587 | e0 || e0 | e1 | e2 | e3 |
588 +----++----+----+----+----+
589 | e1 || e1 | e0 | 0 | 0 |
590 +----++----+----+----+----+
591 | e2 || e2 | 0 | e0 | 0 |
592 +----++----+----+----+----+
593 | e3 || e3 | 0 | 0 | e0 |
594 +----++----+----+----+----+
598 M
= [ [ self
.zero() for j
in range(n
) ]
602 M
[i
][j
] = self
._multiplication
_table
[i
][j
]
606 # Prepend the left "header" column entry Can't do this in
607 # the loop because it messes up the symmetry.
608 M
[i
] = [self
.monomial(i
)] + M
[i
]
610 # Prepend the header row.
611 M
= [["*"] + list(self
.gens())] + M
612 return table(M
, header_row
=True, header_column
=True, frame
=True)
615 def matrix_basis(self
):
617 Return an (often more natural) representation of this algebras
618 basis as an ordered tuple of matrices.
620 Every finite-dimensional Euclidean Jordan Algebra is a, up to
621 Jordan isomorphism, a direct sum of five simple
622 algebras---four of which comprise Hermitian matrices. And the
623 last type of algebra can of course be thought of as `n`-by-`1`
624 column matrices (ambiguusly called column vectors) to avoid
625 special cases. As a result, matrices (and column vectors) are
626 a natural representation format for Euclidean Jordan algebra
629 But, when we construct an algebra from a basis of matrices,
630 those matrix representations are lost in favor of coordinate
631 vectors *with respect to* that basis. We could eventually
632 convert back if we tried hard enough, but having the original
633 representations handy is valuable enough that we simply store
634 them and return them from this method.
636 Why implement this for non-matrix algebras? Avoiding special
637 cases for the :class:`BilinearFormEJA` pays with simplicity in
638 its own right. But mainly, we would like to be able to assume
639 that elements of a :class:`DirectSumEJA` can be displayed
640 nicely, without having to have special classes for direct sums
641 one of whose components was a matrix algebra.
645 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
646 ....: RealSymmetricEJA)
650 sage: J = RealSymmetricEJA(2)
652 Finite family {0: e0, 1: e1, 2: e2}
653 sage: J.matrix_basis()
655 [1 0] [ 0 0.7071067811865475?] [0 0]
656 [0 0], [0.7071067811865475? 0], [0 1]
661 sage: J = JordanSpinEJA(2)
663 Finite family {0: e0, 1: e1}
664 sage: J.matrix_basis()
670 if self
._matrix
_basis
is None:
671 M
= self
.matrix_space()
672 return tuple( M(b
.to_vector()) for b
in self
.basis() )
674 return self
._matrix
_basis
677 def matrix_space(self
):
679 Return the matrix space in which this algebra's elements live, if
680 we think of them as matrices (including column vectors of the
683 Generally this will be an `n`-by-`1` column-vector space,
684 except when the algebra is trivial. There it's `n`-by-`n`
685 (where `n` is zero), to ensure that two elements of the matrix
686 space (empty matrices) can be multiplied.
688 Matrix algebras override this with something more useful.
690 if self
.is_trivial():
691 return MatrixSpace(self
.base_ring(), 0)
692 elif self
._matrix
_basis
is None or len(self
._matrix
_basis
) == 0:
693 return MatrixSpace(self
.base_ring(), self
.dimension(), 1)
695 return self
._matrix
_basis
[0].matrix_space()
701 Return the unit element of this algebra.
705 sage: from mjo.eja.eja_algebra import (HadamardEJA,
710 sage: J = HadamardEJA(5)
712 e0 + e1 + e2 + e3 + e4
716 The identity element acts like the identity::
718 sage: set_random_seed()
719 sage: J = random_eja()
720 sage: x = J.random_element()
721 sage: J.one()*x == x and x*J.one() == x
724 The matrix of the unit element's operator is the identity::
726 sage: set_random_seed()
727 sage: J = random_eja()
728 sage: actual = J.one().operator().matrix()
729 sage: expected = matrix.identity(J.base_ring(), J.dimension())
730 sage: actual == expected
733 Ensure that the cached unit element (often precomputed by
734 hand) agrees with the computed one::
736 sage: set_random_seed()
737 sage: J = random_eja()
738 sage: cached = J.one()
739 sage: J.one.clear_cache()
740 sage: J.one() == cached
744 # We can brute-force compute the matrices of the operators
745 # that correspond to the basis elements of this algebra.
746 # If some linear combination of those basis elements is the
747 # algebra identity, then the same linear combination of
748 # their matrices has to be the identity matrix.
750 # Of course, matrices aren't vectors in sage, so we have to
751 # appeal to the "long vectors" isometry.
752 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
754 # Now we use basic linear algebra to find the coefficients,
755 # of the matrices-as-vectors-linear-combination, which should
756 # work for the original algebra basis too.
757 A
= matrix(self
.base_ring(), oper_vecs
)
759 # We used the isometry on the left-hand side already, but we
760 # still need to do it for the right-hand side. Recall that we
761 # wanted something that summed to the identity matrix.
762 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
764 # Now if there's an identity element in the algebra, this
765 # should work. We solve on the left to avoid having to
766 # transpose the matrix "A".
767 return self
.from_vector(A
.solve_left(b
))
770 def peirce_decomposition(self
, c
):
772 The Peirce decomposition of this algebra relative to the
775 In the future, this can be extended to a complete system of
776 orthogonal idempotents.
780 - ``c`` -- an idempotent of this algebra.
784 A triple (J0, J5, J1) containing two subalgebras and one subspace
787 - ``J0`` -- the algebra on the eigenspace of ``c.operator()``
788 corresponding to the eigenvalue zero.
790 - ``J5`` -- the eigenspace (NOT a subalgebra) of ``c.operator()``
791 corresponding to the eigenvalue one-half.
793 - ``J1`` -- the algebra on the eigenspace of ``c.operator()``
794 corresponding to the eigenvalue one.
796 These are the only possible eigenspaces for that operator, and this
797 algebra is a direct sum of them. The spaces ``J0`` and ``J1`` are
798 orthogonal, and are subalgebras of this algebra with the appropriate
803 sage: from mjo.eja.eja_algebra import random_eja, RealSymmetricEJA
807 The canonical example comes from the symmetric matrices, which
808 decompose into diagonal and off-diagonal parts::
810 sage: J = RealSymmetricEJA(3)
811 sage: C = matrix(QQ, [ [1,0,0],
815 sage: J0,J5,J1 = J.peirce_decomposition(c)
817 Euclidean Jordan algebra of dimension 1...
819 Vector space of degree 6 and dimension 2...
821 Euclidean Jordan algebra of dimension 3...
822 sage: J0.one().to_matrix()
826 sage: orig_df = AA.options.display_format
827 sage: AA.options.display_format = 'radical'
828 sage: J.from_vector(J5.basis()[0]).to_matrix()
832 sage: J.from_vector(J5.basis()[1]).to_matrix()
836 sage: AA.options.display_format = orig_df
837 sage: J1.one().to_matrix()
844 Every algebra decomposes trivially with respect to its identity
847 sage: set_random_seed()
848 sage: J = random_eja()
849 sage: J0,J5,J1 = J.peirce_decomposition(J.one())
850 sage: J0.dimension() == 0 and J5.dimension() == 0
852 sage: J1.superalgebra() == J and J1.dimension() == J.dimension()
855 The decomposition is into eigenspaces, and its components are
856 therefore necessarily orthogonal. Moreover, the identity
857 elements in the two subalgebras are the projections onto their
858 respective subspaces of the superalgebra's identity element::
860 sage: set_random_seed()
861 sage: J = random_eja()
862 sage: x = J.random_element()
863 sage: if not J.is_trivial():
864 ....: while x.is_nilpotent():
865 ....: x = J.random_element()
866 sage: c = x.subalgebra_idempotent()
867 sage: J0,J5,J1 = J.peirce_decomposition(c)
869 sage: for (w,y,z) in zip(J0.basis(), J5.basis(), J1.basis()):
870 ....: w = w.superalgebra_element()
871 ....: y = J.from_vector(y)
872 ....: z = z.superalgebra_element()
873 ....: ipsum += w.inner_product(y).abs()
874 ....: ipsum += w.inner_product(z).abs()
875 ....: ipsum += y.inner_product(z).abs()
878 sage: J1(c) == J1.one()
880 sage: J0(J.one() - c) == J0.one()
884 if not c
.is_idempotent():
885 raise ValueError("element is not idempotent: %s" % c
)
887 from mjo
.eja
.eja_subalgebra
import FiniteDimensionalEuclideanJordanSubalgebra
889 # Default these to what they should be if they turn out to be
890 # trivial, because eigenspaces_left() won't return eigenvalues
891 # corresponding to trivial spaces (e.g. it returns only the
892 # eigenspace corresponding to lambda=1 if you take the
893 # decomposition relative to the identity element).
894 trivial
= FiniteDimensionalEuclideanJordanSubalgebra(self
, ())
895 J0
= trivial
# eigenvalue zero
896 J5
= VectorSpace(self
.base_ring(), 0) # eigenvalue one-half
897 J1
= trivial
# eigenvalue one
899 for (eigval
, eigspace
) in c
.operator().matrix().right_eigenspaces():
900 if eigval
== ~
(self
.base_ring()(2)):
903 gens
= tuple( self
.from_vector(b
) for b
in eigspace
.basis() )
904 subalg
= FiniteDimensionalEuclideanJordanSubalgebra(self
,
912 raise ValueError("unexpected eigenvalue: %s" % eigval
)
917 def random_element(self
, thorough
=False):
919 Return a random element of this algebra.
921 Our algebra superclass method only returns a linear
922 combination of at most two basis elements. We instead
923 want the vector space "random element" method that
924 returns a more diverse selection.
928 - ``thorough`` -- (boolean; default False) whether or not we
929 should generate irrational coefficients for the random
930 element when our base ring is irrational; this slows the
931 algebra operations to a crawl, but any truly random method
935 # For a general base ring... maybe we can trust this to do the
936 # right thing? Unlikely, but.
937 V
= self
.vector_space()
938 v
= V
.random_element()
940 if self
.base_ring() is AA
:
941 # The "random element" method of the algebraic reals is
942 # stupid at the moment, and only returns integers between
943 # -2 and 2, inclusive:
945 # https://trac.sagemath.org/ticket/30875
947 # Instead, we implement our own "random vector" method,
948 # and then coerce that into the algebra. We use the vector
949 # space degree here instead of the dimension because a
950 # subalgebra could (for example) be spanned by only two
951 # vectors, each with five coordinates. We need to
952 # generate all five coordinates.
954 v
*= QQbar
.random_element().real()
956 v
*= QQ
.random_element()
958 return self
.from_vector(V
.coordinate_vector(v
))
960 def random_elements(self
, count
, thorough
=False):
962 Return ``count`` random elements as a tuple.
966 - ``thorough`` -- (boolean; default False) whether or not we
967 should generate irrational coefficients for the random
968 elements when our base ring is irrational; this slows the
969 algebra operations to a crawl, but any truly random method
974 sage: from mjo.eja.eja_algebra import JordanSpinEJA
978 sage: J = JordanSpinEJA(3)
979 sage: x,y,z = J.random_elements(3)
980 sage: all( [ x in J, y in J, z in J ])
982 sage: len( J.random_elements(10) ) == 10
986 return tuple( self
.random_element(thorough
)
987 for idx
in range(count
) )
991 def _charpoly_coefficients(self
):
993 The `r` polynomial coefficients of the "characteristic polynomial
997 R
= self
.coordinate_polynomial_ring()
999 F
= R
.fraction_field()
1002 # From a result in my book, these are the entries of the
1003 # basis representation of L_x.
1004 return sum( vars[k
]*self
.monomial(k
).operator().matrix()[i
,j
]
1007 L_x
= matrix(F
, n
, n
, L_x_i_j
)
1010 if self
.rank
.is_in_cache():
1012 # There's no need to pad the system with redundant
1013 # columns if we *know* they'll be redundant.
1016 # Compute an extra power in case the rank is equal to
1017 # the dimension (otherwise, we would stop at x^(r-1)).
1018 x_powers
= [ (L_x
**k
)*self
.one().to_vector()
1019 for k
in range(n
+1) ]
1020 A
= matrix
.column(F
, x_powers
[:n
])
1021 AE
= A
.extended_echelon_form()
1028 # The theory says that only the first "r" coefficients are
1029 # nonzero, and they actually live in the original polynomial
1030 # ring and not the fraction field. We negate them because
1031 # in the actual characteristic polynomial, they get moved
1032 # to the other side where x^r lives.
1033 return -A_rref
.solve_right(E
*b
).change_ring(R
)[:r
]
1038 Return the rank of this EJA.
1040 This is a cached method because we know the rank a priori for
1041 all of the algebras we can construct. Thus we can avoid the
1042 expensive ``_charpoly_coefficients()`` call unless we truly
1043 need to compute the whole characteristic polynomial.
1047 sage: from mjo.eja.eja_algebra import (HadamardEJA,
1048 ....: JordanSpinEJA,
1049 ....: RealSymmetricEJA,
1050 ....: ComplexHermitianEJA,
1051 ....: QuaternionHermitianEJA,
1056 The rank of the Jordan spin algebra is always two::
1058 sage: JordanSpinEJA(2).rank()
1060 sage: JordanSpinEJA(3).rank()
1062 sage: JordanSpinEJA(4).rank()
1065 The rank of the `n`-by-`n` Hermitian real, complex, or
1066 quaternion matrices is `n`::
1068 sage: RealSymmetricEJA(4).rank()
1070 sage: ComplexHermitianEJA(3).rank()
1072 sage: QuaternionHermitianEJA(2).rank()
1077 Ensure that every EJA that we know how to construct has a
1078 positive integer rank, unless the algebra is trivial in
1079 which case its rank will be zero::
1081 sage: set_random_seed()
1082 sage: J = random_eja()
1086 sage: r > 0 or (r == 0 and J.is_trivial())
1089 Ensure that computing the rank actually works, since the ranks
1090 of all simple algebras are known and will be cached by default::
1092 sage: set_random_seed() # long time
1093 sage: J = random_eja() # long time
1094 sage: caches = J.rank() # long time
1095 sage: J.rank.clear_cache() # long time
1096 sage: J.rank() == cached # long time
1100 return len(self
._charpoly
_coefficients
())
1103 def vector_space(self
):
1105 Return the vector space that underlies this algebra.
1109 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1113 sage: J = RealSymmetricEJA(2)
1114 sage: J.vector_space()
1115 Vector space of dimension 3 over...
1118 return self
.zero().to_vector().parent().ambient_vector_space()
1121 Element
= FiniteDimensionalEuclideanJordanAlgebraElement
1123 class RationalBasisEuclideanJordanAlgebraNg(FiniteDimensionalEuclideanJordanAlgebra
):
1125 New class for algebras whose supplied basis elements have all rational entries.
1129 sage: from mjo.eja.eja_algebra import BilinearFormEJA
1133 The supplied basis is orthonormalized by default::
1135 sage: B = matrix(QQ, [[1, 0, 0], [0, 25, -32], [0, -32, 41]])
1136 sage: J = BilinearFormEJA(B)
1137 sage: J.matrix_basis()
1150 orthonormalize
=True,
1157 vector_basis
= basis
1159 from sage
.structure
.element
import is_Matrix
1160 basis_is_matrices
= False
1164 if is_Matrix(basis
[0]):
1165 basis_is_matrices
= True
1166 vector_basis
= tuple( map(_mat2vec
,basis
) )
1167 degree
= basis
[0].nrows()**2
1169 degree
= basis
[0].degree()
1171 V
= VectorSpace(field
, degree
)
1173 # If we were asked to orthonormalize, and if the orthonormal
1174 # basis is different from the given one, then we also want to
1175 # compute multiplication and inner-product tables for the
1176 # deorthonormalized basis. These can be used later to
1177 # construct a deorthonormalized copy of this algebra over QQ
1178 # in which several operations are much faster.
1179 self
._deortho
_multiplication
_table
= None
1180 self
._deortho
_inner
_product
_table
= None
1183 # Compute the deorthonormalized tables before we orthonormalize
1185 W
= V
.span_of_basis( vector_basis
)
1187 # TODO: use symmetry
1188 self
._deortho
_multiplication
_table
= [ [0 for j
in range(n
)]
1190 self
._deortho
_inner
_product
_table
= [ [0 for j
in range(n
)]
1193 # Note: the Jordan and inner-products are defined in terms
1194 # of the ambient basis. It's important that their arguments
1195 # are in ambient coordinates as well.
1197 for j
in range(i
+1):
1198 # given basis w.r.t. ambient coords
1199 q_i
= vector_basis
[i
]
1200 q_j
= vector_basis
[j
]
1202 if basis_is_matrices
:
1206 elt
= jordan_product(q_i
, q_j
)
1207 ip
= inner_product(q_i
, q_j
)
1209 if basis_is_matrices
:
1210 # do another mat2vec because the multiplication
1211 # table is in terms of vectors
1214 # TODO: use symmetry
1215 elt
= W
.coordinate_vector(elt
)
1216 self
._deortho
_multiplication
_table
[i
][j
] = elt
1217 self
._deortho
_multiplication
_table
[j
][i
] = elt
1218 self
._deortho
_inner
_product
_table
[i
][j
] = ip
1219 self
._deortho
_inner
_product
_table
[j
][i
] = ip
1221 if self
._deortho
_multiplication
_table
is not None:
1222 self
._deortho
_multiplication
_table
= tuple(map(tuple, self
._deortho
_multiplication
_table
))
1223 if self
._deortho
_inner
_product
_table
is not None:
1224 self
._deortho
_inner
_product
_table
= tuple(map(tuple, self
._deortho
_inner
_product
_table
))
1227 from mjo
.eja
.eja_utils
import gram_schmidt
1228 vector_basis
= gram_schmidt(vector_basis
, inner_product
)
1229 W
= V
.span_of_basis( vector_basis
)
1231 # Normalize the "matrix" basis, too!
1232 basis
= vector_basis
1234 if basis_is_matrices
:
1235 from mjo
.eja
.eja_utils
import _vec2mat
1236 basis
= tuple( map(_vec2mat
,basis
) )
1238 W
= V
.span_of_basis( vector_basis
)
1240 # TODO: use symmetry
1241 mult_table
= [ [0 for j
in range(n
)] for i
in range(n
) ]
1242 ip_table
= [ [0 for j
in range(n
)] for i
in range(n
) ]
1244 # Note: the Jordan and inner- products are defined in terms
1245 # of the ambient basis. It's important that their arguments
1246 # are in ambient coordinates as well.
1248 for j
in range(i
+1):
1249 # ortho basis w.r.t. ambient coords
1250 q_i
= vector_basis
[i
]
1251 q_j
= vector_basis
[j
]
1253 if basis_is_matrices
:
1257 elt
= jordan_product(q_i
, q_j
)
1258 ip
= inner_product(q_i
, q_j
)
1260 if basis_is_matrices
:
1261 # do another mat2vec because the multiplication
1262 # table is in terms of vectors
1265 # TODO: use symmetry
1266 elt
= W
.coordinate_vector(elt
)
1267 mult_table
[i
][j
] = elt
1268 mult_table
[j
][i
] = elt
1272 if basis_is_matrices
:
1276 basis
= tuple( x
.column() for x
in basis
)
1278 super().__init
__(field
,
1283 basis
, # matrix basis
1288 def _charpoly_coefficients(self
):
1292 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
1293 ....: JordanSpinEJA)
1297 The returned coefficients should be the same as if we'd never
1298 orthonormalized the basis to begin with::
1300 sage: B = matrix(QQ, [[1, 0, 0],
1302 ....: [0, -32, 41] ])
1303 sage: J1 = BilinearFormEJA(B)
1304 sage: J2 = BilinearFormEJA(B,QQ,orthonormalize=False)
1305 sage: J1._charpoly_coefficients()
1306 (X1^2 - 25*X2^2 + 64*X2*X3 - 41*X3^2, -2*X1)
1307 sage: J2._charpoly_coefficients()
1308 (X1^2 - 25*X2^2 + 64*X2*X3 - 41*X3^2, -2*X1)
1310 The base ring of the resulting polynomial coefficients is what
1311 it should be, and not the rationals (unless the algebra was
1312 already over the rationals)::
1314 sage: J = JordanSpinEJA(3)
1315 sage: J._charpoly_coefficients()
1316 (X1^2 - X2^2 - X3^2, -2*X1)
1317 sage: a0 = J._charpoly_coefficients()[0]
1319 Algebraic Real Field
1320 sage: a0.base_ring()
1321 Algebraic Real Field
1323 if self
.base_ring() is QQ
:
1324 # There's no need to construct *another* algebra over the
1325 # rationals if this one is already over the rationals.
1326 superclass
= super(RationalBasisEuclideanJordanAlgebraNg
, self
)
1327 return superclass
._charpoly
_coefficients
()
1329 # Do the computation over the rationals. The answer will be
1330 # the same, because all we've done is a change of basis.
1331 J
= FiniteDimensionalEuclideanJordanAlgebra(QQ
,
1332 self
._deortho
_multiplication
_table
,
1333 self
._deortho
_inner
_product
_table
)
1334 a
= J
._charpoly
_coefficients
()
1335 return tuple(map(lambda x
: x
.change_ring(self
.base_ring()), a
))
1337 class ConcreteEuclideanJordanAlgebra
:
1339 A class for the Euclidean Jordan algebras that we know by name.
1341 These are the Jordan algebras whose basis, multiplication table,
1342 rank, and so on are known a priori. More to the point, they are
1343 the Euclidean Jordan algebras for which we are able to conjure up
1344 a "random instance."
1348 sage: from mjo.eja.eja_algebra import ConcreteEuclideanJordanAlgebra
1352 Our basis is normalized with respect to the algebra's inner
1353 product, unless we specify otherwise::
1355 sage: set_random_seed()
1356 sage: J = ConcreteEuclideanJordanAlgebra.random_instance()
1357 sage: all( b.norm() == 1 for b in J.gens() )
1360 Since our basis is orthonormal with respect to the algebra's inner
1361 product, and since we know that this algebra is an EJA, any
1362 left-multiplication operator's matrix will be symmetric because
1363 natural->EJA basis representation is an isometry and within the
1364 EJA the operator is self-adjoint by the Jordan axiom::
1366 sage: set_random_seed()
1367 sage: J = ConcreteEuclideanJordanAlgebra.random_instance()
1368 sage: x = J.random_element()
1369 sage: x.operator().is_self_adjoint()
1374 def _max_random_instance_size():
1376 Return an integer "size" that is an upper bound on the size of
1377 this algebra when it is used in a random test
1378 case. Unfortunately, the term "size" is ambiguous -- when
1379 dealing with `R^n` under either the Hadamard or Jordan spin
1380 product, the "size" refers to the dimension `n`. When dealing
1381 with a matrix algebra (real symmetric or complex/quaternion
1382 Hermitian), it refers to the size of the matrix, which is far
1383 less than the dimension of the underlying vector space.
1385 This method must be implemented in each subclass.
1387 raise NotImplementedError
1390 def random_instance(cls
, field
=AA
, **kwargs
):
1392 Return a random instance of this type of algebra.
1394 This method should be implemented in each subclass.
1396 from sage
.misc
.prandom
import choice
1397 eja_class
= choice(cls
.__subclasses
__())
1398 return eja_class
.random_instance(field
)
1401 class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra
):
1403 def __init__(self
, field
, basis
, normalize_basis
=True, **kwargs
):
1405 Compared to the superclass constructor, we take a basis instead of
1406 a multiplication table because the latter can be computed in terms
1407 of the former when the product is known (like it is here).
1409 # Used in this class's fast _charpoly_coefficients() override.
1410 self
._basis
_normalizers
= None
1412 # We're going to loop through this a few times, so now's a good
1413 # time to ensure that it isn't a generator expression.
1414 basis
= tuple(basis
)
1416 algebra_dim
= len(basis
)
1417 degree
= 0 # size of the matrices
1419 degree
= basis
[0].nrows()
1421 if algebra_dim
> 1 and normalize_basis
:
1422 # We'll need sqrt(2) to normalize the basis, and this
1423 # winds up in the multiplication table, so the whole
1424 # algebra needs to be over the field extension.
1425 R
= PolynomialRing(field
, 'z')
1428 if p
.is_irreducible():
1429 field
= field
.extension(p
, 'sqrt2', embedding
=RLF(2).sqrt())
1430 basis
= tuple( s
.change_ring(field
) for s
in basis
)
1431 self
._basis
_normalizers
= tuple(
1432 ~
(self
.matrix_inner_product(s
,s
).sqrt()) for s
in basis
)
1433 basis
= tuple(s
*c
for (s
,c
) in zip(basis
,self
._basis
_normalizers
))
1435 # Now compute the multiplication and inner product tables.
1436 # We have to do this *after* normalizing the basis, because
1437 # scaling affects the answers.
1438 V
= VectorSpace(field
, degree
**2)
1439 W
= V
.span_of_basis( _mat2vec(s
) for s
in basis
)
1440 mult_table
= [[W
.zero() for j
in range(algebra_dim
)]
1441 for i
in range(algebra_dim
)]
1442 ip_table
= [[field
.zero() for j
in range(algebra_dim
)]
1443 for i
in range(algebra_dim
)]
1444 for i
in range(algebra_dim
):
1445 for j
in range(algebra_dim
):
1446 mat_entry
= (basis
[i
]*basis
[j
] + basis
[j
]*basis
[i
])/2
1447 mult_table
[i
][j
] = W
.coordinate_vector(_mat2vec(mat_entry
))
1450 # HACK: ignore the error here if we don't need the
1451 # inner product (as is the case when we construct
1452 # a dummy QQ-algebra for fast charpoly coefficients.
1453 ip_table
[i
][j
] = self
.matrix_inner_product(basis
[i
],
1458 super(MatrixEuclideanJordanAlgebra
, self
).__init
__(field
,
1464 if algebra_dim
== 0:
1465 self
.one
.set_cache(self
.zero())
1467 n
= basis
[0].nrows()
1468 # The identity wrt (A,B) -> (AB + BA)/2 is independent of the
1469 # details of this algebra.
1470 self
.one
.set_cache(self(matrix
.identity(field
,n
)))
1474 def _charpoly_coefficients(self
):
1476 Override the parent method with something that tries to compute
1477 over a faster (non-extension) field.
1479 if self
._basis
_normalizers
is None or self
.base_ring() is QQ
:
1480 # We didn't normalize, or the basis we started with had
1481 # entries in a nice field already. Just compute the thing.
1482 return super(MatrixEuclideanJordanAlgebra
, self
)._charpoly
_coefficients
()
1484 basis
= ( (b
/n
) for (b
,n
) in zip(self
.matrix_basis(),
1485 self
._basis
_normalizers
) )
1487 # Do this over the rationals and convert back at the end.
1488 # Only works because we know the entries of the basis are
1489 # integers. The argument ``check_axioms=False`` is required
1490 # because the trace inner-product method for this
1491 # class is a stub and can't actually be checked.
1492 J
= MatrixEuclideanJordanAlgebra(QQ
,
1494 normalize_basis
=False,
1497 a
= J
._charpoly
_coefficients
()
1499 # Unfortunately, changing the basis does change the
1500 # coefficients of the characteristic polynomial, but since
1501 # these are really the coefficients of the "characteristic
1502 # polynomial of" function, everything is still nice and
1503 # unevaluated. It's therefore "obvious" how scaling the
1504 # basis affects the coordinate variables X1, X2, et
1505 # cetera. Scaling the first basis vector up by "n" adds a
1506 # factor of 1/n into every "X1" term, for example. So here
1507 # we simply undo the basis_normalizer scaling that we
1508 # performed earlier.
1510 # The a[0] access here is safe because trivial algebras
1511 # won't have any basis normalizers and therefore won't
1512 # make it to this "else" branch.
1513 XS
= a
[0].parent().gens()
1514 subs_dict
= { XS
[i
]: self
._basis
_normalizers
[i
]*XS
[i
]
1515 for i
in range(len(XS
)) }
1516 return tuple( a_i
.subs(subs_dict
) for a_i
in a
)
1522 Embed the matrix ``M`` into a space of real matrices.
1524 The matrix ``M`` can have entries in any field at the moment:
1525 the real numbers, complex numbers, or quaternions. And although
1526 they are not a field, we can probably support octonions at some
1527 point, too. This function returns a real matrix that "acts like"
1528 the original with respect to matrix multiplication; i.e.
1530 real_embed(M*N) = real_embed(M)*real_embed(N)
1533 raise NotImplementedError
1537 def real_unembed(M
):
1539 The inverse of :meth:`real_embed`.
1541 raise NotImplementedError
1544 def matrix_inner_product(cls
,X
,Y
):
1545 Xu
= cls
.real_unembed(X
)
1546 Yu
= cls
.real_unembed(Y
)
1547 tr
= (Xu
*Yu
).trace()
1550 # Works in QQ, AA, RDF, et cetera.
1552 except AttributeError:
1553 # A quaternion doesn't have a real() method, but does
1554 # have coefficient_tuple() method that returns the
1555 # coefficients of 1, i, j, and k -- in that order.
1556 return tr
.coefficient_tuple()[0]
1559 class RealMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1563 The identity function, for embedding real matrices into real
1569 def real_unembed(M
):
1571 The identity function, for unembedding real matrices from real
1577 class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra
,
1578 ConcreteEuclideanJordanAlgebra
):
1580 The rank-n simple EJA consisting of real symmetric n-by-n
1581 matrices, the usual symmetric Jordan product, and the trace inner
1582 product. It has dimension `(n^2 + n)/2` over the reals.
1586 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1590 sage: J = RealSymmetricEJA(2)
1591 sage: e0, e1, e2 = J.gens()
1599 In theory, our "field" can be any subfield of the reals::
1601 sage: RealSymmetricEJA(2, RDF)
1602 Euclidean Jordan algebra of dimension 3 over Real Double Field
1603 sage: RealSymmetricEJA(2, RR)
1604 Euclidean Jordan algebra of dimension 3 over Real Field with
1605 53 bits of precision
1609 The dimension of this algebra is `(n^2 + n) / 2`::
1611 sage: set_random_seed()
1612 sage: n_max = RealSymmetricEJA._max_random_instance_size()
1613 sage: n = ZZ.random_element(1, n_max)
1614 sage: J = RealSymmetricEJA(n)
1615 sage: J.dimension() == (n^2 + n)/2
1618 The Jordan multiplication is what we think it is::
1620 sage: set_random_seed()
1621 sage: J = RealSymmetricEJA.random_instance()
1622 sage: x,y = J.random_elements(2)
1623 sage: actual = (x*y).to_matrix()
1624 sage: X = x.to_matrix()
1625 sage: Y = y.to_matrix()
1626 sage: expected = (X*Y + Y*X)/2
1627 sage: actual == expected
1629 sage: J(expected) == x*y
1632 We can change the generator prefix::
1634 sage: RealSymmetricEJA(3, prefix='q').gens()
1635 (q0, q1, q2, q3, q4, q5)
1637 We can construct the (trivial) algebra of rank zero::
1639 sage: RealSymmetricEJA(0)
1640 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1644 def _denormalized_basis(cls
, n
, field
):
1646 Return a basis for the space of real symmetric n-by-n matrices.
1650 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1654 sage: set_random_seed()
1655 sage: n = ZZ.random_element(1,5)
1656 sage: B = RealSymmetricEJA._denormalized_basis(n,QQ)
1657 sage: all( M.is_symmetric() for M in B)
1661 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1665 for j
in range(i
+1):
1666 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
1670 Sij
= Eij
+ Eij
.transpose()
1676 def _max_random_instance_size():
1677 return 4 # Dimension 10
1680 def random_instance(cls
, field
=AA
, **kwargs
):
1682 Return a random instance of this type of algebra.
1684 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1685 return cls(n
, field
, **kwargs
)
1687 def __init__(self
, n
, field
=AA
, **kwargs
):
1688 basis
= self
._denormalized
_basis
(n
, field
)
1689 super(RealSymmetricEJA
, self
).__init
__(field
,
1693 self
.rank
.set_cache(n
)
1696 class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1700 Embed the n-by-n complex matrix ``M`` into the space of real
1701 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1702 bi` to the block matrix ``[[a,b],[-b,a]]``.
1706 sage: from mjo.eja.eja_algebra import \
1707 ....: ComplexMatrixEuclideanJordanAlgebra
1711 sage: F = QuadraticField(-1, 'I')
1712 sage: x1 = F(4 - 2*i)
1713 sage: x2 = F(1 + 2*i)
1716 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1717 sage: ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1726 Embedding is a homomorphism (isomorphism, in fact)::
1728 sage: set_random_seed()
1729 sage: n = ZZ.random_element(3)
1730 sage: F = QuadraticField(-1, 'I')
1731 sage: X = random_matrix(F, n)
1732 sage: Y = random_matrix(F, n)
1733 sage: Xe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X)
1734 sage: Ye = ComplexMatrixEuclideanJordanAlgebra.real_embed(Y)
1735 sage: XYe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1742 raise ValueError("the matrix 'M' must be square")
1744 # We don't need any adjoined elements...
1745 field
= M
.base_ring().base_ring()
1749 a
= z
.list()[0] # real part, I guess
1750 b
= z
.list()[1] # imag part, I guess
1751 blocks
.append(matrix(field
, 2, [[a
,b
],[-b
,a
]]))
1753 return matrix
.block(field
, n
, blocks
)
1757 def real_unembed(M
):
1759 The inverse of _embed_complex_matrix().
1763 sage: from mjo.eja.eja_algebra import \
1764 ....: ComplexMatrixEuclideanJordanAlgebra
1768 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1769 ....: [-2, 1, -4, 3],
1770 ....: [ 9, 10, 11, 12],
1771 ....: [-10, 9, -12, 11] ])
1772 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(A)
1774 [ 10*I + 9 12*I + 11]
1778 Unembedding is the inverse of embedding::
1780 sage: set_random_seed()
1781 sage: F = QuadraticField(-1, 'I')
1782 sage: M = random_matrix(F, 3)
1783 sage: Me = ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1784 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1790 raise ValueError("the matrix 'M' must be square")
1791 if not n
.mod(2).is_zero():
1792 raise ValueError("the matrix 'M' must be a complex embedding")
1794 # If "M" was normalized, its base ring might have roots
1795 # adjoined and they can stick around after unembedding.
1796 field
= M
.base_ring()
1797 R
= PolynomialRing(field
, 'z')
1800 # Sage doesn't know how to embed AA into QQbar, i.e. how
1801 # to adjoin sqrt(-1) to AA.
1804 F
= field
.extension(z
**2 + 1, 'I', embedding
=CLF(-1).sqrt())
1807 # Go top-left to bottom-right (reading order), converting every
1808 # 2-by-2 block we see to a single complex element.
1810 for k
in range(n
/2):
1811 for j
in range(n
/2):
1812 submat
= M
[2*k
:2*k
+2,2*j
:2*j
+2]
1813 if submat
[0,0] != submat
[1,1]:
1814 raise ValueError('bad on-diagonal submatrix')
1815 if submat
[0,1] != -submat
[1,0]:
1816 raise ValueError('bad off-diagonal submatrix')
1817 z
= submat
[0,0] + submat
[0,1]*i
1820 return matrix(F
, n
/2, elements
)
1824 def matrix_inner_product(cls
,X
,Y
):
1826 Compute a matrix inner product in this algebra directly from
1831 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1835 This gives the same answer as the slow, default method implemented
1836 in :class:`MatrixEuclideanJordanAlgebra`::
1838 sage: set_random_seed()
1839 sage: J = ComplexHermitianEJA.random_instance()
1840 sage: x,y = J.random_elements(2)
1841 sage: Xe = x.to_matrix()
1842 sage: Ye = y.to_matrix()
1843 sage: X = ComplexHermitianEJA.real_unembed(Xe)
1844 sage: Y = ComplexHermitianEJA.real_unembed(Ye)
1845 sage: expected = (X*Y).trace().real()
1846 sage: actual = ComplexHermitianEJA.matrix_inner_product(Xe,Ye)
1847 sage: actual == expected
1851 return RealMatrixEuclideanJordanAlgebra
.matrix_inner_product(X
,Y
)/2
1854 class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra
,
1855 ConcreteEuclideanJordanAlgebra
):
1857 The rank-n simple EJA consisting of complex Hermitian n-by-n
1858 matrices over the real numbers, the usual symmetric Jordan product,
1859 and the real-part-of-trace inner product. It has dimension `n^2` over
1864 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1868 In theory, our "field" can be any subfield of the reals::
1870 sage: ComplexHermitianEJA(2, RDF)
1871 Euclidean Jordan algebra of dimension 4 over Real Double Field
1872 sage: ComplexHermitianEJA(2, RR)
1873 Euclidean Jordan algebra of dimension 4 over Real Field with
1874 53 bits of precision
1878 The dimension of this algebra is `n^2`::
1880 sage: set_random_seed()
1881 sage: n_max = ComplexHermitianEJA._max_random_instance_size()
1882 sage: n = ZZ.random_element(1, n_max)
1883 sage: J = ComplexHermitianEJA(n)
1884 sage: J.dimension() == n^2
1887 The Jordan multiplication is what we think it is::
1889 sage: set_random_seed()
1890 sage: J = ComplexHermitianEJA.random_instance()
1891 sage: x,y = J.random_elements(2)
1892 sage: actual = (x*y).to_matrix()
1893 sage: X = x.to_matrix()
1894 sage: Y = y.to_matrix()
1895 sage: expected = (X*Y + Y*X)/2
1896 sage: actual == expected
1898 sage: J(expected) == x*y
1901 We can change the generator prefix::
1903 sage: ComplexHermitianEJA(2, prefix='z').gens()
1906 We can construct the (trivial) algebra of rank zero::
1908 sage: ComplexHermitianEJA(0)
1909 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1914 def _denormalized_basis(cls
, n
, field
):
1916 Returns a basis for the space of complex Hermitian n-by-n matrices.
1918 Why do we embed these? Basically, because all of numerical linear
1919 algebra assumes that you're working with vectors consisting of `n`
1920 entries from a field and scalars from the same field. There's no way
1921 to tell SageMath that (for example) the vectors contain complex
1922 numbers, while the scalar field is real.
1926 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1930 sage: set_random_seed()
1931 sage: n = ZZ.random_element(1,5)
1932 sage: field = QuadraticField(2, 'sqrt2')
1933 sage: B = ComplexHermitianEJA._denormalized_basis(n, field)
1934 sage: all( M.is_symmetric() for M in B)
1938 R
= PolynomialRing(field
, 'z')
1940 F
= field
.extension(z
**2 + 1, 'I')
1943 # This is like the symmetric case, but we need to be careful:
1945 # * We want conjugate-symmetry, not just symmetry.
1946 # * The diagonal will (as a result) be real.
1950 for j
in range(i
+1):
1951 Eij
= matrix(F
, n
, lambda k
,l
: k
==i
and l
==j
)
1953 Sij
= cls
.real_embed(Eij
)
1956 # The second one has a minus because it's conjugated.
1957 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
1959 Sij_imag
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
1962 # Since we embedded these, we can drop back to the "field" that we
1963 # started with instead of the complex extension "F".
1964 return ( s
.change_ring(field
) for s
in S
)
1967 def __init__(self
, n
, field
=AA
, **kwargs
):
1968 basis
= self
._denormalized
_basis
(n
,field
)
1969 super(ComplexHermitianEJA
,self
).__init
__(field
,
1973 self
.rank
.set_cache(n
)
1976 def _max_random_instance_size():
1977 return 3 # Dimension 9
1980 def random_instance(cls
, field
=AA
, **kwargs
):
1982 Return a random instance of this type of algebra.
1984 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1985 return cls(n
, field
, **kwargs
)
1987 class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1991 Embed the n-by-n quaternion matrix ``M`` into the space of real
1992 matrices of size 4n-by-4n by first sending each quaternion entry `z
1993 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
1994 c+di],[-c + di, a-bi]]`, and then embedding those into a real
1999 sage: from mjo.eja.eja_algebra import \
2000 ....: QuaternionMatrixEuclideanJordanAlgebra
2004 sage: Q = QuaternionAlgebra(QQ,-1,-1)
2005 sage: i,j,k = Q.gens()
2006 sage: x = 1 + 2*i + 3*j + 4*k
2007 sage: M = matrix(Q, 1, [[x]])
2008 sage: QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
2014 Embedding is a homomorphism (isomorphism, in fact)::
2016 sage: set_random_seed()
2017 sage: n = ZZ.random_element(2)
2018 sage: Q = QuaternionAlgebra(QQ,-1,-1)
2019 sage: X = random_matrix(Q, n)
2020 sage: Y = random_matrix(Q, n)
2021 sage: Xe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X)
2022 sage: Ye = QuaternionMatrixEuclideanJordanAlgebra.real_embed(Y)
2023 sage: XYe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X*Y)
2028 quaternions
= M
.base_ring()
2031 raise ValueError("the matrix 'M' must be square")
2033 F
= QuadraticField(-1, 'I')
2038 t
= z
.coefficient_tuple()
2043 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
2044 [-c
+ d
*i
, a
- b
*i
]])
2045 realM
= ComplexMatrixEuclideanJordanAlgebra
.real_embed(cplxM
)
2046 blocks
.append(realM
)
2048 # We should have real entries by now, so use the realest field
2049 # we've got for the return value.
2050 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
2055 def real_unembed(M
):
2057 The inverse of _embed_quaternion_matrix().
2061 sage: from mjo.eja.eja_algebra import \
2062 ....: QuaternionMatrixEuclideanJordanAlgebra
2066 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
2067 ....: [-2, 1, -4, 3],
2068 ....: [-3, 4, 1, -2],
2069 ....: [-4, -3, 2, 1]])
2070 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(M)
2071 [1 + 2*i + 3*j + 4*k]
2075 Unembedding is the inverse of embedding::
2077 sage: set_random_seed()
2078 sage: Q = QuaternionAlgebra(QQ, -1, -1)
2079 sage: M = random_matrix(Q, 3)
2080 sage: Me = QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
2081 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
2087 raise ValueError("the matrix 'M' must be square")
2088 if not n
.mod(4).is_zero():
2089 raise ValueError("the matrix 'M' must be a quaternion embedding")
2091 # Use the base ring of the matrix to ensure that its entries can be
2092 # multiplied by elements of the quaternion algebra.
2093 field
= M
.base_ring()
2094 Q
= QuaternionAlgebra(field
,-1,-1)
2097 # Go top-left to bottom-right (reading order), converting every
2098 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
2101 for l
in range(n
/4):
2102 for m
in range(n
/4):
2103 submat
= ComplexMatrixEuclideanJordanAlgebra
.real_unembed(
2104 M
[4*l
:4*l
+4,4*m
:4*m
+4] )
2105 if submat
[0,0] != submat
[1,1].conjugate():
2106 raise ValueError('bad on-diagonal submatrix')
2107 if submat
[0,1] != -submat
[1,0].conjugate():
2108 raise ValueError('bad off-diagonal submatrix')
2109 z
= submat
[0,0].real()
2110 z
+= submat
[0,0].imag()*i
2111 z
+= submat
[0,1].real()*j
2112 z
+= submat
[0,1].imag()*k
2115 return matrix(Q
, n
/4, elements
)
2119 def matrix_inner_product(cls
,X
,Y
):
2121 Compute a matrix inner product in this algebra directly from
2126 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2130 This gives the same answer as the slow, default method implemented
2131 in :class:`MatrixEuclideanJordanAlgebra`::
2133 sage: set_random_seed()
2134 sage: J = QuaternionHermitianEJA.random_instance()
2135 sage: x,y = J.random_elements(2)
2136 sage: Xe = x.to_matrix()
2137 sage: Ye = y.to_matrix()
2138 sage: X = QuaternionHermitianEJA.real_unembed(Xe)
2139 sage: Y = QuaternionHermitianEJA.real_unembed(Ye)
2140 sage: expected = (X*Y).trace().coefficient_tuple()[0]
2141 sage: actual = QuaternionHermitianEJA.matrix_inner_product(Xe,Ye)
2142 sage: actual == expected
2146 return RealMatrixEuclideanJordanAlgebra
.matrix_inner_product(X
,Y
)/4
2149 class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra
,
2150 ConcreteEuclideanJordanAlgebra
):
2152 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
2153 matrices, the usual symmetric Jordan product, and the
2154 real-part-of-trace inner product. It has dimension `2n^2 - n` over
2159 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2163 In theory, our "field" can be any subfield of the reals::
2165 sage: QuaternionHermitianEJA(2, RDF)
2166 Euclidean Jordan algebra of dimension 6 over Real Double Field
2167 sage: QuaternionHermitianEJA(2, RR)
2168 Euclidean Jordan algebra of dimension 6 over Real Field with
2169 53 bits of precision
2173 The dimension of this algebra is `2*n^2 - n`::
2175 sage: set_random_seed()
2176 sage: n_max = QuaternionHermitianEJA._max_random_instance_size()
2177 sage: n = ZZ.random_element(1, n_max)
2178 sage: J = QuaternionHermitianEJA(n)
2179 sage: J.dimension() == 2*(n^2) - n
2182 The Jordan multiplication is what we think it is::
2184 sage: set_random_seed()
2185 sage: J = QuaternionHermitianEJA.random_instance()
2186 sage: x,y = J.random_elements(2)
2187 sage: actual = (x*y).to_matrix()
2188 sage: X = x.to_matrix()
2189 sage: Y = y.to_matrix()
2190 sage: expected = (X*Y + Y*X)/2
2191 sage: actual == expected
2193 sage: J(expected) == x*y
2196 We can change the generator prefix::
2198 sage: QuaternionHermitianEJA(2, prefix='a').gens()
2199 (a0, a1, a2, a3, a4, a5)
2201 We can construct the (trivial) algebra of rank zero::
2203 sage: QuaternionHermitianEJA(0)
2204 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2208 def _denormalized_basis(cls
, n
, field
):
2210 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
2212 Why do we embed these? Basically, because all of numerical
2213 linear algebra assumes that you're working with vectors consisting
2214 of `n` entries from a field and scalars from the same field. There's
2215 no way to tell SageMath that (for example) the vectors contain
2216 complex numbers, while the scalar field is real.
2220 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2224 sage: set_random_seed()
2225 sage: n = ZZ.random_element(1,5)
2226 sage: B = QuaternionHermitianEJA._denormalized_basis(n,QQ)
2227 sage: all( M.is_symmetric() for M in B )
2231 Q
= QuaternionAlgebra(QQ
,-1,-1)
2234 # This is like the symmetric case, but we need to be careful:
2236 # * We want conjugate-symmetry, not just symmetry.
2237 # * The diagonal will (as a result) be real.
2241 for j
in range(i
+1):
2242 Eij
= matrix(Q
, n
, lambda k
,l
: k
==i
and l
==j
)
2244 Sij
= cls
.real_embed(Eij
)
2247 # The second, third, and fourth ones have a minus
2248 # because they're conjugated.
2249 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
2251 Sij_I
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
2253 Sij_J
= cls
.real_embed(J
*Eij
- J
*Eij
.transpose())
2255 Sij_K
= cls
.real_embed(K
*Eij
- K
*Eij
.transpose())
2258 # Since we embedded these, we can drop back to the "field" that we
2259 # started with instead of the quaternion algebra "Q".
2260 return ( s
.change_ring(field
) for s
in S
)
2263 def __init__(self
, n
, field
=AA
, **kwargs
):
2264 basis
= self
._denormalized
_basis
(n
,field
)
2265 super(QuaternionHermitianEJA
,self
).__init
__(field
,
2269 self
.rank
.set_cache(n
)
2272 def _max_random_instance_size():
2274 The maximum rank of a random QuaternionHermitianEJA.
2276 return 2 # Dimension 6
2279 def random_instance(cls
, field
=AA
, **kwargs
):
2281 Return a random instance of this type of algebra.
2283 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2284 return cls(n
, field
, **kwargs
)
2287 class HadamardEJA(RationalBasisEuclideanJordanAlgebraNg
,
2288 ConcreteEuclideanJordanAlgebra
):
2290 Return the Euclidean Jordan Algebra corresponding to the set
2291 `R^n` under the Hadamard product.
2293 Note: this is nothing more than the Cartesian product of ``n``
2294 copies of the spin algebra. Once Cartesian product algebras
2295 are implemented, this can go.
2299 sage: from mjo.eja.eja_algebra import HadamardEJA
2303 This multiplication table can be verified by hand::
2305 sage: J = HadamardEJA(3)
2306 sage: e0,e1,e2 = J.gens()
2322 We can change the generator prefix::
2324 sage: HadamardEJA(3, prefix='r').gens()
2328 def __init__(self
, n
, field
=AA
, **kwargs
):
2329 V
= VectorSpace(field
, n
)
2332 def jordan_product(x
,y
):
2333 return V([ xi
*yi
for (xi
,yi
) in zip(x
,y
) ])
2334 def inner_product(x
,y
):
2335 return x
.inner_product(y
)
2337 super(HadamardEJA
, self
).__init
__(field
,
2342 self
.rank
.set_cache(n
)
2345 self
.one
.set_cache( self
.zero() )
2347 self
.one
.set_cache( sum(self
.gens()) )
2350 def _max_random_instance_size():
2352 The maximum dimension of a random HadamardEJA.
2357 def random_instance(cls
, field
=AA
, **kwargs
):
2359 Return a random instance of this type of algebra.
2361 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2362 return cls(n
, field
, **kwargs
)
2365 class BilinearFormEJA(RationalBasisEuclideanJordanAlgebraNg
,
2366 ConcreteEuclideanJordanAlgebra
):
2368 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2369 with the half-trace inner product and jordan product ``x*y =
2370 (<Bx,y>,y_bar>, x0*y_bar + y0*x_bar)`` where `B = 1 \times B22` is
2371 a symmetric positive-definite "bilinear form" matrix. Its
2372 dimension is the size of `B`, and it has rank two in dimensions
2373 larger than two. It reduces to the ``JordanSpinEJA`` when `B` is
2374 the identity matrix of order ``n``.
2376 We insist that the one-by-one upper-left identity block of `B` be
2377 passed in as well so that we can be passed a matrix of size zero
2378 to construct a trivial algebra.
2382 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
2383 ....: JordanSpinEJA)
2387 When no bilinear form is specified, the identity matrix is used,
2388 and the resulting algebra is the Jordan spin algebra::
2390 sage: B = matrix.identity(AA,3)
2391 sage: J0 = BilinearFormEJA(B)
2392 sage: J1 = JordanSpinEJA(3)
2393 sage: J0.multiplication_table() == J0.multiplication_table()
2396 An error is raised if the matrix `B` does not correspond to a
2397 positive-definite bilinear form::
2399 sage: B = matrix.random(QQ,2,3)
2400 sage: J = BilinearFormEJA(B)
2401 Traceback (most recent call last):
2403 ValueError: bilinear form is not positive-definite
2404 sage: B = matrix.zero(QQ,3)
2405 sage: J = BilinearFormEJA(B)
2406 Traceback (most recent call last):
2408 ValueError: bilinear form is not positive-definite
2412 We can create a zero-dimensional algebra::
2414 sage: B = matrix.identity(AA,0)
2415 sage: J = BilinearFormEJA(B)
2419 We can check the multiplication condition given in the Jordan, von
2420 Neumann, and Wigner paper (and also discussed on my "On the
2421 symmetry..." paper). Note that this relies heavily on the standard
2422 choice of basis, as does anything utilizing the bilinear form
2423 matrix. We opt not to orthonormalize the basis, because if we
2424 did, we would have to normalize the `s_{i}` in a similar manner::
2426 sage: set_random_seed()
2427 sage: n = ZZ.random_element(5)
2428 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2429 sage: B11 = matrix.identity(QQ,1)
2430 sage: B22 = M.transpose()*M
2431 sage: B = block_matrix(2,2,[ [B11,0 ],
2433 sage: J = BilinearFormEJA(B, orthonormalize=False)
2434 sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
2435 sage: V = J.vector_space()
2436 sage: sis = [ J( V([0] + (M.inverse()*ei).list()).column() )
2437 ....: for ei in eis ]
2438 sage: actual = [ sis[i]*sis[j]
2439 ....: for i in range(n-1)
2440 ....: for j in range(n-1) ]
2441 sage: expected = [ J.one() if i == j else J.zero()
2442 ....: for i in range(n-1)
2443 ....: for j in range(n-1) ]
2444 sage: actual == expected
2447 def __init__(self
, B
, field
=AA
, **kwargs
):
2448 if not B
.is_positive_definite():
2449 raise ValueError("bilinear form is not positive-definite")
2452 V
= VectorSpace(field
, n
)
2454 def inner_product(x
,y
):
2455 return (B
*x
).inner_product(y
)
2457 def jordan_product(x
,y
):
2462 z0
= inner_product(x
,y
)
2463 zbar
= y0
*xbar
+ x0
*ybar
2464 return V([z0
] + zbar
.list())
2466 super(BilinearFormEJA
, self
).__init
__(field
,
2472 # The rank of this algebra is two, unless we're in a
2473 # one-dimensional ambient space (because the rank is bounded
2474 # by the ambient dimension).
2475 self
.rank
.set_cache(min(n
,2))
2478 self
.one
.set_cache( self
.zero() )
2480 self
.one
.set_cache( self
.monomial(0) )
2483 def _max_random_instance_size():
2485 The maximum dimension of a random BilinearFormEJA.
2490 def random_instance(cls
, field
=AA
, **kwargs
):
2492 Return a random instance of this algebra.
2494 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2496 B
= matrix
.identity(field
, n
)
2497 return cls(B
, field
, **kwargs
)
2499 B11
= matrix
.identity(field
,1)
2500 M
= matrix
.random(field
, n
-1)
2501 I
= matrix
.identity(field
, n
-1)
2502 alpha
= field
.zero()
2503 while alpha
.is_zero():
2504 alpha
= field
.random_element().abs()
2505 B22
= M
.transpose()*M
+ alpha
*I
2507 from sage
.matrix
.special
import block_matrix
2508 B
= block_matrix(2,2, [ [B11
, ZZ(0) ],
2511 return cls(B
, field
, **kwargs
)
2514 class JordanSpinEJA(BilinearFormEJA
):
2516 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2517 with the usual inner product and jordan product ``x*y =
2518 (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2523 sage: from mjo.eja.eja_algebra import JordanSpinEJA
2527 This multiplication table can be verified by hand::
2529 sage: J = JordanSpinEJA(4)
2530 sage: e0,e1,e2,e3 = J.gens()
2546 We can change the generator prefix::
2548 sage: JordanSpinEJA(2, prefix='B').gens()
2553 Ensure that we have the usual inner product on `R^n`::
2555 sage: set_random_seed()
2556 sage: J = JordanSpinEJA.random_instance()
2557 sage: x,y = J.random_elements(2)
2558 sage: actual = x.inner_product(y)
2559 sage: expected = x.to_vector().inner_product(y.to_vector())
2560 sage: actual == expected
2564 def __init__(self
, n
, field
=AA
, **kwargs
):
2565 # This is a special case of the BilinearFormEJA with the identity
2566 # matrix as its bilinear form.
2567 B
= matrix
.identity(field
, n
)
2568 super(JordanSpinEJA
, self
).__init
__(B
, field
, **kwargs
)
2571 def _max_random_instance_size():
2573 The maximum dimension of a random JordanSpinEJA.
2578 def random_instance(cls
, field
=AA
, **kwargs
):
2580 Return a random instance of this type of algebra.
2582 Needed here to override the implementation for ``BilinearFormEJA``.
2584 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2585 return cls(n
, field
, **kwargs
)
2588 class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra
,
2589 ConcreteEuclideanJordanAlgebra
):
2591 The trivial Euclidean Jordan algebra consisting of only a zero element.
2595 sage: from mjo.eja.eja_algebra import TrivialEJA
2599 sage: J = TrivialEJA()
2606 sage: 7*J.one()*12*J.one()
2608 sage: J.one().inner_product(J.one())
2610 sage: J.one().norm()
2612 sage: J.one().subalgebra_generated_by()
2613 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2618 def __init__(self
, field
=AA
, **kwargs
):
2621 super(TrivialEJA
, self
).__init
__(field
,
2626 # The rank is zero using my definition, namely the dimension of the
2627 # largest subalgebra generated by any element.
2628 self
.rank
.set_cache(0)
2629 self
.one
.set_cache( self
.zero() )
2632 def random_instance(cls
, field
=AA
, **kwargs
):
2633 # We don't take a "size" argument so the superclass method is
2634 # inappropriate for us.
2635 return cls(field
, **kwargs
)
2637 class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra
):
2639 The external (orthogonal) direct sum of two other Euclidean Jordan
2640 algebras. Essentially the Cartesian product of its two factors.
2641 Every Euclidean Jordan algebra decomposes into an orthogonal
2642 direct sum of simple Euclidean Jordan algebras, so no generality
2643 is lost by providing only this construction.
2647 sage: from mjo.eja.eja_algebra import (random_eja,
2649 ....: RealSymmetricEJA,
2654 sage: J1 = HadamardEJA(2)
2655 sage: J2 = RealSymmetricEJA(3)
2656 sage: J = DirectSumEJA(J1,J2)
2664 The external direct sum construction is only valid when the two factors
2665 have the same base ring; an error is raised otherwise::
2667 sage: set_random_seed()
2668 sage: J1 = random_eja(AA)
2669 sage: J2 = random_eja(QQ)
2670 sage: J = DirectSumEJA(J1,J2)
2671 Traceback (most recent call last):
2673 ValueError: algebras must share the same base field
2676 def __init__(self
, J1
, J2
, **kwargs
):
2677 if J1
.base_ring() != J2
.base_ring():
2678 raise ValueError("algebras must share the same base field")
2679 field
= J1
.base_ring()
2681 self
._factors
= (J1
, J2
)
2685 V
= VectorSpace(field
, n
)
2686 mult_table
= [ [ V
.zero() for j
in range(n
) ]
2690 p
= (J1
.monomial(i
)*J1
.monomial(j
)).to_vector()
2691 mult_table
[i
][j
] = V(p
.list() + [field
.zero()]*n2
)
2695 p
= (J2
.monomial(i
)*J2
.monomial(j
)).to_vector()
2696 mult_table
[n1
+i
][n1
+j
] = V([field
.zero()]*n1
+ p
.list())
2698 # TODO: build the IP table here from the two constituent IP
2699 # matrices (it'll be block diagonal, I think).
2701 super(DirectSumEJA
, self
).__init
__(field
,
2706 self
.rank
.set_cache(J1
.rank() + J2
.rank())
2711 Return the pair of this algebra's factors.
2715 sage: from mjo.eja.eja_algebra import (HadamardEJA,
2716 ....: JordanSpinEJA,
2721 sage: J1 = HadamardEJA(2,QQ)
2722 sage: J2 = JordanSpinEJA(3,QQ)
2723 sage: J = DirectSumEJA(J1,J2)
2725 (Euclidean Jordan algebra of dimension 2 over Rational Field,
2726 Euclidean Jordan algebra of dimension 3 over Rational Field)
2729 return self
._factors
2731 def projections(self
):
2733 Return a pair of projections onto this algebra's factors.
2737 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
2738 ....: ComplexHermitianEJA,
2743 sage: J1 = JordanSpinEJA(2)
2744 sage: J2 = ComplexHermitianEJA(2)
2745 sage: J = DirectSumEJA(J1,J2)
2746 sage: (pi_left, pi_right) = J.projections()
2747 sage: J.one().to_vector()
2749 sage: pi_left(J.one()).to_vector()
2751 sage: pi_right(J.one()).to_vector()
2755 (J1
,J2
) = self
.factors()
2758 V_basis
= self
.vector_space().basis()
2759 # Need to specify the dimensions explicitly so that we don't
2760 # wind up with a zero-by-zero matrix when we want e.g. a
2761 # zero-by-two matrix (important for composing things).
2762 P1
= matrix(self
.base_ring(), m
, m
+n
, V_basis
[:m
])
2763 P2
= matrix(self
.base_ring(), n
, m
+n
, V_basis
[m
:])
2764 pi_left
= FiniteDimensionalEuclideanJordanAlgebraOperator(self
,J1
,P1
)
2765 pi_right
= FiniteDimensionalEuclideanJordanAlgebraOperator(self
,J2
,P2
)
2766 return (pi_left
, pi_right
)
2768 def inclusions(self
):
2770 Return the pair of inclusion maps from our factors into us.
2774 sage: from mjo.eja.eja_algebra import (random_eja,
2775 ....: JordanSpinEJA,
2776 ....: RealSymmetricEJA,
2781 sage: J1 = JordanSpinEJA(3)
2782 sage: J2 = RealSymmetricEJA(2)
2783 sage: J = DirectSumEJA(J1,J2)
2784 sage: (iota_left, iota_right) = J.inclusions()
2785 sage: iota_left(J1.zero()) == J.zero()
2787 sage: iota_right(J2.zero()) == J.zero()
2789 sage: J1.one().to_vector()
2791 sage: iota_left(J1.one()).to_vector()
2793 sage: J2.one().to_vector()
2795 sage: iota_right(J2.one()).to_vector()
2797 sage: J.one().to_vector()
2802 Composing a projection with the corresponding inclusion should
2803 produce the identity map, and mismatching them should produce
2806 sage: set_random_seed()
2807 sage: J1 = random_eja()
2808 sage: J2 = random_eja()
2809 sage: J = DirectSumEJA(J1,J2)
2810 sage: (iota_left, iota_right) = J.inclusions()
2811 sage: (pi_left, pi_right) = J.projections()
2812 sage: pi_left*iota_left == J1.one().operator()
2814 sage: pi_right*iota_right == J2.one().operator()
2816 sage: (pi_left*iota_right).is_zero()
2818 sage: (pi_right*iota_left).is_zero()
2822 (J1
,J2
) = self
.factors()
2825 V_basis
= self
.vector_space().basis()
2826 # Need to specify the dimensions explicitly so that we don't
2827 # wind up with a zero-by-zero matrix when we want e.g. a
2828 # two-by-zero matrix (important for composing things).
2829 I1
= matrix
.column(self
.base_ring(), m
, m
+n
, V_basis
[:m
])
2830 I2
= matrix
.column(self
.base_ring(), n
, m
+n
, V_basis
[m
:])
2831 iota_left
= FiniteDimensionalEuclideanJordanAlgebraOperator(J1
,self
,I1
)
2832 iota_right
= FiniteDimensionalEuclideanJordanAlgebraOperator(J2
,self
,I2
)
2833 return (iota_left
, iota_right
)
2835 def inner_product(self
, x
, y
):
2837 The standard Cartesian inner-product.
2839 We project ``x`` and ``y`` onto our factors, and add up the
2840 inner-products from the subalgebras.
2845 sage: from mjo.eja.eja_algebra import (HadamardEJA,
2846 ....: QuaternionHermitianEJA,
2851 sage: J1 = HadamardEJA(3,QQ)
2852 sage: J2 = QuaternionHermitianEJA(2,QQ,normalize_basis=False)
2853 sage: J = DirectSumEJA(J1,J2)
2858 sage: x1.inner_product(x2)
2860 sage: y1.inner_product(y2)
2862 sage: J.one().inner_product(J.one())
2866 (pi_left
, pi_right
) = self
.projections()
2872 return (x1
.inner_product(y1
) + x2
.inner_product(y2
))
2876 random_eja
= ConcreteEuclideanJordanAlgebra
.random_instance