2 Euclidean Jordan Algebras. These are formally-real Jordan Algebras;
3 specifically those where u^2 + v^2 = 0 implies that u = v = 0. They
4 are used in optimization, and have some additional nice methods beyond
5 what can be supported in a general Jordan Algebra.
10 sage: from mjo.eja.eja_algebra import random_eja
15 Euclidean Jordan algebra of dimension...
19 from itertools
import repeat
21 from sage
.algebras
.quatalg
.quaternion_algebra
import QuaternionAlgebra
22 from sage
.categories
.magmatic_algebras
import MagmaticAlgebras
23 from sage
.combinat
.free_module
import CombinatorialFreeModule
24 from sage
.matrix
.constructor
import matrix
25 from sage
.matrix
.matrix_space
import MatrixSpace
26 from sage
.misc
.cachefunc
import cached_method
27 from sage
.misc
.table
import table
28 from sage
.modules
.free_module
import FreeModule
, VectorSpace
29 from sage
.rings
.all
import (ZZ
, QQ
, AA
, QQbar
, RR
, RLF
, CLF
,
32 from mjo
.eja
.eja_element
import FiniteDimensionalEuclideanJordanAlgebraElement
33 from mjo
.eja
.eja_operator
import FiniteDimensionalEuclideanJordanAlgebraOperator
34 from mjo
.eja
.eja_utils
import _mat2vec
36 class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule
):
38 The lowest-level class for representing a Euclidean Jordan algebra.
40 def _coerce_map_from_base_ring(self
):
42 Disable the map from the base ring into the algebra.
44 Performing a nonsense conversion like this automatically
45 is counterpedagogical. The fallback is to try the usual
46 element constructor, which should also fail.
50 sage: from mjo.eja.eja_algebra import random_eja
54 sage: set_random_seed()
55 sage: J = random_eja()
57 Traceback (most recent call last):
59 ValueError: not an element of this algebra
76 * field -- the scalar field for this algebra (must be real)
78 * multiplication_table -- the multiplication table for this
79 algebra's implicit basis. Only the lower-triangular portion
80 of the table is used, since the multiplication is assumed
85 sage: from mjo.eja.eja_algebra import (
86 ....: FiniteDimensionalEuclideanJordanAlgebra,
92 By definition, Jordan multiplication commutes::
94 sage: set_random_seed()
95 sage: J = random_eja()
96 sage: x,y = J.random_elements(2)
100 An error is raised if the Jordan product is not commutative::
102 sage: JP = ((1,2),(0,0))
103 sage: IP = ((1,0),(0,1))
104 sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,JP,IP)
105 Traceback (most recent call last):
107 ValueError: Jordan product is not commutative
109 An error is raised if the inner-product is not commutative::
111 sage: JP = ((1,0),(0,1))
112 sage: IP = ((1,2),(0,0))
113 sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,JP,IP)
114 Traceback (most recent call last):
116 ValueError: inner-product is not commutative
120 The ``field`` we're given must be real with ``check_field=True``::
122 sage: JordanSpinEJA(2, field=QQbar)
123 Traceback (most recent call last):
125 ValueError: scalar field is not real
126 sage: JordanSpinEJA(2, field=QQbar, check_field=False)
127 Euclidean Jordan algebra of dimension 2 over Algebraic Field
129 The multiplication table must be square with ``check_axioms=True``::
131 sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,((),()),((1,),))
132 Traceback (most recent call last):
134 ValueError: multiplication table is not square
136 The multiplication and inner-product tables must be the same
137 size (and in particular, the inner-product table must also be
138 square) with ``check_axioms=True``::
140 sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,((1,),),(()))
141 Traceback (most recent call last):
143 ValueError: multiplication and inner-product tables are
145 sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,((1,),),((1,2),))
146 Traceback (most recent call last):
148 ValueError: multiplication and inner-product tables are
153 if not field
.is_subring(RR
):
154 # Note: this does return true for the real algebraic
155 # field, the rationals, and any quadratic field where
156 # we've specified a real embedding.
157 raise ValueError("scalar field is not real")
160 # The multiplication and inner-product tables should be square
161 # if the user wants us to verify them. And we verify them as
162 # soon as possible, because we want to exploit their symmetry.
163 n
= len(multiplication_table
)
165 if not all( len(l
) == n
for l
in multiplication_table
):
166 raise ValueError("multiplication table is not square")
168 # If the multiplication table is square, we can check if
169 # the inner-product table is square by comparing it to the
170 # multiplication table's dimensions.
171 msg
= "multiplication and inner-product tables are different sizes"
172 if not len(inner_product_table
) == n
:
173 raise ValueError(msg
)
175 if not all( len(l
) == n
for l
in inner_product_table
):
176 raise ValueError(msg
)
178 # Check commutativity of the Jordan product (symmetry of
179 # the multiplication table) and the commutativity of the
180 # inner-product (symmetry of the inner-product table)
181 # first if we're going to check them at all.. This has to
182 # be done before we define product_on_basis(), because
183 # that method assumes that self._multiplication_table is
184 # symmetric. And it has to be done before we build
185 # self._inner_product_matrix, because the process used to
186 # construct it assumes symmetry as well.
187 if not all( multiplication_table
[j
][i
]
188 == multiplication_table
[i
][j
]
190 for j
in range(i
+1) ):
191 raise ValueError("Jordan product is not commutative")
193 if not all( inner_product_table
[j
][i
]
194 == inner_product_table
[i
][j
]
196 for j
in range(i
+1) ):
197 raise ValueError("inner-product is not commutative")
199 self
._matrix
_basis
= matrix_basis
202 category
= MagmaticAlgebras(field
).FiniteDimensional()
203 category
= category
.WithBasis().Unital()
205 fda
= super(FiniteDimensionalEuclideanJordanAlgebra
, self
)
210 self
.print_options(bracket
='')
212 # The multiplication table we're given is necessarily in terms
213 # of vectors, because we don't have an algebra yet for
214 # anything to be an element of. However, it's faster in the
215 # long run to have the multiplication table be in terms of
216 # algebra elements. We do this after calling the superclass
217 # constructor so that from_vector() knows what to do.
219 # Note: we take advantage of symmetry here, and only store
220 # the lower-triangular portion of the table.
221 self
._multiplication
_table
= [ [ self
.vector_space().zero()
222 for j
in range(i
+1) ]
227 elt
= self
.from_vector(multiplication_table
[i
][j
])
228 self
._multiplication
_table
[i
][j
] = elt
230 self
._multiplication
_table
= tuple(map(tuple, self
._multiplication
_table
))
232 # Save our inner product as a matrix, since the efficiency of
233 # matrix multiplication will usually outweigh the fact that we
234 # have to store a redundant upper- or lower-triangular part.
235 # Pre-cache the fact that these are Hermitian (real symmetric,
236 # in fact) in case some e.g. matrix multiplication routine can
237 # take advantage of it.
238 ip_matrix_constructor
= lambda i
,j
: inner_product_table
[i
][j
] if j
<= i
else inner_product_table
[j
][i
]
239 self
._inner
_product
_matrix
= matrix(field
, n
, ip_matrix_constructor
)
240 self
._inner
_product
_matrix
._cache
= {'hermitian': True}
241 self
._inner
_product
_matrix
.set_immutable()
244 if not self
._is
_jordanian
():
245 raise ValueError("Jordan identity does not hold")
246 if not self
._inner
_product
_is
_associative
():
247 raise ValueError("inner product is not associative")
249 def _element_constructor_(self
, elt
):
251 Construct an element of this algebra from its vector or matrix
254 This gets called only after the parent element _call_ method
255 fails to find a coercion for the argument.
259 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
261 ....: RealSymmetricEJA)
265 The identity in `S^n` is converted to the identity in the EJA::
267 sage: J = RealSymmetricEJA(3)
268 sage: I = matrix.identity(QQ,3)
269 sage: J(I) == J.one()
272 This skew-symmetric matrix can't be represented in the EJA::
274 sage: J = RealSymmetricEJA(3)
275 sage: A = matrix(QQ,3, lambda i,j: i-j)
277 Traceback (most recent call last):
279 ValueError: not an element of this algebra
283 Ensure that we can convert any element of the two non-matrix
284 simple algebras (whose matrix representations are columns)
285 back and forth faithfully::
287 sage: set_random_seed()
288 sage: J = HadamardEJA.random_instance()
289 sage: x = J.random_element()
290 sage: J(x.to_vector().column()) == x
292 sage: J = JordanSpinEJA.random_instance()
293 sage: x = J.random_element()
294 sage: J(x.to_vector().column()) == x
298 msg
= "not an element of this algebra"
300 # The superclass implementation of random_element()
301 # needs to be able to coerce "0" into the algebra.
303 elif elt
in self
.base_ring():
304 # Ensure that no base ring -> algebra coercion is performed
305 # by this method. There's some stupidity in sage that would
306 # otherwise propagate to this method; for example, sage thinks
307 # that the integer 3 belongs to the space of 2-by-2 matrices.
308 raise ValueError(msg
)
310 if elt
not in self
.matrix_space():
311 raise ValueError(msg
)
313 # Thanks for nothing! Matrix spaces aren't vector spaces in
314 # Sage, so we have to figure out its matrix-basis coordinates
315 # ourselves. We use the basis space's ring instead of the
316 # element's ring because the basis space might be an algebraic
317 # closure whereas the base ring of the 3-by-3 identity matrix
318 # could be QQ instead of QQbar.
320 # We pass check=False because the matrix basis is "guaranteed"
321 # to be linearly independent... right? Ha ha.
322 V
= VectorSpace(self
.base_ring(), elt
.nrows()*elt
.ncols())
323 W
= V
.span_of_basis( (_mat2vec(s
) for s
in self
.matrix_basis()),
327 coords
= W
.coordinate_vector(_mat2vec(elt
))
328 except ArithmeticError: # vector is not in free module
329 raise ValueError(msg
)
331 return self
.from_vector(coords
)
335 Return a string representation of ``self``.
339 sage: from mjo.eja.eja_algebra import JordanSpinEJA
343 Ensure that it says what we think it says::
345 sage: JordanSpinEJA(2, field=AA)
346 Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
347 sage: JordanSpinEJA(3, field=RDF)
348 Euclidean Jordan algebra of dimension 3 over Real Double Field
351 fmt
= "Euclidean Jordan algebra of dimension {} over {}"
352 return fmt
.format(self
.dimension(), self
.base_ring())
354 def product_on_basis(self
, i
, j
):
355 # We only stored the lower-triangular portion of the
356 # multiplication table.
358 return self
._multiplication
_table
[i
][j
]
360 return self
._multiplication
_table
[j
][i
]
362 def _is_commutative(self
):
364 Whether or not this algebra's multiplication table is commutative.
366 This method should of course always return ``True``, unless
367 this algebra was constructed with ``check_axioms=False`` and
368 passed an invalid multiplication table.
370 return all( self
.product_on_basis(i
,j
) == self
.product_on_basis(i
,j
)
371 for i
in range(self
.dimension())
372 for j
in range(self
.dimension()) )
374 def _is_jordanian(self
):
376 Whether or not this algebra's multiplication table respects the
377 Jordan identity `(x^{2})(xy) = x(x^{2}y)`.
379 We only check one arrangement of `x` and `y`, so for a
380 ``True`` result to be truly true, you should also check
381 :meth:`_is_commutative`. This method should of course always
382 return ``True``, unless this algebra was constructed with
383 ``check_axioms=False`` and passed an invalid multiplication table.
385 return all( (self
.monomial(i
)**2)*(self
.monomial(i
)*self
.monomial(j
))
387 (self
.monomial(i
))*((self
.monomial(i
)**2)*self
.monomial(j
))
388 for i
in range(self
.dimension())
389 for j
in range(self
.dimension()) )
391 def _inner_product_is_associative(self
):
393 Return whether or not this algebra's inner product `B` is
394 associative; that is, whether or not `B(xy,z) = B(x,yz)`.
396 This method should of course always return ``True``, unless
397 this algebra was constructed with ``check_axioms=False`` and
398 passed an invalid multiplication table.
401 # Used to check whether or not something is zero in an inexact
402 # ring. This number is sufficient to allow the construction of
403 # QuaternionHermitianEJA(2, field=RDF) with check_axioms=True.
406 for i
in range(self
.dimension()):
407 for j
in range(self
.dimension()):
408 for k
in range(self
.dimension()):
412 diff
= (x
*y
).inner_product(z
) - x
.inner_product(y
*z
)
414 if self
.base_ring().is_exact():
418 if diff
.abs() > epsilon
:
424 def characteristic_polynomial_of(self
):
426 Return the algebra's "characteristic polynomial of" function,
427 which is itself a multivariate polynomial that, when evaluated
428 at the coordinates of some algebra element, returns that
429 element's characteristic polynomial.
431 The resulting polynomial has `n+1` variables, where `n` is the
432 dimension of this algebra. The first `n` variables correspond to
433 the coordinates of an algebra element: when evaluated at the
434 coordinates of an algebra element with respect to a certain
435 basis, the result is a univariate polynomial (in the one
436 remaining variable ``t``), namely the characteristic polynomial
441 sage: from mjo.eja.eja_algebra import JordanSpinEJA, TrivialEJA
445 The characteristic polynomial in the spin algebra is given in
446 Alizadeh, Example 11.11::
448 sage: J = JordanSpinEJA(3)
449 sage: p = J.characteristic_polynomial_of(); p
450 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
451 sage: xvec = J.one().to_vector()
455 By definition, the characteristic polynomial is a monic
456 degree-zero polynomial in a rank-zero algebra. Note that
457 Cayley-Hamilton is indeed satisfied since the polynomial
458 ``1`` evaluates to the identity element of the algebra on
461 sage: J = TrivialEJA()
462 sage: J.characteristic_polynomial_of()
469 # The list of coefficient polynomials a_0, a_1, a_2, ..., a_(r-1).
470 a
= self
._charpoly
_coefficients
()
472 # We go to a bit of trouble here to reorder the
473 # indeterminates, so that it's easier to evaluate the
474 # characteristic polynomial at x's coordinates and get back
475 # something in terms of t, which is what we want.
476 S
= PolynomialRing(self
.base_ring(),'t')
480 S
= PolynomialRing(S
, R
.variable_names())
483 return (t
**r
+ sum( a
[k
]*(t
**k
) for k
in range(r
) ))
485 def coordinate_polynomial_ring(self
):
487 The multivariate polynomial ring in which this algebra's
488 :meth:`characteristic_polynomial_of` lives.
492 sage: from mjo.eja.eja_algebra import (HadamardEJA,
493 ....: RealSymmetricEJA)
497 sage: J = HadamardEJA(2)
498 sage: J.coordinate_polynomial_ring()
499 Multivariate Polynomial Ring in X1, X2...
500 sage: J = RealSymmetricEJA(3,field=QQ,orthonormalize=False)
501 sage: J.coordinate_polynomial_ring()
502 Multivariate Polynomial Ring in X1, X2, X3, X4, X5, X6...
505 var_names
= tuple( "X%d" % z
for z
in range(1, self
.dimension()+1) )
506 return PolynomialRing(self
.base_ring(), var_names
)
508 def inner_product(self
, x
, y
):
510 The inner product associated with this Euclidean Jordan algebra.
512 Defaults to the trace inner product, but can be overridden by
513 subclasses if they are sure that the necessary properties are
518 sage: from mjo.eja.eja_algebra import (random_eja,
520 ....: BilinearFormEJA)
524 Our inner product is "associative," which means the following for
525 a symmetric bilinear form::
527 sage: set_random_seed()
528 sage: J = random_eja()
529 sage: x,y,z = J.random_elements(3)
530 sage: (x*y).inner_product(z) == y.inner_product(x*z)
535 Ensure that this is the usual inner product for the algebras
538 sage: set_random_seed()
539 sage: J = HadamardEJA.random_instance()
540 sage: x,y = J.random_elements(2)
541 sage: actual = x.inner_product(y)
542 sage: expected = x.to_vector().inner_product(y.to_vector())
543 sage: actual == expected
546 Ensure that this is one-half of the trace inner-product in a
547 BilinearFormEJA that isn't just the reals (when ``n`` isn't
548 one). This is in Faraut and Koranyi, and also my "On the
551 sage: set_random_seed()
552 sage: J = BilinearFormEJA.random_instance()
553 sage: n = J.dimension()
554 sage: x = J.random_element()
555 sage: y = J.random_element()
556 sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
559 B
= self
._inner
_product
_matrix
560 return (B
*x
.to_vector()).inner_product(y
.to_vector())
563 def is_trivial(self
):
565 Return whether or not this algebra is trivial.
567 A trivial algebra contains only the zero element.
571 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
576 sage: J = ComplexHermitianEJA(3)
582 sage: J = TrivialEJA()
587 return self
.dimension() == 0
590 def multiplication_table(self
):
592 Return a visual representation of this algebra's multiplication
593 table (on basis elements).
597 sage: from mjo.eja.eja_algebra import JordanSpinEJA
601 sage: J = JordanSpinEJA(4)
602 sage: J.multiplication_table()
603 +----++----+----+----+----+
604 | * || e0 | e1 | e2 | e3 |
605 +====++====+====+====+====+
606 | e0 || e0 | e1 | e2 | e3 |
607 +----++----+----+----+----+
608 | e1 || e1 | e0 | 0 | 0 |
609 +----++----+----+----+----+
610 | e2 || e2 | 0 | e0 | 0 |
611 +----++----+----+----+----+
612 | e3 || e3 | 0 | 0 | e0 |
613 +----++----+----+----+----+
617 # Prepend the header row.
618 M
= [["*"] + list(self
.gens())]
620 # And to each subsequent row, prepend an entry that belongs to
621 # the left-side "header column."
622 M
+= [ [self
.monomial(i
)] + [ self
.product_on_basis(i
,j
)
626 return table(M
, header_row
=True, header_column
=True, frame
=True)
629 def matrix_basis(self
):
631 Return an (often more natural) representation of this algebras
632 basis as an ordered tuple of matrices.
634 Every finite-dimensional Euclidean Jordan Algebra is a, up to
635 Jordan isomorphism, a direct sum of five simple
636 algebras---four of which comprise Hermitian matrices. And the
637 last type of algebra can of course be thought of as `n`-by-`1`
638 column matrices (ambiguusly called column vectors) to avoid
639 special cases. As a result, matrices (and column vectors) are
640 a natural representation format for Euclidean Jordan algebra
643 But, when we construct an algebra from a basis of matrices,
644 those matrix representations are lost in favor of coordinate
645 vectors *with respect to* that basis. We could eventually
646 convert back if we tried hard enough, but having the original
647 representations handy is valuable enough that we simply store
648 them and return them from this method.
650 Why implement this for non-matrix algebras? Avoiding special
651 cases for the :class:`BilinearFormEJA` pays with simplicity in
652 its own right. But mainly, we would like to be able to assume
653 that elements of a :class:`DirectSumEJA` can be displayed
654 nicely, without having to have special classes for direct sums
655 one of whose components was a matrix algebra.
659 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
660 ....: RealSymmetricEJA)
664 sage: J = RealSymmetricEJA(2)
666 Finite family {0: e0, 1: e1, 2: e2}
667 sage: J.matrix_basis()
669 [1 0] [ 0 0.7071067811865475?] [0 0]
670 [0 0], [0.7071067811865475? 0], [0 1]
675 sage: J = JordanSpinEJA(2)
677 Finite family {0: e0, 1: e1}
678 sage: J.matrix_basis()
684 if self
._matrix
_basis
is None:
685 M
= self
.matrix_space()
686 return tuple( M(b
.to_vector()) for b
in self
.basis() )
688 return self
._matrix
_basis
691 def matrix_space(self
):
693 Return the matrix space in which this algebra's elements live, if
694 we think of them as matrices (including column vectors of the
697 Generally this will be an `n`-by-`1` column-vector space,
698 except when the algebra is trivial. There it's `n`-by-`n`
699 (where `n` is zero), to ensure that two elements of the matrix
700 space (empty matrices) can be multiplied.
702 Matrix algebras override this with something more useful.
704 if self
.is_trivial():
705 return MatrixSpace(self
.base_ring(), 0)
706 elif self
._matrix
_basis
is None or len(self
._matrix
_basis
) == 0:
707 return MatrixSpace(self
.base_ring(), self
.dimension(), 1)
709 return self
._matrix
_basis
[0].matrix_space()
715 Return the unit element of this algebra.
719 sage: from mjo.eja.eja_algebra import (HadamardEJA,
724 sage: J = HadamardEJA(5)
726 e0 + e1 + e2 + e3 + e4
730 The identity element acts like the identity::
732 sage: set_random_seed()
733 sage: J = random_eja()
734 sage: x = J.random_element()
735 sage: J.one()*x == x and x*J.one() == x
738 The matrix of the unit element's operator is the identity::
740 sage: set_random_seed()
741 sage: J = random_eja()
742 sage: actual = J.one().operator().matrix()
743 sage: expected = matrix.identity(J.base_ring(), J.dimension())
744 sage: actual == expected
747 Ensure that the cached unit element (often precomputed by
748 hand) agrees with the computed one::
750 sage: set_random_seed()
751 sage: J = random_eja()
752 sage: cached = J.one()
753 sage: J.one.clear_cache()
754 sage: J.one() == cached
758 # We can brute-force compute the matrices of the operators
759 # that correspond to the basis elements of this algebra.
760 # If some linear combination of those basis elements is the
761 # algebra identity, then the same linear combination of
762 # their matrices has to be the identity matrix.
764 # Of course, matrices aren't vectors in sage, so we have to
765 # appeal to the "long vectors" isometry.
766 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
768 # Now we use basic linear algebra to find the coefficients,
769 # of the matrices-as-vectors-linear-combination, which should
770 # work for the original algebra basis too.
771 A
= matrix(self
.base_ring(), oper_vecs
)
773 # We used the isometry on the left-hand side already, but we
774 # still need to do it for the right-hand side. Recall that we
775 # wanted something that summed to the identity matrix.
776 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
778 # Now if there's an identity element in the algebra, this
779 # should work. We solve on the left to avoid having to
780 # transpose the matrix "A".
781 return self
.from_vector(A
.solve_left(b
))
784 def peirce_decomposition(self
, c
):
786 The Peirce decomposition of this algebra relative to the
789 In the future, this can be extended to a complete system of
790 orthogonal idempotents.
794 - ``c`` -- an idempotent of this algebra.
798 A triple (J0, J5, J1) containing two subalgebras and one subspace
801 - ``J0`` -- the algebra on the eigenspace of ``c.operator()``
802 corresponding to the eigenvalue zero.
804 - ``J5`` -- the eigenspace (NOT a subalgebra) of ``c.operator()``
805 corresponding to the eigenvalue one-half.
807 - ``J1`` -- the algebra on the eigenspace of ``c.operator()``
808 corresponding to the eigenvalue one.
810 These are the only possible eigenspaces for that operator, and this
811 algebra is a direct sum of them. The spaces ``J0`` and ``J1`` are
812 orthogonal, and are subalgebras of this algebra with the appropriate
817 sage: from mjo.eja.eja_algebra import random_eja, RealSymmetricEJA
821 The canonical example comes from the symmetric matrices, which
822 decompose into diagonal and off-diagonal parts::
824 sage: J = RealSymmetricEJA(3)
825 sage: C = matrix(QQ, [ [1,0,0],
829 sage: J0,J5,J1 = J.peirce_decomposition(c)
831 Euclidean Jordan algebra of dimension 1...
833 Vector space of degree 6 and dimension 2...
835 Euclidean Jordan algebra of dimension 3...
836 sage: J0.one().to_matrix()
840 sage: orig_df = AA.options.display_format
841 sage: AA.options.display_format = 'radical'
842 sage: J.from_vector(J5.basis()[0]).to_matrix()
846 sage: J.from_vector(J5.basis()[1]).to_matrix()
850 sage: AA.options.display_format = orig_df
851 sage: J1.one().to_matrix()
858 Every algebra decomposes trivially with respect to its identity
861 sage: set_random_seed()
862 sage: J = random_eja()
863 sage: J0,J5,J1 = J.peirce_decomposition(J.one())
864 sage: J0.dimension() == 0 and J5.dimension() == 0
866 sage: J1.superalgebra() == J and J1.dimension() == J.dimension()
869 The decomposition is into eigenspaces, and its components are
870 therefore necessarily orthogonal. Moreover, the identity
871 elements in the two subalgebras are the projections onto their
872 respective subspaces of the superalgebra's identity element::
874 sage: set_random_seed()
875 sage: J = random_eja()
876 sage: x = J.random_element()
877 sage: if not J.is_trivial():
878 ....: while x.is_nilpotent():
879 ....: x = J.random_element()
880 sage: c = x.subalgebra_idempotent()
881 sage: J0,J5,J1 = J.peirce_decomposition(c)
883 sage: for (w,y,z) in zip(J0.basis(), J5.basis(), J1.basis()):
884 ....: w = w.superalgebra_element()
885 ....: y = J.from_vector(y)
886 ....: z = z.superalgebra_element()
887 ....: ipsum += w.inner_product(y).abs()
888 ....: ipsum += w.inner_product(z).abs()
889 ....: ipsum += y.inner_product(z).abs()
892 sage: J1(c) == J1.one()
894 sage: J0(J.one() - c) == J0.one()
898 if not c
.is_idempotent():
899 raise ValueError("element is not idempotent: %s" % c
)
901 from mjo
.eja
.eja_subalgebra
import FiniteDimensionalEuclideanJordanSubalgebra
903 # Default these to what they should be if they turn out to be
904 # trivial, because eigenspaces_left() won't return eigenvalues
905 # corresponding to trivial spaces (e.g. it returns only the
906 # eigenspace corresponding to lambda=1 if you take the
907 # decomposition relative to the identity element).
908 trivial
= FiniteDimensionalEuclideanJordanSubalgebra(self
, ())
909 J0
= trivial
# eigenvalue zero
910 J5
= VectorSpace(self
.base_ring(), 0) # eigenvalue one-half
911 J1
= trivial
# eigenvalue one
913 for (eigval
, eigspace
) in c
.operator().matrix().right_eigenspaces():
914 if eigval
== ~
(self
.base_ring()(2)):
917 gens
= tuple( self
.from_vector(b
) for b
in eigspace
.basis() )
918 subalg
= FiniteDimensionalEuclideanJordanSubalgebra(self
,
926 raise ValueError("unexpected eigenvalue: %s" % eigval
)
931 def random_element(self
, thorough
=False):
933 Return a random element of this algebra.
935 Our algebra superclass method only returns a linear
936 combination of at most two basis elements. We instead
937 want the vector space "random element" method that
938 returns a more diverse selection.
942 - ``thorough`` -- (boolean; default False) whether or not we
943 should generate irrational coefficients for the random
944 element when our base ring is irrational; this slows the
945 algebra operations to a crawl, but any truly random method
949 # For a general base ring... maybe we can trust this to do the
950 # right thing? Unlikely, but.
951 V
= self
.vector_space()
952 v
= V
.random_element()
954 if self
.base_ring() is AA
:
955 # The "random element" method of the algebraic reals is
956 # stupid at the moment, and only returns integers between
957 # -2 and 2, inclusive:
959 # https://trac.sagemath.org/ticket/30875
961 # Instead, we implement our own "random vector" method,
962 # and then coerce that into the algebra. We use the vector
963 # space degree here instead of the dimension because a
964 # subalgebra could (for example) be spanned by only two
965 # vectors, each with five coordinates. We need to
966 # generate all five coordinates.
968 v
*= QQbar
.random_element().real()
970 v
*= QQ
.random_element()
972 return self
.from_vector(V
.coordinate_vector(v
))
974 def random_elements(self
, count
, thorough
=False):
976 Return ``count`` random elements as a tuple.
980 - ``thorough`` -- (boolean; default False) whether or not we
981 should generate irrational coefficients for the random
982 elements when our base ring is irrational; this slows the
983 algebra operations to a crawl, but any truly random method
988 sage: from mjo.eja.eja_algebra import JordanSpinEJA
992 sage: J = JordanSpinEJA(3)
993 sage: x,y,z = J.random_elements(3)
994 sage: all( [ x in J, y in J, z in J ])
996 sage: len( J.random_elements(10) ) == 10
1000 return tuple( self
.random_element(thorough
)
1001 for idx
in range(count
) )
1005 def _charpoly_coefficients(self
):
1007 The `r` polynomial coefficients of the "characteristic polynomial
1010 n
= self
.dimension()
1011 R
= self
.coordinate_polynomial_ring()
1013 F
= R
.fraction_field()
1016 # From a result in my book, these are the entries of the
1017 # basis representation of L_x.
1018 return sum( vars[k
]*self
.monomial(k
).operator().matrix()[i
,j
]
1021 L_x
= matrix(F
, n
, n
, L_x_i_j
)
1024 if self
.rank
.is_in_cache():
1026 # There's no need to pad the system with redundant
1027 # columns if we *know* they'll be redundant.
1030 # Compute an extra power in case the rank is equal to
1031 # the dimension (otherwise, we would stop at x^(r-1)).
1032 x_powers
= [ (L_x
**k
)*self
.one().to_vector()
1033 for k
in range(n
+1) ]
1034 A
= matrix
.column(F
, x_powers
[:n
])
1035 AE
= A
.extended_echelon_form()
1042 # The theory says that only the first "r" coefficients are
1043 # nonzero, and they actually live in the original polynomial
1044 # ring and not the fraction field. We negate them because
1045 # in the actual characteristic polynomial, they get moved
1046 # to the other side where x^r lives.
1047 return -A_rref
.solve_right(E
*b
).change_ring(R
)[:r
]
1052 Return the rank of this EJA.
1054 This is a cached method because we know the rank a priori for
1055 all of the algebras we can construct. Thus we can avoid the
1056 expensive ``_charpoly_coefficients()`` call unless we truly
1057 need to compute the whole characteristic polynomial.
1061 sage: from mjo.eja.eja_algebra import (HadamardEJA,
1062 ....: JordanSpinEJA,
1063 ....: RealSymmetricEJA,
1064 ....: ComplexHermitianEJA,
1065 ....: QuaternionHermitianEJA,
1070 The rank of the Jordan spin algebra is always two::
1072 sage: JordanSpinEJA(2).rank()
1074 sage: JordanSpinEJA(3).rank()
1076 sage: JordanSpinEJA(4).rank()
1079 The rank of the `n`-by-`n` Hermitian real, complex, or
1080 quaternion matrices is `n`::
1082 sage: RealSymmetricEJA(4).rank()
1084 sage: ComplexHermitianEJA(3).rank()
1086 sage: QuaternionHermitianEJA(2).rank()
1091 Ensure that every EJA that we know how to construct has a
1092 positive integer rank, unless the algebra is trivial in
1093 which case its rank will be zero::
1095 sage: set_random_seed()
1096 sage: J = random_eja()
1100 sage: r > 0 or (r == 0 and J.is_trivial())
1103 Ensure that computing the rank actually works, since the ranks
1104 of all simple algebras are known and will be cached by default::
1106 sage: set_random_seed() # long time
1107 sage: J = random_eja() # long time
1108 sage: caches = J.rank() # long time
1109 sage: J.rank.clear_cache() # long time
1110 sage: J.rank() == cached # long time
1114 return len(self
._charpoly
_coefficients
())
1117 def vector_space(self
):
1119 Return the vector space that underlies this algebra.
1123 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1127 sage: J = RealSymmetricEJA(2)
1128 sage: J.vector_space()
1129 Vector space of dimension 3 over...
1132 return self
.zero().to_vector().parent().ambient_vector_space()
1135 Element
= FiniteDimensionalEuclideanJordanAlgebraElement
1137 class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra
):
1139 New class for algebras whose supplied basis elements have all rational entries.
1143 sage: from mjo.eja.eja_algebra import BilinearFormEJA
1147 The supplied basis is orthonormalized by default::
1149 sage: B = matrix(QQ, [[1, 0, 0], [0, 25, -32], [0, -32, 41]])
1150 sage: J = BilinearFormEJA(B)
1151 sage: J.matrix_basis()
1164 orthonormalize
=True,
1171 # Abuse the check_field parameter to check that the entries of
1172 # out basis (in ambient coordinates) are in the field QQ.
1173 if not all( all(b_i
in QQ
for b_i
in b
.list()) for b
in basis
):
1174 raise TypeError("basis not rational")
1176 # Temporary(?) hack to ensure that the matrix and vector bases
1177 # are over the same ring.
1178 basis
= tuple( b
.change_ring(field
) for b
in basis
)
1181 vector_basis
= basis
1183 from sage
.structure
.element
import is_Matrix
1184 basis_is_matrices
= False
1188 if is_Matrix(basis
[0]):
1189 basis_is_matrices
= True
1190 from mjo
.eja
.eja_utils
import _vec2mat
1191 vector_basis
= tuple( map(_mat2vec
,basis
) )
1192 degree
= basis
[0].nrows()**2
1194 degree
= basis
[0].degree()
1196 V
= VectorSpace(field
, degree
)
1198 # Save a copy of an algebra with the original, rational basis
1199 # and over QQ where computations are fast.
1200 self
._rational
_algebra
= None
1203 # There's no point in constructing the extra algebra if this
1204 # one is already rational.
1206 # Note: the same Jordan and inner-products work here,
1207 # because they are necessarily defined with respect to
1208 # ambient coordinates and not any particular basis.
1209 self
._rational
_algebra
= RationalBasisEuclideanJordanAlgebra(
1214 orthonormalize
=False,
1221 # Compute the deorthonormalized tables before we orthonormalize
1222 # the given basis. The "check" parameter here guarantees that
1223 # the basis is linearly-independent.
1224 W
= V
.span_of_basis( vector_basis
, check
=check_axioms
)
1226 # Note: the Jordan and inner-products are defined in terms
1227 # of the ambient basis. It's important that their arguments
1228 # are in ambient coordinates as well.
1230 for j
in range(i
+1):
1231 # given basis w.r.t. ambient coords
1232 q_i
= vector_basis
[i
]
1233 q_j
= vector_basis
[j
]
1235 if basis_is_matrices
:
1239 elt
= jordan_product(q_i
, q_j
)
1240 ip
= inner_product(q_i
, q_j
)
1242 if basis_is_matrices
:
1243 # do another mat2vec because the multiplication
1244 # table is in terms of vectors
1247 # We overwrite the name "vector_basis" in a second, but never modify it
1248 # in place, to this effectively makes a copy of it.
1249 deortho_vector_basis
= vector_basis
1250 self
._deortho
_matrix
= None
1253 from mjo
.eja
.eja_utils
import gram_schmidt
1254 if basis_is_matrices
:
1255 vector_ip
= lambda x
,y
: inner_product(_vec2mat(x
), _vec2mat(y
))
1256 vector_basis
= gram_schmidt(vector_basis
, vector_ip
)
1258 vector_basis
= gram_schmidt(vector_basis
, inner_product
)
1260 # Normalize the "matrix" basis, too!
1261 basis
= vector_basis
1263 if basis_is_matrices
:
1264 basis
= tuple( map(_vec2mat
,basis
) )
1266 W
= V
.span_of_basis( vector_basis
, check
=check_axioms
)
1268 # Now "W" is the vector space of our algebra coordinates. The
1269 # variables "X1", "X2",... refer to the entries of vectors in
1270 # W. Thus to convert back and forth between the orthonormal
1271 # coordinates and the given ones, we need to stick the original
1273 U
= V
.span_of_basis( deortho_vector_basis
, check
=check_axioms
)
1274 self
._deortho
_matrix
= matrix( U
.coordinate_vector(q
)
1275 for q
in vector_basis
)
1277 # If the superclass constructor is going to verify the
1278 # symmetry of this table, it has better at least be
1281 mult_table
= [ [0 for j
in range(n
)] for i
in range(n
) ]
1282 ip_table
= [ [0 for j
in range(n
)] for i
in range(n
) ]
1284 mult_table
= [ [0 for j
in range(i
+1)] for i
in range(n
) ]
1285 ip_table
= [ [0 for j
in range(i
+1)] for i
in range(n
) ]
1287 # Note: the Jordan and inner-products are defined in terms
1288 # of the ambient basis. It's important that their arguments
1289 # are in ambient coordinates as well.
1291 for j
in range(i
+1):
1292 # ortho basis w.r.t. ambient coords
1293 q_i
= vector_basis
[i
]
1294 q_j
= vector_basis
[j
]
1296 if basis_is_matrices
:
1300 elt
= jordan_product(q_i
, q_j
)
1301 ip
= inner_product(q_i
, q_j
)
1303 if basis_is_matrices
:
1304 # do another mat2vec because the multiplication
1305 # table is in terms of vectors
1308 elt
= W
.coordinate_vector(elt
)
1309 mult_table
[i
][j
] = elt
1312 # The tables are square if we're verifying that they
1314 mult_table
[j
][i
] = elt
1317 if basis_is_matrices
:
1321 basis
= tuple( x
.column() for x
in basis
)
1323 super().__init
__(field
,
1328 basis
, # matrix basis
1333 def _charpoly_coefficients(self
):
1337 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
1338 ....: JordanSpinEJA)
1342 The base ring of the resulting polynomial coefficients is what
1343 it should be, and not the rationals (unless the algebra was
1344 already over the rationals)::
1346 sage: J = JordanSpinEJA(3)
1347 sage: J._charpoly_coefficients()
1348 (X1^2 - X2^2 - X3^2, -2*X1)
1349 sage: a0 = J._charpoly_coefficients()[0]
1351 Algebraic Real Field
1352 sage: a0.base_ring()
1353 Algebraic Real Field
1356 if self
._rational
_algebra
is None:
1357 # There's no need to construct *another* algebra over the
1358 # rationals if this one is already over the
1359 # rationals. Likewise, if we never orthonormalized our
1360 # basis, we might as well just use the given one.
1361 superclass
= super(RationalBasisEuclideanJordanAlgebra
, self
)
1362 return superclass
._charpoly
_coefficients
()
1364 # Do the computation over the rationals. The answer will be
1365 # the same, because all we've done is a change of basis.
1366 # Then, change back from QQ to our real base ring
1367 a
= ( a_i
.change_ring(self
.base_ring())
1368 for a_i
in self
._rational
_algebra
._charpoly
_coefficients
() )
1370 # Now convert the coordinate variables back to the
1371 # deorthonormalized ones.
1372 R
= self
.coordinate_polynomial_ring()
1373 from sage
.modules
.free_module_element
import vector
1374 X
= vector(R
, R
.gens())
1375 BX
= self
._deortho
_matrix
*X
1377 subs_dict
= { X[i]: BX[i] for i in range(len(X)) }
1378 return tuple( a_i
.subs(subs_dict
) for a_i
in a
)
1380 class ConcreteEuclideanJordanAlgebra(RationalBasisEuclideanJordanAlgebra
):
1382 A class for the Euclidean Jordan algebras that we know by name.
1384 These are the Jordan algebras whose basis, multiplication table,
1385 rank, and so on are known a priori. More to the point, they are
1386 the Euclidean Jordan algebras for which we are able to conjure up
1387 a "random instance."
1391 sage: from mjo.eja.eja_algebra import ConcreteEuclideanJordanAlgebra
1395 Our basis is normalized with respect to the algebra's inner
1396 product, unless we specify otherwise::
1398 sage: set_random_seed()
1399 sage: J = ConcreteEuclideanJordanAlgebra.random_instance()
1400 sage: all( b.norm() == 1 for b in J.gens() )
1403 Since our basis is orthonormal with respect to the algebra's inner
1404 product, and since we know that this algebra is an EJA, any
1405 left-multiplication operator's matrix will be symmetric because
1406 natural->EJA basis representation is an isometry and within the
1407 EJA the operator is self-adjoint by the Jordan axiom::
1409 sage: set_random_seed()
1410 sage: J = ConcreteEuclideanJordanAlgebra.random_instance()
1411 sage: x = J.random_element()
1412 sage: x.operator().is_self_adjoint()
1417 def _max_random_instance_size():
1419 Return an integer "size" that is an upper bound on the size of
1420 this algebra when it is used in a random test
1421 case. Unfortunately, the term "size" is ambiguous -- when
1422 dealing with `R^n` under either the Hadamard or Jordan spin
1423 product, the "size" refers to the dimension `n`. When dealing
1424 with a matrix algebra (real symmetric or complex/quaternion
1425 Hermitian), it refers to the size of the matrix, which is far
1426 less than the dimension of the underlying vector space.
1428 This method must be implemented in each subclass.
1430 raise NotImplementedError
1433 def random_instance(cls
, *args
, **kwargs
):
1435 Return a random instance of this type of algebra.
1437 This method should be implemented in each subclass.
1439 from sage
.misc
.prandom
import choice
1440 eja_class
= choice(cls
.__subclasses
__())
1442 # These all bubble up to the RationalBasisEuclideanJordanAlgebra
1443 # superclass constructor, so any (kw)args valid there are also
1445 return eja_class
.random_instance(*args
, **kwargs
)
1448 class MatrixEuclideanJordanAlgebra
:
1450 def dimension_over_reals():
1452 The dimension of this matrix's base ring over the reals.
1454 The reals are dimension one over themselves, obviously; that's
1455 just `\mathbb{R}^{1}`. Likewise, the complex numbers `a + bi`
1456 have dimension two. Finally, the quaternions have dimension
1457 four over the reals.
1459 This is used to determine the size of the matrix returned from
1460 :meth:`real_embed`, among other things.
1462 raise NotImplementedError
1465 def real_embed(cls
,M
):
1467 Embed the matrix ``M`` into a space of real matrices.
1469 The matrix ``M`` can have entries in any field at the moment:
1470 the real numbers, complex numbers, or quaternions. And although
1471 they are not a field, we can probably support octonions at some
1472 point, too. This function returns a real matrix that "acts like"
1473 the original with respect to matrix multiplication; i.e.
1475 real_embed(M*N) = real_embed(M)*real_embed(N)
1478 if M
.ncols() != M
.nrows():
1479 raise ValueError("the matrix 'M' must be square")
1484 def real_unembed(cls
,M
):
1486 The inverse of :meth:`real_embed`.
1488 if M
.ncols() != M
.nrows():
1489 raise ValueError("the matrix 'M' must be square")
1490 if not ZZ(M
.nrows()).mod(cls
.dimension_over_reals()).is_zero():
1491 raise ValueError("the matrix 'M' must be a real embedding")
1495 def jordan_product(X
,Y
):
1496 return (X
*Y
+ Y
*X
)/2
1499 def trace_inner_product(cls
,X
,Y
):
1501 Compute the trace inner-product of two real-embeddings.
1505 sage: from mjo.eja.eja_algebra import (RealSymmetricEJA,
1506 ....: ComplexHermitianEJA,
1507 ....: QuaternionHermitianEJA)
1511 This gives the same answer as it would if we computed the trace
1512 from the unembedded (original) matrices::
1514 sage: set_random_seed()
1515 sage: J = RealSymmetricEJA.random_instance()
1516 sage: x,y = J.random_elements(2)
1517 sage: Xe = x.to_matrix()
1518 sage: Ye = y.to_matrix()
1519 sage: X = J.real_unembed(Xe)
1520 sage: Y = J.real_unembed(Ye)
1521 sage: expected = (X*Y).trace()
1522 sage: actual = J.trace_inner_product(Xe,Ye)
1523 sage: actual == expected
1528 sage: set_random_seed()
1529 sage: J = ComplexHermitianEJA.random_instance()
1530 sage: x,y = J.random_elements(2)
1531 sage: Xe = x.to_matrix()
1532 sage: Ye = y.to_matrix()
1533 sage: X = J.real_unembed(Xe)
1534 sage: Y = J.real_unembed(Ye)
1535 sage: expected = (X*Y).trace().real()
1536 sage: actual = J.trace_inner_product(Xe,Ye)
1537 sage: actual == expected
1542 sage: set_random_seed()
1543 sage: J = QuaternionHermitianEJA.random_instance()
1544 sage: x,y = J.random_elements(2)
1545 sage: Xe = x.to_matrix()
1546 sage: Ye = y.to_matrix()
1547 sage: X = J.real_unembed(Xe)
1548 sage: Y = J.real_unembed(Ye)
1549 sage: expected = (X*Y).trace().coefficient_tuple()[0]
1550 sage: actual = J.trace_inner_product(Xe,Ye)
1551 sage: actual == expected
1555 Xu
= cls
.real_unembed(X
)
1556 Yu
= cls
.real_unembed(Y
)
1557 tr
= (Xu
*Yu
).trace()
1560 # Works in QQ, AA, RDF, et cetera.
1562 except AttributeError:
1563 # A quaternion doesn't have a real() method, but does
1564 # have coefficient_tuple() method that returns the
1565 # coefficients of 1, i, j, and k -- in that order.
1566 return tr
.coefficient_tuple()[0]
1569 class RealMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1571 def dimension_over_reals():
1575 class RealSymmetricEJA(ConcreteEuclideanJordanAlgebra
,
1576 RealMatrixEuclideanJordanAlgebra
):
1578 The rank-n simple EJA consisting of real symmetric n-by-n
1579 matrices, the usual symmetric Jordan product, and the trace inner
1580 product. It has dimension `(n^2 + n)/2` over the reals.
1584 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1588 sage: J = RealSymmetricEJA(2)
1589 sage: e0, e1, e2 = J.gens()
1597 In theory, our "field" can be any subfield of the reals::
1599 sage: RealSymmetricEJA(2, field=RDF)
1600 Euclidean Jordan algebra of dimension 3 over Real Double Field
1601 sage: RealSymmetricEJA(2, field=RR)
1602 Euclidean Jordan algebra of dimension 3 over Real Field with
1603 53 bits of precision
1607 The dimension of this algebra is `(n^2 + n) / 2`::
1609 sage: set_random_seed()
1610 sage: n_max = RealSymmetricEJA._max_random_instance_size()
1611 sage: n = ZZ.random_element(1, n_max)
1612 sage: J = RealSymmetricEJA(n)
1613 sage: J.dimension() == (n^2 + n)/2
1616 The Jordan multiplication is what we think it is::
1618 sage: set_random_seed()
1619 sage: J = RealSymmetricEJA.random_instance()
1620 sage: x,y = J.random_elements(2)
1621 sage: actual = (x*y).to_matrix()
1622 sage: X = x.to_matrix()
1623 sage: Y = y.to_matrix()
1624 sage: expected = (X*Y + Y*X)/2
1625 sage: actual == expected
1627 sage: J(expected) == x*y
1630 We can change the generator prefix::
1632 sage: RealSymmetricEJA(3, prefix='q').gens()
1633 (q0, q1, q2, q3, q4, q5)
1635 We can construct the (trivial) algebra of rank zero::
1637 sage: RealSymmetricEJA(0)
1638 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1642 def _denormalized_basis(cls
, n
):
1644 Return a basis for the space of real symmetric n-by-n matrices.
1648 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1652 sage: set_random_seed()
1653 sage: n = ZZ.random_element(1,5)
1654 sage: B = RealSymmetricEJA._denormalized_basis(n)
1655 sage: all( M.is_symmetric() for M in B)
1659 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1663 for j
in range(i
+1):
1664 Eij
= matrix(ZZ
, n
, lambda k
,l
: k
==i
and l
==j
)
1668 Sij
= Eij
+ Eij
.transpose()
1674 def _max_random_instance_size():
1675 return 4 # Dimension 10
1678 def random_instance(cls
, **kwargs
):
1680 Return a random instance of this type of algebra.
1682 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1683 return cls(n
, **kwargs
)
1685 def __init__(self
, n
, **kwargs
):
1686 # We know this is a valid EJA, but will double-check
1687 # if the user passes check_axioms=True.
1688 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
1690 super(RealSymmetricEJA
, self
).__init
__(self
._denormalized
_basis
(n
),
1691 self
.jordan_product
,
1692 self
.trace_inner_product
,
1695 # TODO: this could be factored out somehow, but is left here
1696 # because the MatrixEuclideanJordanAlgebra is not presently
1697 # a subclass of the FDEJA class that defines rank() and one().
1698 self
.rank
.set_cache(n
)
1699 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
1700 self
.one
.set_cache(self(idV
))
1704 class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1706 def dimension_over_reals():
1710 def real_embed(cls
,M
):
1712 Embed the n-by-n complex matrix ``M`` into the space of real
1713 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1714 bi` to the block matrix ``[[a,b],[-b,a]]``.
1718 sage: from mjo.eja.eja_algebra import \
1719 ....: ComplexMatrixEuclideanJordanAlgebra
1723 sage: F = QuadraticField(-1, 'I')
1724 sage: x1 = F(4 - 2*i)
1725 sage: x2 = F(1 + 2*i)
1728 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1729 sage: ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1738 Embedding is a homomorphism (isomorphism, in fact)::
1740 sage: set_random_seed()
1741 sage: n = ZZ.random_element(3)
1742 sage: F = QuadraticField(-1, 'I')
1743 sage: X = random_matrix(F, n)
1744 sage: Y = random_matrix(F, n)
1745 sage: Xe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X)
1746 sage: Ye = ComplexMatrixEuclideanJordanAlgebra.real_embed(Y)
1747 sage: XYe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1752 super(ComplexMatrixEuclideanJordanAlgebra
,cls
).real_embed(M
)
1755 # We don't need any adjoined elements...
1756 field
= M
.base_ring().base_ring()
1760 a
= z
.list()[0] # real part, I guess
1761 b
= z
.list()[1] # imag part, I guess
1762 blocks
.append(matrix(field
, 2, [[a
,b
],[-b
,a
]]))
1764 return matrix
.block(field
, n
, blocks
)
1768 def real_unembed(cls
,M
):
1770 The inverse of _embed_complex_matrix().
1774 sage: from mjo.eja.eja_algebra import \
1775 ....: ComplexMatrixEuclideanJordanAlgebra
1779 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1780 ....: [-2, 1, -4, 3],
1781 ....: [ 9, 10, 11, 12],
1782 ....: [-10, 9, -12, 11] ])
1783 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(A)
1785 [ 10*I + 9 12*I + 11]
1789 Unembedding is the inverse of embedding::
1791 sage: set_random_seed()
1792 sage: F = QuadraticField(-1, 'I')
1793 sage: M = random_matrix(F, 3)
1794 sage: Me = ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1795 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1799 super(ComplexMatrixEuclideanJordanAlgebra
,cls
).real_unembed(M
)
1801 d
= cls
.dimension_over_reals()
1803 # If "M" was normalized, its base ring might have roots
1804 # adjoined and they can stick around after unembedding.
1805 field
= M
.base_ring()
1806 R
= PolynomialRing(field
, 'z')
1809 # Sage doesn't know how to adjoin the complex "i" (the root of
1810 # x^2 + 1) to a field in a general way. Here, we just enumerate
1811 # all of the cases that I have cared to support so far.
1813 # Sage doesn't know how to embed AA into QQbar, i.e. how
1814 # to adjoin sqrt(-1) to AA.
1816 elif not field
.is_exact():
1818 F
= field
.complex_field()
1820 # Works for QQ and... maybe some other fields.
1821 F
= field
.extension(z
**2 + 1, 'I', embedding
=CLF(-1).sqrt())
1824 # Go top-left to bottom-right (reading order), converting every
1825 # 2-by-2 block we see to a single complex element.
1827 for k
in range(n
/d
):
1828 for j
in range(n
/d
):
1829 submat
= M
[d
*k
:d
*k
+d
,d
*j
:d
*j
+d
]
1830 if submat
[0,0] != submat
[1,1]:
1831 raise ValueError('bad on-diagonal submatrix')
1832 if submat
[0,1] != -submat
[1,0]:
1833 raise ValueError('bad off-diagonal submatrix')
1834 z
= submat
[0,0] + submat
[0,1]*i
1837 return matrix(F
, n
/d
, elements
)
1840 class ComplexHermitianEJA(ConcreteEuclideanJordanAlgebra
,
1841 ComplexMatrixEuclideanJordanAlgebra
):
1843 The rank-n simple EJA consisting of complex Hermitian n-by-n
1844 matrices over the real numbers, the usual symmetric Jordan product,
1845 and the real-part-of-trace inner product. It has dimension `n^2` over
1850 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1854 In theory, our "field" can be any subfield of the reals::
1856 sage: ComplexHermitianEJA(2, field=RDF)
1857 Euclidean Jordan algebra of dimension 4 over Real Double Field
1858 sage: ComplexHermitianEJA(2, field=RR)
1859 Euclidean Jordan algebra of dimension 4 over Real Field with
1860 53 bits of precision
1864 The dimension of this algebra is `n^2`::
1866 sage: set_random_seed()
1867 sage: n_max = ComplexHermitianEJA._max_random_instance_size()
1868 sage: n = ZZ.random_element(1, n_max)
1869 sage: J = ComplexHermitianEJA(n)
1870 sage: J.dimension() == n^2
1873 The Jordan multiplication is what we think it is::
1875 sage: set_random_seed()
1876 sage: J = ComplexHermitianEJA.random_instance()
1877 sage: x,y = J.random_elements(2)
1878 sage: actual = (x*y).to_matrix()
1879 sage: X = x.to_matrix()
1880 sage: Y = y.to_matrix()
1881 sage: expected = (X*Y + Y*X)/2
1882 sage: actual == expected
1884 sage: J(expected) == x*y
1887 We can change the generator prefix::
1889 sage: ComplexHermitianEJA(2, prefix='z').gens()
1892 We can construct the (trivial) algebra of rank zero::
1894 sage: ComplexHermitianEJA(0)
1895 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1900 def _denormalized_basis(cls
, n
):
1902 Returns a basis for the space of complex Hermitian n-by-n matrices.
1904 Why do we embed these? Basically, because all of numerical linear
1905 algebra assumes that you're working with vectors consisting of `n`
1906 entries from a field and scalars from the same field. There's no way
1907 to tell SageMath that (for example) the vectors contain complex
1908 numbers, while the scalar field is real.
1912 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1916 sage: set_random_seed()
1917 sage: n = ZZ.random_element(1,5)
1918 sage: field = QuadraticField(2, 'sqrt2')
1919 sage: B = ComplexHermitianEJA._denormalized_basis(n)
1920 sage: all( M.is_symmetric() for M in B)
1925 R
= PolynomialRing(field
, 'z')
1927 F
= field
.extension(z
**2 + 1, 'I')
1930 # This is like the symmetric case, but we need to be careful:
1932 # * We want conjugate-symmetry, not just symmetry.
1933 # * The diagonal will (as a result) be real.
1937 for j
in range(i
+1):
1938 Eij
= matrix(F
, n
, lambda k
,l
: k
==i
and l
==j
)
1940 Sij
= cls
.real_embed(Eij
)
1943 # The second one has a minus because it's conjugated.
1944 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
1946 Sij_imag
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
1949 # Since we embedded these, we can drop back to the "field" that we
1950 # started with instead of the complex extension "F".
1951 return tuple( s
.change_ring(field
) for s
in S
)
1954 def __init__(self
, n
, **kwargs
):
1955 # We know this is a valid EJA, but will double-check
1956 # if the user passes check_axioms=True.
1957 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
1959 super(ComplexHermitianEJA
, self
).__init
__(self
._denormalized
_basis
(n
),
1960 self
.jordan_product
,
1961 self
.trace_inner_product
,
1963 # TODO: this could be factored out somehow, but is left here
1964 # because the MatrixEuclideanJordanAlgebra is not presently
1965 # a subclass of the FDEJA class that defines rank() and one().
1966 self
.rank
.set_cache(n
)
1967 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
1968 self
.one
.set_cache(self(idV
))
1971 def _max_random_instance_size():
1972 return 3 # Dimension 9
1975 def random_instance(cls
, **kwargs
):
1977 Return a random instance of this type of algebra.
1979 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1980 return cls(n
, **kwargs
)
1982 class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1984 def dimension_over_reals():
1988 def real_embed(cls
,M
):
1990 Embed the n-by-n quaternion matrix ``M`` into the space of real
1991 matrices of size 4n-by-4n by first sending each quaternion entry `z
1992 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
1993 c+di],[-c + di, a-bi]]`, and then embedding those into a real
1998 sage: from mjo.eja.eja_algebra import \
1999 ....: QuaternionMatrixEuclideanJordanAlgebra
2003 sage: Q = QuaternionAlgebra(QQ,-1,-1)
2004 sage: i,j,k = Q.gens()
2005 sage: x = 1 + 2*i + 3*j + 4*k
2006 sage: M = matrix(Q, 1, [[x]])
2007 sage: QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
2013 Embedding is a homomorphism (isomorphism, in fact)::
2015 sage: set_random_seed()
2016 sage: n = ZZ.random_element(2)
2017 sage: Q = QuaternionAlgebra(QQ,-1,-1)
2018 sage: X = random_matrix(Q, n)
2019 sage: Y = random_matrix(Q, n)
2020 sage: Xe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X)
2021 sage: Ye = QuaternionMatrixEuclideanJordanAlgebra.real_embed(Y)
2022 sage: XYe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X*Y)
2027 super(QuaternionMatrixEuclideanJordanAlgebra
,cls
).real_embed(M
)
2028 quaternions
= M
.base_ring()
2031 F
= QuadraticField(-1, 'I')
2036 t
= z
.coefficient_tuple()
2041 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
2042 [-c
+ d
*i
, a
- b
*i
]])
2043 realM
= ComplexMatrixEuclideanJordanAlgebra
.real_embed(cplxM
)
2044 blocks
.append(realM
)
2046 # We should have real entries by now, so use the realest field
2047 # we've got for the return value.
2048 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
2053 def real_unembed(cls
,M
):
2055 The inverse of _embed_quaternion_matrix().
2059 sage: from mjo.eja.eja_algebra import \
2060 ....: QuaternionMatrixEuclideanJordanAlgebra
2064 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
2065 ....: [-2, 1, -4, 3],
2066 ....: [-3, 4, 1, -2],
2067 ....: [-4, -3, 2, 1]])
2068 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(M)
2069 [1 + 2*i + 3*j + 4*k]
2073 Unembedding is the inverse of embedding::
2075 sage: set_random_seed()
2076 sage: Q = QuaternionAlgebra(QQ, -1, -1)
2077 sage: M = random_matrix(Q, 3)
2078 sage: Me = QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
2079 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
2083 super(QuaternionMatrixEuclideanJordanAlgebra
,cls
).real_unembed(M
)
2085 d
= cls
.dimension_over_reals()
2087 # Use the base ring of the matrix to ensure that its entries can be
2088 # multiplied by elements of the quaternion algebra.
2089 field
= M
.base_ring()
2090 Q
= QuaternionAlgebra(field
,-1,-1)
2093 # Go top-left to bottom-right (reading order), converting every
2094 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
2097 for l
in range(n
/d
):
2098 for m
in range(n
/d
):
2099 submat
= ComplexMatrixEuclideanJordanAlgebra
.real_unembed(
2100 M
[d
*l
:d
*l
+d
,d
*m
:d
*m
+d
] )
2101 if submat
[0,0] != submat
[1,1].conjugate():
2102 raise ValueError('bad on-diagonal submatrix')
2103 if submat
[0,1] != -submat
[1,0].conjugate():
2104 raise ValueError('bad off-diagonal submatrix')
2105 z
= submat
[0,0].real()
2106 z
+= submat
[0,0].imag()*i
2107 z
+= submat
[0,1].real()*j
2108 z
+= submat
[0,1].imag()*k
2111 return matrix(Q
, n
/d
, elements
)
2114 class QuaternionHermitianEJA(ConcreteEuclideanJordanAlgebra
,
2115 QuaternionMatrixEuclideanJordanAlgebra
):
2117 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
2118 matrices, the usual symmetric Jordan product, and the
2119 real-part-of-trace inner product. It has dimension `2n^2 - n` over
2124 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2128 In theory, our "field" can be any subfield of the reals::
2130 sage: QuaternionHermitianEJA(2, field=RDF)
2131 Euclidean Jordan algebra of dimension 6 over Real Double Field
2132 sage: QuaternionHermitianEJA(2, field=RR)
2133 Euclidean Jordan algebra of dimension 6 over Real Field with
2134 53 bits of precision
2138 The dimension of this algebra is `2*n^2 - n`::
2140 sage: set_random_seed()
2141 sage: n_max = QuaternionHermitianEJA._max_random_instance_size()
2142 sage: n = ZZ.random_element(1, n_max)
2143 sage: J = QuaternionHermitianEJA(n)
2144 sage: J.dimension() == 2*(n^2) - n
2147 The Jordan multiplication is what we think it is::
2149 sage: set_random_seed()
2150 sage: J = QuaternionHermitianEJA.random_instance()
2151 sage: x,y = J.random_elements(2)
2152 sage: actual = (x*y).to_matrix()
2153 sage: X = x.to_matrix()
2154 sage: Y = y.to_matrix()
2155 sage: expected = (X*Y + Y*X)/2
2156 sage: actual == expected
2158 sage: J(expected) == x*y
2161 We can change the generator prefix::
2163 sage: QuaternionHermitianEJA(2, prefix='a').gens()
2164 (a0, a1, a2, a3, a4, a5)
2166 We can construct the (trivial) algebra of rank zero::
2168 sage: QuaternionHermitianEJA(0)
2169 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2173 def _denormalized_basis(cls
, n
):
2175 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
2177 Why do we embed these? Basically, because all of numerical
2178 linear algebra assumes that you're working with vectors consisting
2179 of `n` entries from a field and scalars from the same field. There's
2180 no way to tell SageMath that (for example) the vectors contain
2181 complex numbers, while the scalar field is real.
2185 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2189 sage: set_random_seed()
2190 sage: n = ZZ.random_element(1,5)
2191 sage: B = QuaternionHermitianEJA._denormalized_basis(n)
2192 sage: all( M.is_symmetric() for M in B )
2197 Q
= QuaternionAlgebra(QQ
,-1,-1)
2200 # This is like the symmetric case, but we need to be careful:
2202 # * We want conjugate-symmetry, not just symmetry.
2203 # * The diagonal will (as a result) be real.
2207 for j
in range(i
+1):
2208 Eij
= matrix(Q
, n
, lambda k
,l
: k
==i
and l
==j
)
2210 Sij
= cls
.real_embed(Eij
)
2213 # The second, third, and fourth ones have a minus
2214 # because they're conjugated.
2215 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
2217 Sij_I
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
2219 Sij_J
= cls
.real_embed(J
*Eij
- J
*Eij
.transpose())
2221 Sij_K
= cls
.real_embed(K
*Eij
- K
*Eij
.transpose())
2224 # Since we embedded these, we can drop back to the "field" that we
2225 # started with instead of the quaternion algebra "Q".
2226 return tuple( s
.change_ring(field
) for s
in S
)
2229 def __init__(self
, n
, **kwargs
):
2230 # We know this is a valid EJA, but will double-check
2231 # if the user passes check_axioms=True.
2232 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2234 super(QuaternionHermitianEJA
, self
).__init
__(self
._denormalized
_basis
(n
),
2235 self
.jordan_product
,
2236 self
.trace_inner_product
,
2238 # TODO: this could be factored out somehow, but is left here
2239 # because the MatrixEuclideanJordanAlgebra is not presently
2240 # a subclass of the FDEJA class that defines rank() and one().
2241 self
.rank
.set_cache(n
)
2242 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
2243 self
.one
.set_cache(self(idV
))
2247 def _max_random_instance_size():
2249 The maximum rank of a random QuaternionHermitianEJA.
2251 return 2 # Dimension 6
2254 def random_instance(cls
, **kwargs
):
2256 Return a random instance of this type of algebra.
2258 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2259 return cls(n
, **kwargs
)
2262 class HadamardEJA(ConcreteEuclideanJordanAlgebra
):
2264 Return the Euclidean Jordan Algebra corresponding to the set
2265 `R^n` under the Hadamard product.
2267 Note: this is nothing more than the Cartesian product of ``n``
2268 copies of the spin algebra. Once Cartesian product algebras
2269 are implemented, this can go.
2273 sage: from mjo.eja.eja_algebra import HadamardEJA
2277 This multiplication table can be verified by hand::
2279 sage: J = HadamardEJA(3)
2280 sage: e0,e1,e2 = J.gens()
2296 We can change the generator prefix::
2298 sage: HadamardEJA(3, prefix='r').gens()
2302 def __init__(self
, n
, **kwargs
):
2303 def jordan_product(x
,y
):
2305 return P(tuple( xi
*yi
for (xi
,yi
) in zip(x
,y
) ))
2306 def inner_product(x
,y
):
2307 return x
.inner_product(y
)
2309 # New defaults for keyword arguments. Don't orthonormalize
2310 # because our basis is already orthonormal with respect to our
2311 # inner-product. Don't check the axioms, because we know this
2312 # is a valid EJA... but do double-check if the user passes
2313 # check_axioms=True. Note: we DON'T override the "check_field"
2314 # default here, because the user can pass in a field!
2315 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2316 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2319 standard_basis
= FreeModule(ZZ
, n
).basis()
2320 super(HadamardEJA
, self
).__init
__(standard_basis
,
2324 self
.rank
.set_cache(n
)
2327 self
.one
.set_cache( self
.zero() )
2329 self
.one
.set_cache( sum(self
.gens()) )
2332 def _max_random_instance_size():
2334 The maximum dimension of a random HadamardEJA.
2339 def random_instance(cls
, **kwargs
):
2341 Return a random instance of this type of algebra.
2343 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2344 return cls(n
, **kwargs
)
2347 class BilinearFormEJA(ConcreteEuclideanJordanAlgebra
):
2349 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2350 with the half-trace inner product and jordan product ``x*y =
2351 (<Bx,y>,y_bar>, x0*y_bar + y0*x_bar)`` where `B = 1 \times B22` is
2352 a symmetric positive-definite "bilinear form" matrix. Its
2353 dimension is the size of `B`, and it has rank two in dimensions
2354 larger than two. It reduces to the ``JordanSpinEJA`` when `B` is
2355 the identity matrix of order ``n``.
2357 We insist that the one-by-one upper-left identity block of `B` be
2358 passed in as well so that we can be passed a matrix of size zero
2359 to construct a trivial algebra.
2363 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
2364 ....: JordanSpinEJA)
2368 When no bilinear form is specified, the identity matrix is used,
2369 and the resulting algebra is the Jordan spin algebra::
2371 sage: B = matrix.identity(AA,3)
2372 sage: J0 = BilinearFormEJA(B)
2373 sage: J1 = JordanSpinEJA(3)
2374 sage: J0.multiplication_table() == J0.multiplication_table()
2377 An error is raised if the matrix `B` does not correspond to a
2378 positive-definite bilinear form::
2380 sage: B = matrix.random(QQ,2,3)
2381 sage: J = BilinearFormEJA(B)
2382 Traceback (most recent call last):
2384 ValueError: bilinear form is not positive-definite
2385 sage: B = matrix.zero(QQ,3)
2386 sage: J = BilinearFormEJA(B)
2387 Traceback (most recent call last):
2389 ValueError: bilinear form is not positive-definite
2393 We can create a zero-dimensional algebra::
2395 sage: B = matrix.identity(AA,0)
2396 sage: J = BilinearFormEJA(B)
2400 We can check the multiplication condition given in the Jordan, von
2401 Neumann, and Wigner paper (and also discussed on my "On the
2402 symmetry..." paper). Note that this relies heavily on the standard
2403 choice of basis, as does anything utilizing the bilinear form
2404 matrix. We opt not to orthonormalize the basis, because if we
2405 did, we would have to normalize the `s_{i}` in a similar manner::
2407 sage: set_random_seed()
2408 sage: n = ZZ.random_element(5)
2409 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2410 sage: B11 = matrix.identity(QQ,1)
2411 sage: B22 = M.transpose()*M
2412 sage: B = block_matrix(2,2,[ [B11,0 ],
2414 sage: J = BilinearFormEJA(B, orthonormalize=False)
2415 sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
2416 sage: V = J.vector_space()
2417 sage: sis = [ J( V([0] + (M.inverse()*ei).list()).column() )
2418 ....: for ei in eis ]
2419 sage: actual = [ sis[i]*sis[j]
2420 ....: for i in range(n-1)
2421 ....: for j in range(n-1) ]
2422 sage: expected = [ J.one() if i == j else J.zero()
2423 ....: for i in range(n-1)
2424 ....: for j in range(n-1) ]
2425 sage: actual == expected
2428 def __init__(self
, B
, **kwargs
):
2429 if not B
.is_positive_definite():
2430 raise ValueError("bilinear form is not positive-definite")
2432 def inner_product(x
,y
):
2433 return (B
*x
).inner_product(y
)
2435 def jordan_product(x
,y
):
2441 z0
= inner_product(x
,y
)
2442 zbar
= y0
*xbar
+ x0
*ybar
2443 return P((z0
,) + tuple(zbar
))
2445 # We know this is a valid EJA, but will double-check
2446 # if the user passes check_axioms=True.
2447 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2450 standard_basis
= FreeModule(ZZ
, n
).basis()
2451 super(BilinearFormEJA
, self
).__init
__(standard_basis
,
2456 # The rank of this algebra is two, unless we're in a
2457 # one-dimensional ambient space (because the rank is bounded
2458 # by the ambient dimension).
2459 self
.rank
.set_cache(min(n
,2))
2462 self
.one
.set_cache( self
.zero() )
2464 self
.one
.set_cache( self
.monomial(0) )
2467 def _max_random_instance_size():
2469 The maximum dimension of a random BilinearFormEJA.
2474 def random_instance(cls
, **kwargs
):
2476 Return a random instance of this algebra.
2478 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2480 B
= matrix
.identity(ZZ
, n
)
2481 return cls(B
, **kwargs
)
2483 B11
= matrix
.identity(ZZ
, 1)
2484 M
= matrix
.random(ZZ
, n
-1)
2485 I
= matrix
.identity(ZZ
, n
-1)
2487 while alpha
.is_zero():
2488 alpha
= ZZ
.random_element().abs()
2489 B22
= M
.transpose()*M
+ alpha
*I
2491 from sage
.matrix
.special
import block_matrix
2492 B
= block_matrix(2,2, [ [B11
, ZZ(0) ],
2495 return cls(B
, **kwargs
)
2498 class JordanSpinEJA(BilinearFormEJA
):
2500 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2501 with the usual inner product and jordan product ``x*y =
2502 (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2507 sage: from mjo.eja.eja_algebra import JordanSpinEJA
2511 This multiplication table can be verified by hand::
2513 sage: J = JordanSpinEJA(4)
2514 sage: e0,e1,e2,e3 = J.gens()
2530 We can change the generator prefix::
2532 sage: JordanSpinEJA(2, prefix='B').gens()
2537 Ensure that we have the usual inner product on `R^n`::
2539 sage: set_random_seed()
2540 sage: J = JordanSpinEJA.random_instance()
2541 sage: x,y = J.random_elements(2)
2542 sage: actual = x.inner_product(y)
2543 sage: expected = x.to_vector().inner_product(y.to_vector())
2544 sage: actual == expected
2548 def __init__(self
, n
, **kwargs
):
2549 # This is a special case of the BilinearFormEJA with the
2550 # identity matrix as its bilinear form.
2551 B
= matrix
.identity(ZZ
, n
)
2553 # Don't orthonormalize because our basis is already
2554 # orthonormal with respect to our inner-product.
2555 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2557 # But also don't pass check_field=False here, because the user
2558 # can pass in a field!
2559 super(JordanSpinEJA
, self
).__init
__(B
, **kwargs
)
2562 def _max_random_instance_size():
2564 The maximum dimension of a random JordanSpinEJA.
2569 def random_instance(cls
, **kwargs
):
2571 Return a random instance of this type of algebra.
2573 Needed here to override the implementation for ``BilinearFormEJA``.
2575 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2576 return cls(n
, **kwargs
)
2579 class TrivialEJA(ConcreteEuclideanJordanAlgebra
):
2581 The trivial Euclidean Jordan algebra consisting of only a zero element.
2585 sage: from mjo.eja.eja_algebra import TrivialEJA
2589 sage: J = TrivialEJA()
2596 sage: 7*J.one()*12*J.one()
2598 sage: J.one().inner_product(J.one())
2600 sage: J.one().norm()
2602 sage: J.one().subalgebra_generated_by()
2603 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2608 def __init__(self
, **kwargs
):
2609 jordan_product
= lambda x
,y
: x
2610 inner_product
= lambda x
,y
: 0
2613 # New defaults for keyword arguments
2614 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2615 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2617 super(TrivialEJA
, self
).__init
__(basis
,
2621 # The rank is zero using my definition, namely the dimension of the
2622 # largest subalgebra generated by any element.
2623 self
.rank
.set_cache(0)
2624 self
.one
.set_cache( self
.zero() )
2627 def random_instance(cls
, **kwargs
):
2628 # We don't take a "size" argument so the superclass method is
2629 # inappropriate for us.
2630 return cls(**kwargs
)
2632 class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra
):
2634 The external (orthogonal) direct sum of two other Euclidean Jordan
2635 algebras. Essentially the Cartesian product of its two factors.
2636 Every Euclidean Jordan algebra decomposes into an orthogonal
2637 direct sum of simple Euclidean Jordan algebras, so no generality
2638 is lost by providing only this construction.
2642 sage: from mjo.eja.eja_algebra import (random_eja,
2644 ....: RealSymmetricEJA,
2649 sage: J1 = HadamardEJA(2)
2650 sage: J2 = RealSymmetricEJA(3)
2651 sage: J = DirectSumEJA(J1,J2)
2659 The external direct sum construction is only valid when the two factors
2660 have the same base ring; an error is raised otherwise::
2662 sage: set_random_seed()
2663 sage: J1 = random_eja(field=AA)
2664 sage: J2 = random_eja(field=QQ,orthonormalize=False)
2665 sage: J = DirectSumEJA(J1,J2)
2666 Traceback (most recent call last):
2668 ValueError: algebras must share the same base field
2671 def __init__(self
, J1
, J2
, **kwargs
):
2672 if J1
.base_ring() != J2
.base_ring():
2673 raise ValueError("algebras must share the same base field")
2674 field
= J1
.base_ring()
2676 self
._factors
= (J1
, J2
)
2680 V
= VectorSpace(field
, n
)
2681 mult_table
= [ [ V
.zero() for j
in range(i
+1) ]
2684 for j
in range(i
+1):
2685 p
= (J1
.monomial(i
)*J1
.monomial(j
)).to_vector()
2686 mult_table
[i
][j
] = V(p
.list() + [field
.zero()]*n2
)
2689 for j
in range(i
+1):
2690 p
= (J2
.monomial(i
)*J2
.monomial(j
)).to_vector()
2691 mult_table
[n1
+i
][n1
+j
] = V([field
.zero()]*n1
+ p
.list())
2693 # TODO: build the IP table here from the two constituent IP
2694 # matrices (it'll be block diagonal, I think).
2695 ip_table
= [ [ field
.zero() for j
in range(i
+1) ]
2697 super(DirectSumEJA
, self
).__init
__(field
,
2702 self
.rank
.set_cache(J1
.rank() + J2
.rank())
2707 Return the pair of this algebra's factors.
2711 sage: from mjo.eja.eja_algebra import (HadamardEJA,
2712 ....: JordanSpinEJA,
2717 sage: J1 = HadamardEJA(2, field=QQ)
2718 sage: J2 = JordanSpinEJA(3, field=QQ)
2719 sage: J = DirectSumEJA(J1,J2)
2721 (Euclidean Jordan algebra of dimension 2 over Rational Field,
2722 Euclidean Jordan algebra of dimension 3 over Rational Field)
2725 return self
._factors
2727 def projections(self
):
2729 Return a pair of projections onto this algebra's factors.
2733 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
2734 ....: ComplexHermitianEJA,
2739 sage: J1 = JordanSpinEJA(2)
2740 sage: J2 = ComplexHermitianEJA(2)
2741 sage: J = DirectSumEJA(J1,J2)
2742 sage: (pi_left, pi_right) = J.projections()
2743 sage: J.one().to_vector()
2745 sage: pi_left(J.one()).to_vector()
2747 sage: pi_right(J.one()).to_vector()
2751 (J1
,J2
) = self
.factors()
2754 V_basis
= self
.vector_space().basis()
2755 # Need to specify the dimensions explicitly so that we don't
2756 # wind up with a zero-by-zero matrix when we want e.g. a
2757 # zero-by-two matrix (important for composing things).
2758 P1
= matrix(self
.base_ring(), m
, m
+n
, V_basis
[:m
])
2759 P2
= matrix(self
.base_ring(), n
, m
+n
, V_basis
[m
:])
2760 pi_left
= FiniteDimensionalEuclideanJordanAlgebraOperator(self
,J1
,P1
)
2761 pi_right
= FiniteDimensionalEuclideanJordanAlgebraOperator(self
,J2
,P2
)
2762 return (pi_left
, pi_right
)
2764 def inclusions(self
):
2766 Return the pair of inclusion maps from our factors into us.
2770 sage: from mjo.eja.eja_algebra import (random_eja,
2771 ....: JordanSpinEJA,
2772 ....: RealSymmetricEJA,
2777 sage: J1 = JordanSpinEJA(3)
2778 sage: J2 = RealSymmetricEJA(2)
2779 sage: J = DirectSumEJA(J1,J2)
2780 sage: (iota_left, iota_right) = J.inclusions()
2781 sage: iota_left(J1.zero()) == J.zero()
2783 sage: iota_right(J2.zero()) == J.zero()
2785 sage: J1.one().to_vector()
2787 sage: iota_left(J1.one()).to_vector()
2789 sage: J2.one().to_vector()
2791 sage: iota_right(J2.one()).to_vector()
2793 sage: J.one().to_vector()
2798 Composing a projection with the corresponding inclusion should
2799 produce the identity map, and mismatching them should produce
2802 sage: set_random_seed()
2803 sage: J1 = random_eja()
2804 sage: J2 = random_eja()
2805 sage: J = DirectSumEJA(J1,J2)
2806 sage: (iota_left, iota_right) = J.inclusions()
2807 sage: (pi_left, pi_right) = J.projections()
2808 sage: pi_left*iota_left == J1.one().operator()
2810 sage: pi_right*iota_right == J2.one().operator()
2812 sage: (pi_left*iota_right).is_zero()
2814 sage: (pi_right*iota_left).is_zero()
2818 (J1
,J2
) = self
.factors()
2821 V_basis
= self
.vector_space().basis()
2822 # Need to specify the dimensions explicitly so that we don't
2823 # wind up with a zero-by-zero matrix when we want e.g. a
2824 # two-by-zero matrix (important for composing things).
2825 I1
= matrix
.column(self
.base_ring(), m
, m
+n
, V_basis
[:m
])
2826 I2
= matrix
.column(self
.base_ring(), n
, m
+n
, V_basis
[m
:])
2827 iota_left
= FiniteDimensionalEuclideanJordanAlgebraOperator(J1
,self
,I1
)
2828 iota_right
= FiniteDimensionalEuclideanJordanAlgebraOperator(J2
,self
,I2
)
2829 return (iota_left
, iota_right
)
2831 def inner_product(self
, x
, y
):
2833 The standard Cartesian inner-product.
2835 We project ``x`` and ``y`` onto our factors, and add up the
2836 inner-products from the subalgebras.
2841 sage: from mjo.eja.eja_algebra import (HadamardEJA,
2842 ....: QuaternionHermitianEJA,
2847 sage: J1 = HadamardEJA(3,field=QQ)
2848 sage: J2 = QuaternionHermitianEJA(2,field=QQ,orthonormalize=False)
2849 sage: J = DirectSumEJA(J1,J2)
2854 sage: x1.inner_product(x2)
2856 sage: y1.inner_product(y2)
2858 sage: J.one().inner_product(J.one())
2862 (pi_left
, pi_right
) = self
.projections()
2868 return (x1
.inner_product(y1
) + x2
.inner_product(y2
))
2872 random_eja
= ConcreteEuclideanJordanAlgebra
.random_instance