2 Euclidean Jordan Algebras. These are formally-real Jordan Algebras;
3 specifically those where u^2 + v^2 = 0 implies that u = v = 0. They
4 are used in optimization, and have some additional nice methods beyond
5 what can be supported in a general Jordan Algebra.
10 sage: from mjo.eja.eja_algebra import random_eja
15 Euclidean Jordan algebra of dimension...
19 from itertools
import repeat
21 from sage
.algebras
.quatalg
.quaternion_algebra
import QuaternionAlgebra
22 from sage
.categories
.magmatic_algebras
import MagmaticAlgebras
23 from sage
.combinat
.free_module
import CombinatorialFreeModule
24 from sage
.matrix
.constructor
import matrix
25 from sage
.matrix
.matrix_space
import MatrixSpace
26 from sage
.misc
.cachefunc
import cached_method
27 from sage
.misc
.table
import table
28 from sage
.modules
.free_module
import FreeModule
, VectorSpace
29 from sage
.rings
.all
import (ZZ
, QQ
, AA
, QQbar
, RR
, RLF
, CLF
,
32 from mjo
.eja
.eja_element
import FiniteDimensionalEuclideanJordanAlgebraElement
33 from mjo
.eja
.eja_operator
import FiniteDimensionalEuclideanJordanAlgebraOperator
34 from mjo
.eja
.eja_utils
import _mat2vec
36 class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule
):
38 The lowest-level class for representing a Euclidean Jordan algebra.
40 def _coerce_map_from_base_ring(self
):
42 Disable the map from the base ring into the algebra.
44 Performing a nonsense conversion like this automatically
45 is counterpedagogical. The fallback is to try the usual
46 element constructor, which should also fail.
50 sage: from mjo.eja.eja_algebra import random_eja
54 sage: set_random_seed()
55 sage: J = random_eja()
57 Traceback (most recent call last):
59 ValueError: not an element of this algebra
76 * field -- the scalar field for this algebra (must be real)
78 * multiplication_table -- the multiplication table for this
79 algebra's implicit basis. Only the lower-triangular portion
80 of the table is used, since the multiplication is assumed
85 sage: from mjo.eja.eja_algebra import (
86 ....: FiniteDimensionalEuclideanJordanAlgebra,
92 By definition, Jordan multiplication commutes::
94 sage: set_random_seed()
95 sage: J = random_eja()
96 sage: x,y = J.random_elements(2)
100 An error is raised if the Jordan product is not commutative::
102 sage: JP = ((1,2),(0,0))
103 sage: IP = ((1,0),(0,1))
104 sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,JP,IP)
105 Traceback (most recent call last):
107 ValueError: Jordan product is not commutative
109 An error is raised if the inner-product is not commutative::
111 sage: JP = ((1,0),(0,1))
112 sage: IP = ((1,2),(0,0))
113 sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,JP,IP)
114 Traceback (most recent call last):
116 ValueError: inner-product is not commutative
120 The ``field`` we're given must be real with ``check_field=True``::
122 sage: JordanSpinEJA(2,QQbar)
123 Traceback (most recent call last):
125 ValueError: scalar field is not real
127 The multiplication table must be square with ``check_axioms=True``::
129 sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,((),()),((1,),))
130 Traceback (most recent call last):
132 ValueError: multiplication table is not square
134 The multiplication and inner-product tables must be the same
135 size (and in particular, the inner-product table must also be
136 square) with ``check_axioms=True``::
138 sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,((1,),),(()))
139 Traceback (most recent call last):
141 ValueError: multiplication and inner-product tables are
143 sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,((1,),),((1,2),))
144 Traceback (most recent call last):
146 ValueError: multiplication and inner-product tables are
151 if not field
.is_subring(RR
):
152 # Note: this does return true for the real algebraic
153 # field, the rationals, and any quadratic field where
154 # we've specified a real embedding.
155 raise ValueError("scalar field is not real")
158 # The multiplication and inner-product tables should be square
159 # if the user wants us to verify them. And we verify them as
160 # soon as possible, because we want to exploit their symmetry.
161 n
= len(multiplication_table
)
163 if not all( len(l
) == n
for l
in multiplication_table
):
164 raise ValueError("multiplication table is not square")
166 # If the multiplication table is square, we can check if
167 # the inner-product table is square by comparing it to the
168 # multiplication table's dimensions.
169 msg
= "multiplication and inner-product tables are different sizes"
170 if not len(inner_product_table
) == n
:
171 raise ValueError(msg
)
173 if not all( len(l
) == n
for l
in inner_product_table
):
174 raise ValueError(msg
)
176 # Check commutativity of the Jordan product (symmetry of
177 # the multiplication table) and the commutativity of the
178 # inner-product (symmetry of the inner-product table)
179 # first if we're going to check them at all.. This has to
180 # be done before we define product_on_basis(), because
181 # that method assumes that self._multiplication_table is
182 # symmetric. And it has to be done before we build
183 # self._inner_product_matrix, because the process used to
184 # construct it assumes symmetry as well.
185 if not all( multiplication_table
[j
][i
]
186 == multiplication_table
[i
][j
]
188 for j
in range(i
+1) ):
189 raise ValueError("Jordan product is not commutative")
191 if not all( inner_product_table
[j
][i
]
192 == inner_product_table
[i
][j
]
194 for j
in range(i
+1) ):
195 raise ValueError("inner-product is not commutative")
197 self
._matrix
_basis
= matrix_basis
200 category
= MagmaticAlgebras(field
).FiniteDimensional()
201 category
= category
.WithBasis().Unital()
203 fda
= super(FiniteDimensionalEuclideanJordanAlgebra
, self
)
208 self
.print_options(bracket
='')
210 # The multiplication table we're given is necessarily in terms
211 # of vectors, because we don't have an algebra yet for
212 # anything to be an element of. However, it's faster in the
213 # long run to have the multiplication table be in terms of
214 # algebra elements. We do this after calling the superclass
215 # constructor so that from_vector() knows what to do.
217 # Note: we take advantage of symmetry here, and only store
218 # the lower-triangular portion of the table.
219 self
._multiplication
_table
= [ [ self
.vector_space().zero()
220 for j
in range(i
+1) ]
225 elt
= self
.from_vector(multiplication_table
[i
][j
])
226 self
._multiplication
_table
[i
][j
] = elt
228 self
._multiplication
_table
= tuple(map(tuple, self
._multiplication
_table
))
230 # Save our inner product as a matrix, since the efficiency of
231 # matrix multiplication will usually outweigh the fact that we
232 # have to store a redundant upper- or lower-triangular part.
233 # Pre-cache the fact that these are Hermitian (real symmetric,
234 # in fact) in case some e.g. matrix multiplication routine can
235 # take advantage of it.
236 ip_matrix_constructor
= lambda i
,j
: inner_product_table
[i
][j
] if j
<= i
else inner_product_table
[j
][i
]
237 self
._inner
_product
_matrix
= matrix(field
, n
, ip_matrix_constructor
)
238 self
._inner
_product
_matrix
._cache
= {'hermitian': True}
239 self
._inner
_product
_matrix
.set_immutable()
242 if not self
._is
_jordanian
():
243 raise ValueError("Jordan identity does not hold")
244 if not self
._inner
_product
_is
_associative
():
245 raise ValueError("inner product is not associative")
247 def _element_constructor_(self
, elt
):
249 Construct an element of this algebra from its vector or matrix
252 This gets called only after the parent element _call_ method
253 fails to find a coercion for the argument.
257 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
259 ....: RealSymmetricEJA)
263 The identity in `S^n` is converted to the identity in the EJA::
265 sage: J = RealSymmetricEJA(3)
266 sage: I = matrix.identity(QQ,3)
267 sage: J(I) == J.one()
270 This skew-symmetric matrix can't be represented in the EJA::
272 sage: J = RealSymmetricEJA(3)
273 sage: A = matrix(QQ,3, lambda i,j: i-j)
275 Traceback (most recent call last):
277 ValueError: not an element of this algebra
281 Ensure that we can convert any element of the two non-matrix
282 simple algebras (whose matrix representations are columns)
283 back and forth faithfully::
285 sage: set_random_seed()
286 sage: J = HadamardEJA.random_instance()
287 sage: x = J.random_element()
288 sage: J(x.to_vector().column()) == x
290 sage: J = JordanSpinEJA.random_instance()
291 sage: x = J.random_element()
292 sage: J(x.to_vector().column()) == x
295 msg
= "not an element of this algebra"
297 # The superclass implementation of random_element()
298 # needs to be able to coerce "0" into the algebra.
300 elif elt
in self
.base_ring():
301 # Ensure that no base ring -> algebra coercion is performed
302 # by this method. There's some stupidity in sage that would
303 # otherwise propagate to this method; for example, sage thinks
304 # that the integer 3 belongs to the space of 2-by-2 matrices.
305 raise ValueError(msg
)
307 if elt
not in self
.matrix_space():
308 raise ValueError(msg
)
310 # Thanks for nothing! Matrix spaces aren't vector spaces in
311 # Sage, so we have to figure out its matrix-basis coordinates
312 # ourselves. We use the basis space's ring instead of the
313 # element's ring because the basis space might be an algebraic
314 # closure whereas the base ring of the 3-by-3 identity matrix
315 # could be QQ instead of QQbar.
316 V
= VectorSpace(self
.base_ring(), elt
.nrows()*elt
.ncols())
317 W
= V
.span_of_basis( _mat2vec(s
) for s
in self
.matrix_basis() )
320 coords
= W
.coordinate_vector(_mat2vec(elt
))
321 except ArithmeticError: # vector is not in free module
322 raise ValueError(msg
)
324 return self
.from_vector(coords
)
328 Return a string representation of ``self``.
332 sage: from mjo.eja.eja_algebra import JordanSpinEJA
336 Ensure that it says what we think it says::
338 sage: JordanSpinEJA(2, field=AA)
339 Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
340 sage: JordanSpinEJA(3, field=RDF)
341 Euclidean Jordan algebra of dimension 3 over Real Double Field
344 fmt
= "Euclidean Jordan algebra of dimension {} over {}"
345 return fmt
.format(self
.dimension(), self
.base_ring())
347 def product_on_basis(self
, i
, j
):
348 # We only stored the lower-triangular portion of the
349 # multiplication table.
351 return self
._multiplication
_table
[i
][j
]
353 return self
._multiplication
_table
[j
][i
]
355 def _is_commutative(self
):
357 Whether or not this algebra's multiplication table is commutative.
359 This method should of course always return ``True``, unless
360 this algebra was constructed with ``check_axioms=False`` and
361 passed an invalid multiplication table.
363 return all( self
.product_on_basis(i
,j
) == self
.product_on_basis(i
,j
)
364 for i
in range(self
.dimension())
365 for j
in range(self
.dimension()) )
367 def _is_jordanian(self
):
369 Whether or not this algebra's multiplication table respects the
370 Jordan identity `(x^{2})(xy) = x(x^{2}y)`.
372 We only check one arrangement of `x` and `y`, so for a
373 ``True`` result to be truly true, you should also check
374 :meth:`_is_commutative`. This method should of course always
375 return ``True``, unless this algebra was constructed with
376 ``check_axioms=False`` and passed an invalid multiplication table.
378 return all( (self
.monomial(i
)**2)*(self
.monomial(i
)*self
.monomial(j
))
380 (self
.monomial(i
))*((self
.monomial(i
)**2)*self
.monomial(j
))
381 for i
in range(self
.dimension())
382 for j
in range(self
.dimension()) )
384 def _inner_product_is_associative(self
):
386 Return whether or not this algebra's inner product `B` is
387 associative; that is, whether or not `B(xy,z) = B(x,yz)`.
389 This method should of course always return ``True``, unless
390 this algebra was constructed with ``check_axioms=False`` and
391 passed an invalid multiplication table.
394 # Used to check whether or not something is zero in an inexact
395 # ring. This number is sufficient to allow the construction of
396 # QuaternionHermitianEJA(2, RDF) with check_axioms=True.
399 for i
in range(self
.dimension()):
400 for j
in range(self
.dimension()):
401 for k
in range(self
.dimension()):
405 diff
= (x
*y
).inner_product(z
) - x
.inner_product(y
*z
)
407 if self
.base_ring().is_exact():
411 if diff
.abs() > epsilon
:
417 def characteristic_polynomial_of(self
):
419 Return the algebra's "characteristic polynomial of" function,
420 which is itself a multivariate polynomial that, when evaluated
421 at the coordinates of some algebra element, returns that
422 element's characteristic polynomial.
424 The resulting polynomial has `n+1` variables, where `n` is the
425 dimension of this algebra. The first `n` variables correspond to
426 the coordinates of an algebra element: when evaluated at the
427 coordinates of an algebra element with respect to a certain
428 basis, the result is a univariate polynomial (in the one
429 remaining variable ``t``), namely the characteristic polynomial
434 sage: from mjo.eja.eja_algebra import JordanSpinEJA, TrivialEJA
438 The characteristic polynomial in the spin algebra is given in
439 Alizadeh, Example 11.11::
441 sage: J = JordanSpinEJA(3)
442 sage: p = J.characteristic_polynomial_of(); p
443 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
444 sage: xvec = J.one().to_vector()
448 By definition, the characteristic polynomial is a monic
449 degree-zero polynomial in a rank-zero algebra. Note that
450 Cayley-Hamilton is indeed satisfied since the polynomial
451 ``1`` evaluates to the identity element of the algebra on
454 sage: J = TrivialEJA()
455 sage: J.characteristic_polynomial_of()
462 # The list of coefficient polynomials a_0, a_1, a_2, ..., a_(r-1).
463 a
= self
._charpoly
_coefficients
()
465 # We go to a bit of trouble here to reorder the
466 # indeterminates, so that it's easier to evaluate the
467 # characteristic polynomial at x's coordinates and get back
468 # something in terms of t, which is what we want.
469 S
= PolynomialRing(self
.base_ring(),'t')
473 S
= PolynomialRing(S
, R
.variable_names())
476 return (t
**r
+ sum( a
[k
]*(t
**k
) for k
in range(r
) ))
478 def coordinate_polynomial_ring(self
):
480 The multivariate polynomial ring in which this algebra's
481 :meth:`characteristic_polynomial_of` lives.
485 sage: from mjo.eja.eja_algebra import (HadamardEJA,
486 ....: RealSymmetricEJA)
490 sage: J = HadamardEJA(2)
491 sage: J.coordinate_polynomial_ring()
492 Multivariate Polynomial Ring in X1, X2...
493 sage: J = RealSymmetricEJA(3,QQ,orthonormalize=False)
494 sage: J.coordinate_polynomial_ring()
495 Multivariate Polynomial Ring in X1, X2, X3, X4, X5, X6...
498 var_names
= tuple( "X%d" % z
for z
in range(1, self
.dimension()+1) )
499 return PolynomialRing(self
.base_ring(), var_names
)
501 def inner_product(self
, x
, y
):
503 The inner product associated with this Euclidean Jordan algebra.
505 Defaults to the trace inner product, but can be overridden by
506 subclasses if they are sure that the necessary properties are
511 sage: from mjo.eja.eja_algebra import (random_eja,
513 ....: BilinearFormEJA)
517 Our inner product is "associative," which means the following for
518 a symmetric bilinear form::
520 sage: set_random_seed()
521 sage: J = random_eja()
522 sage: x,y,z = J.random_elements(3)
523 sage: (x*y).inner_product(z) == y.inner_product(x*z)
528 Ensure that this is the usual inner product for the algebras
531 sage: set_random_seed()
532 sage: J = HadamardEJA.random_instance()
533 sage: x,y = J.random_elements(2)
534 sage: actual = x.inner_product(y)
535 sage: expected = x.to_vector().inner_product(y.to_vector())
536 sage: actual == expected
539 Ensure that this is one-half of the trace inner-product in a
540 BilinearFormEJA that isn't just the reals (when ``n`` isn't
541 one). This is in Faraut and Koranyi, and also my "On the
544 sage: set_random_seed()
545 sage: J = BilinearFormEJA.random_instance()
546 sage: n = J.dimension()
547 sage: x = J.random_element()
548 sage: y = J.random_element()
549 sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
552 B
= self
._inner
_product
_matrix
553 return (B
*x
.to_vector()).inner_product(y
.to_vector())
556 def is_trivial(self
):
558 Return whether or not this algebra is trivial.
560 A trivial algebra contains only the zero element.
564 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
569 sage: J = ComplexHermitianEJA(3)
575 sage: J = TrivialEJA()
580 return self
.dimension() == 0
583 def multiplication_table(self
):
585 Return a visual representation of this algebra's multiplication
586 table (on basis elements).
590 sage: from mjo.eja.eja_algebra import JordanSpinEJA
594 sage: J = JordanSpinEJA(4)
595 sage: J.multiplication_table()
596 +----++----+----+----+----+
597 | * || e0 | e1 | e2 | e3 |
598 +====++====+====+====+====+
599 | e0 || e0 | e1 | e2 | e3 |
600 +----++----+----+----+----+
601 | e1 || e1 | e0 | 0 | 0 |
602 +----++----+----+----+----+
603 | e2 || e2 | 0 | e0 | 0 |
604 +----++----+----+----+----+
605 | e3 || e3 | 0 | 0 | e0 |
606 +----++----+----+----+----+
610 M
= [ [ self
.zero() for j
in range(n
) ]
614 M
[i
][j
] = self
._multiplication
_table
[i
][j
]
618 # Prepend the left "header" column entry Can't do this in
619 # the loop because it messes up the symmetry.
620 M
[i
] = [self
.monomial(i
)] + M
[i
]
622 # Prepend the header row.
623 M
= [["*"] + list(self
.gens())] + M
624 return table(M
, header_row
=True, header_column
=True, frame
=True)
627 def matrix_basis(self
):
629 Return an (often more natural) representation of this algebras
630 basis as an ordered tuple of matrices.
632 Every finite-dimensional Euclidean Jordan Algebra is a, up to
633 Jordan isomorphism, a direct sum of five simple
634 algebras---four of which comprise Hermitian matrices. And the
635 last type of algebra can of course be thought of as `n`-by-`1`
636 column matrices (ambiguusly called column vectors) to avoid
637 special cases. As a result, matrices (and column vectors) are
638 a natural representation format for Euclidean Jordan algebra
641 But, when we construct an algebra from a basis of matrices,
642 those matrix representations are lost in favor of coordinate
643 vectors *with respect to* that basis. We could eventually
644 convert back if we tried hard enough, but having the original
645 representations handy is valuable enough that we simply store
646 them and return them from this method.
648 Why implement this for non-matrix algebras? Avoiding special
649 cases for the :class:`BilinearFormEJA` pays with simplicity in
650 its own right. But mainly, we would like to be able to assume
651 that elements of a :class:`DirectSumEJA` can be displayed
652 nicely, without having to have special classes for direct sums
653 one of whose components was a matrix algebra.
657 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
658 ....: RealSymmetricEJA)
662 sage: J = RealSymmetricEJA(2)
664 Finite family {0: e0, 1: e1, 2: e2}
665 sage: J.matrix_basis()
667 [1 0] [ 0 0.7071067811865475?] [0 0]
668 [0 0], [0.7071067811865475? 0], [0 1]
673 sage: J = JordanSpinEJA(2)
675 Finite family {0: e0, 1: e1}
676 sage: J.matrix_basis()
682 if self
._matrix
_basis
is None:
683 M
= self
.matrix_space()
684 return tuple( M(b
.to_vector()) for b
in self
.basis() )
686 return self
._matrix
_basis
689 def matrix_space(self
):
691 Return the matrix space in which this algebra's elements live, if
692 we think of them as matrices (including column vectors of the
695 Generally this will be an `n`-by-`1` column-vector space,
696 except when the algebra is trivial. There it's `n`-by-`n`
697 (where `n` is zero), to ensure that two elements of the matrix
698 space (empty matrices) can be multiplied.
700 Matrix algebras override this with something more useful.
702 if self
.is_trivial():
703 return MatrixSpace(self
.base_ring(), 0)
704 elif self
._matrix
_basis
is None or len(self
._matrix
_basis
) == 0:
705 return MatrixSpace(self
.base_ring(), self
.dimension(), 1)
707 return self
._matrix
_basis
[0].matrix_space()
713 Return the unit element of this algebra.
717 sage: from mjo.eja.eja_algebra import (HadamardEJA,
722 sage: J = HadamardEJA(5)
724 e0 + e1 + e2 + e3 + e4
728 The identity element acts like the identity::
730 sage: set_random_seed()
731 sage: J = random_eja()
732 sage: x = J.random_element()
733 sage: J.one()*x == x and x*J.one() == x
736 The matrix of the unit element's operator is the identity::
738 sage: set_random_seed()
739 sage: J = random_eja()
740 sage: actual = J.one().operator().matrix()
741 sage: expected = matrix.identity(J.base_ring(), J.dimension())
742 sage: actual == expected
745 Ensure that the cached unit element (often precomputed by
746 hand) agrees with the computed one::
748 sage: set_random_seed()
749 sage: J = random_eja()
750 sage: cached = J.one()
751 sage: J.one.clear_cache()
752 sage: J.one() == cached
756 # We can brute-force compute the matrices of the operators
757 # that correspond to the basis elements of this algebra.
758 # If some linear combination of those basis elements is the
759 # algebra identity, then the same linear combination of
760 # their matrices has to be the identity matrix.
762 # Of course, matrices aren't vectors in sage, so we have to
763 # appeal to the "long vectors" isometry.
764 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
766 # Now we use basic linear algebra to find the coefficients,
767 # of the matrices-as-vectors-linear-combination, which should
768 # work for the original algebra basis too.
769 A
= matrix(self
.base_ring(), oper_vecs
)
771 # We used the isometry on the left-hand side already, but we
772 # still need to do it for the right-hand side. Recall that we
773 # wanted something that summed to the identity matrix.
774 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
776 # Now if there's an identity element in the algebra, this
777 # should work. We solve on the left to avoid having to
778 # transpose the matrix "A".
779 return self
.from_vector(A
.solve_left(b
))
782 def peirce_decomposition(self
, c
):
784 The Peirce decomposition of this algebra relative to the
787 In the future, this can be extended to a complete system of
788 orthogonal idempotents.
792 - ``c`` -- an idempotent of this algebra.
796 A triple (J0, J5, J1) containing two subalgebras and one subspace
799 - ``J0`` -- the algebra on the eigenspace of ``c.operator()``
800 corresponding to the eigenvalue zero.
802 - ``J5`` -- the eigenspace (NOT a subalgebra) of ``c.operator()``
803 corresponding to the eigenvalue one-half.
805 - ``J1`` -- the algebra on the eigenspace of ``c.operator()``
806 corresponding to the eigenvalue one.
808 These are the only possible eigenspaces for that operator, and this
809 algebra is a direct sum of them. The spaces ``J0`` and ``J1`` are
810 orthogonal, and are subalgebras of this algebra with the appropriate
815 sage: from mjo.eja.eja_algebra import random_eja, RealSymmetricEJA
819 The canonical example comes from the symmetric matrices, which
820 decompose into diagonal and off-diagonal parts::
822 sage: J = RealSymmetricEJA(3)
823 sage: C = matrix(QQ, [ [1,0,0],
827 sage: J0,J5,J1 = J.peirce_decomposition(c)
829 Euclidean Jordan algebra of dimension 1...
831 Vector space of degree 6 and dimension 2...
833 Euclidean Jordan algebra of dimension 3...
834 sage: J0.one().to_matrix()
838 sage: orig_df = AA.options.display_format
839 sage: AA.options.display_format = 'radical'
840 sage: J.from_vector(J5.basis()[0]).to_matrix()
844 sage: J.from_vector(J5.basis()[1]).to_matrix()
848 sage: AA.options.display_format = orig_df
849 sage: J1.one().to_matrix()
856 Every algebra decomposes trivially with respect to its identity
859 sage: set_random_seed()
860 sage: J = random_eja()
861 sage: J0,J5,J1 = J.peirce_decomposition(J.one())
862 sage: J0.dimension() == 0 and J5.dimension() == 0
864 sage: J1.superalgebra() == J and J1.dimension() == J.dimension()
867 The decomposition is into eigenspaces, and its components are
868 therefore necessarily orthogonal. Moreover, the identity
869 elements in the two subalgebras are the projections onto their
870 respective subspaces of the superalgebra's identity element::
872 sage: set_random_seed()
873 sage: J = random_eja()
874 sage: x = J.random_element()
875 sage: if not J.is_trivial():
876 ....: while x.is_nilpotent():
877 ....: x = J.random_element()
878 sage: c = x.subalgebra_idempotent()
879 sage: J0,J5,J1 = J.peirce_decomposition(c)
881 sage: for (w,y,z) in zip(J0.basis(), J5.basis(), J1.basis()):
882 ....: w = w.superalgebra_element()
883 ....: y = J.from_vector(y)
884 ....: z = z.superalgebra_element()
885 ....: ipsum += w.inner_product(y).abs()
886 ....: ipsum += w.inner_product(z).abs()
887 ....: ipsum += y.inner_product(z).abs()
890 sage: J1(c) == J1.one()
892 sage: J0(J.one() - c) == J0.one()
896 if not c
.is_idempotent():
897 raise ValueError("element is not idempotent: %s" % c
)
899 from mjo
.eja
.eja_subalgebra
import FiniteDimensionalEuclideanJordanSubalgebra
901 # Default these to what they should be if they turn out to be
902 # trivial, because eigenspaces_left() won't return eigenvalues
903 # corresponding to trivial spaces (e.g. it returns only the
904 # eigenspace corresponding to lambda=1 if you take the
905 # decomposition relative to the identity element).
906 trivial
= FiniteDimensionalEuclideanJordanSubalgebra(self
, ())
907 J0
= trivial
# eigenvalue zero
908 J5
= VectorSpace(self
.base_ring(), 0) # eigenvalue one-half
909 J1
= trivial
# eigenvalue one
911 for (eigval
, eigspace
) in c
.operator().matrix().right_eigenspaces():
912 if eigval
== ~
(self
.base_ring()(2)):
915 gens
= tuple( self
.from_vector(b
) for b
in eigspace
.basis() )
916 subalg
= FiniteDimensionalEuclideanJordanSubalgebra(self
,
924 raise ValueError("unexpected eigenvalue: %s" % eigval
)
929 def random_element(self
, thorough
=False):
931 Return a random element of this algebra.
933 Our algebra superclass method only returns a linear
934 combination of at most two basis elements. We instead
935 want the vector space "random element" method that
936 returns a more diverse selection.
940 - ``thorough`` -- (boolean; default False) whether or not we
941 should generate irrational coefficients for the random
942 element when our base ring is irrational; this slows the
943 algebra operations to a crawl, but any truly random method
947 # For a general base ring... maybe we can trust this to do the
948 # right thing? Unlikely, but.
949 V
= self
.vector_space()
950 v
= V
.random_element()
952 if self
.base_ring() is AA
:
953 # The "random element" method of the algebraic reals is
954 # stupid at the moment, and only returns integers between
955 # -2 and 2, inclusive:
957 # https://trac.sagemath.org/ticket/30875
959 # Instead, we implement our own "random vector" method,
960 # and then coerce that into the algebra. We use the vector
961 # space degree here instead of the dimension because a
962 # subalgebra could (for example) be spanned by only two
963 # vectors, each with five coordinates. We need to
964 # generate all five coordinates.
966 v
*= QQbar
.random_element().real()
968 v
*= QQ
.random_element()
970 return self
.from_vector(V
.coordinate_vector(v
))
972 def random_elements(self
, count
, thorough
=False):
974 Return ``count`` random elements as a tuple.
978 - ``thorough`` -- (boolean; default False) whether or not we
979 should generate irrational coefficients for the random
980 elements when our base ring is irrational; this slows the
981 algebra operations to a crawl, but any truly random method
986 sage: from mjo.eja.eja_algebra import JordanSpinEJA
990 sage: J = JordanSpinEJA(3)
991 sage: x,y,z = J.random_elements(3)
992 sage: all( [ x in J, y in J, z in J ])
994 sage: len( J.random_elements(10) ) == 10
998 return tuple( self
.random_element(thorough
)
999 for idx
in range(count
) )
1003 def _charpoly_coefficients(self
):
1005 The `r` polynomial coefficients of the "characteristic polynomial
1008 n
= self
.dimension()
1009 R
= self
.coordinate_polynomial_ring()
1011 F
= R
.fraction_field()
1014 # From a result in my book, these are the entries of the
1015 # basis representation of L_x.
1016 return sum( vars[k
]*self
.monomial(k
).operator().matrix()[i
,j
]
1019 L_x
= matrix(F
, n
, n
, L_x_i_j
)
1022 if self
.rank
.is_in_cache():
1024 # There's no need to pad the system with redundant
1025 # columns if we *know* they'll be redundant.
1028 # Compute an extra power in case the rank is equal to
1029 # the dimension (otherwise, we would stop at x^(r-1)).
1030 x_powers
= [ (L_x
**k
)*self
.one().to_vector()
1031 for k
in range(n
+1) ]
1032 A
= matrix
.column(F
, x_powers
[:n
])
1033 AE
= A
.extended_echelon_form()
1040 # The theory says that only the first "r" coefficients are
1041 # nonzero, and they actually live in the original polynomial
1042 # ring and not the fraction field. We negate them because
1043 # in the actual characteristic polynomial, they get moved
1044 # to the other side where x^r lives.
1045 return -A_rref
.solve_right(E
*b
).change_ring(R
)[:r
]
1050 Return the rank of this EJA.
1052 This is a cached method because we know the rank a priori for
1053 all of the algebras we can construct. Thus we can avoid the
1054 expensive ``_charpoly_coefficients()`` call unless we truly
1055 need to compute the whole characteristic polynomial.
1059 sage: from mjo.eja.eja_algebra import (HadamardEJA,
1060 ....: JordanSpinEJA,
1061 ....: RealSymmetricEJA,
1062 ....: ComplexHermitianEJA,
1063 ....: QuaternionHermitianEJA,
1068 The rank of the Jordan spin algebra is always two::
1070 sage: JordanSpinEJA(2).rank()
1072 sage: JordanSpinEJA(3).rank()
1074 sage: JordanSpinEJA(4).rank()
1077 The rank of the `n`-by-`n` Hermitian real, complex, or
1078 quaternion matrices is `n`::
1080 sage: RealSymmetricEJA(4).rank()
1082 sage: ComplexHermitianEJA(3).rank()
1084 sage: QuaternionHermitianEJA(2).rank()
1089 Ensure that every EJA that we know how to construct has a
1090 positive integer rank, unless the algebra is trivial in
1091 which case its rank will be zero::
1093 sage: set_random_seed()
1094 sage: J = random_eja()
1098 sage: r > 0 or (r == 0 and J.is_trivial())
1101 Ensure that computing the rank actually works, since the ranks
1102 of all simple algebras are known and will be cached by default::
1104 sage: set_random_seed() # long time
1105 sage: J = random_eja() # long time
1106 sage: caches = J.rank() # long time
1107 sage: J.rank.clear_cache() # long time
1108 sage: J.rank() == cached # long time
1112 return len(self
._charpoly
_coefficients
())
1115 def vector_space(self
):
1117 Return the vector space that underlies this algebra.
1121 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1125 sage: J = RealSymmetricEJA(2)
1126 sage: J.vector_space()
1127 Vector space of dimension 3 over...
1130 return self
.zero().to_vector().parent().ambient_vector_space()
1133 Element
= FiniteDimensionalEuclideanJordanAlgebraElement
1135 class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra
):
1137 New class for algebras whose supplied basis elements have all rational entries.
1141 sage: from mjo.eja.eja_algebra import BilinearFormEJA
1145 The supplied basis is orthonormalized by default::
1147 sage: B = matrix(QQ, [[1, 0, 0], [0, 25, -32], [0, -32, 41]])
1148 sage: J = BilinearFormEJA(B)
1149 sage: J.matrix_basis()
1162 orthonormalize
=True,
1169 # Abuse the check_field parameter to check that the entries of
1170 # out basis (in ambient coordinates) are in the field QQ.
1171 if not all( all(b_i
in QQ
for b_i
in b
.list()) for b
in basis
):
1172 raise TypeError("basis not rational")
1175 vector_basis
= basis
1177 from sage
.structure
.element
import is_Matrix
1178 basis_is_matrices
= False
1182 if is_Matrix(basis
[0]):
1183 basis_is_matrices
= True
1184 from mjo
.eja
.eja_utils
import _vec2mat
1185 vector_basis
= tuple( map(_mat2vec
,basis
) )
1186 degree
= basis
[0].nrows()**2
1188 degree
= basis
[0].degree()
1190 V
= VectorSpace(field
, degree
)
1192 # If we were asked to orthonormalize, and if the orthonormal
1193 # basis is different from the given one, then we also want to
1194 # compute multiplication and inner-product tables for the
1195 # deorthonormalized basis. These can be used later to
1196 # construct a deorthonormalized copy of this algebra over QQ
1197 # in which several operations are much faster.
1198 self
._rational
_algebra
= None
1201 if self
.base_ring() is not QQ
:
1202 # There's no point in constructing the extra algebra if this
1203 # one is already rational. If the original basis is rational
1204 # but normalization would make it irrational, then this whole
1205 # constructor will just fail anyway as it tries to stick an
1206 # irrational number into a rational algebra.
1208 # Note: the same Jordan and inner-products work here,
1209 # because they are necessarily defined with respect to
1210 # ambient coordinates and not any particular basis.
1211 self
._rational
_algebra
= RationalBasisEuclideanJordanAlgebra(
1216 orthonormalize
=False,
1222 # Compute the deorthonormalized tables before we orthonormalize
1224 W
= V
.span_of_basis( vector_basis
)
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 W
= V
.span_of_basis( vector_basis
)
1262 # Normalize the "matrix" basis, too!
1263 basis
= vector_basis
1265 if basis_is_matrices
:
1266 basis
= tuple( map(_vec2mat
,basis
) )
1268 W
= V
.span_of_basis( vector_basis
)
1270 # Now "W" is the vector space of our algebra coordinates. The
1271 # variables "X1", "X2",... refer to the entries of vectors in
1272 # W. Thus to convert back and forth between the orthonormal
1273 # coordinates and the given ones, we need to stick the original
1275 U
= V
.span_of_basis( deortho_vector_basis
)
1276 self
._deortho
_matrix
= matrix( U
.coordinate_vector(q
)
1277 for q
in vector_basis
)
1279 # If the superclass constructor is going to verify the
1280 # symmetry of this table, it has better at least be
1283 mult_table
= [ [0 for j
in range(n
)] for i
in range(n
) ]
1284 ip_table
= [ [0 for j
in range(n
)] for i
in range(n
) ]
1286 mult_table
= [ [0 for j
in range(i
+1)] for i
in range(n
) ]
1287 ip_table
= [ [0 for j
in range(i
+1)] for i
in range(n
) ]
1289 # Note: the Jordan and inner-products are defined in terms
1290 # of the ambient basis. It's important that their arguments
1291 # are in ambient coordinates as well.
1293 for j
in range(i
+1):
1294 # ortho basis w.r.t. ambient coords
1295 q_i
= vector_basis
[i
]
1296 q_j
= vector_basis
[j
]
1298 if basis_is_matrices
:
1302 elt
= jordan_product(q_i
, q_j
)
1303 ip
= inner_product(q_i
, q_j
)
1305 if basis_is_matrices
:
1306 # do another mat2vec because the multiplication
1307 # table is in terms of vectors
1310 elt
= W
.coordinate_vector(elt
)
1311 mult_table
[i
][j
] = elt
1314 # The tables are square if we're verifying that they
1316 mult_table
[j
][i
] = elt
1319 if basis_is_matrices
:
1323 basis
= tuple( x
.column() for x
in basis
)
1325 super().__init
__(field
,
1330 basis
, # matrix basis
1335 def _charpoly_coefficients(self
):
1339 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
1340 ....: JordanSpinEJA)
1344 The base ring of the resulting polynomial coefficients is what
1345 it should be, and not the rationals (unless the algebra was
1346 already over the rationals)::
1348 sage: J = JordanSpinEJA(3)
1349 sage: J._charpoly_coefficients()
1350 (X1^2 - X2^2 - X3^2, -2*X1)
1351 sage: a0 = J._charpoly_coefficients()[0]
1353 Algebraic Real Field
1354 sage: a0.base_ring()
1355 Algebraic Real Field
1358 if self
.base_ring() is QQ
or self
._rational
_algebra
is None:
1359 # There's no need to construct *another* algebra over the
1360 # rationals if this one is already over the
1361 # rationals. Likewise, if we never orthonormalized our
1362 # basis, we might as well just use the given one.
1363 superclass
= super(RationalBasisEuclideanJordanAlgebra
, self
)
1364 return superclass
._charpoly
_coefficients
()
1366 # Do the computation over the rationals. The answer will be
1367 # the same, because all we've done is a change of basis.
1368 # Then, change back from QQ to our real base ring
1369 a
= ( a_i
.change_ring(self
.base_ring())
1370 for a_i
in self
._rational
_algebra
._charpoly
_coefficients
() )
1372 # Now convert the coordinate variables back to the
1373 # deorthonormalized ones.
1374 R
= self
.coordinate_polynomial_ring()
1375 from sage
.modules
.free_module_element
import vector
1376 X
= vector(R
, R
.gens())
1377 BX
= self
._deortho
_matrix
*X
1379 subs_dict
= { X[i]: BX[i] for i in range(len(X)) }
1380 return tuple( a_i
.subs(subs_dict
) for a_i
in a
)
1382 class ConcreteEuclideanJordanAlgebra(RationalBasisEuclideanJordanAlgebra
):
1384 A class for the Euclidean Jordan algebras that we know by name.
1386 These are the Jordan algebras whose basis, multiplication table,
1387 rank, and so on are known a priori. More to the point, they are
1388 the Euclidean Jordan algebras for which we are able to conjure up
1389 a "random instance."
1393 sage: from mjo.eja.eja_algebra import ConcreteEuclideanJordanAlgebra
1397 Our basis is normalized with respect to the algebra's inner
1398 product, unless we specify otherwise::
1400 sage: set_random_seed()
1401 sage: J = ConcreteEuclideanJordanAlgebra.random_instance()
1402 sage: all( b.norm() == 1 for b in J.gens() )
1405 Since our basis is orthonormal with respect to the algebra's inner
1406 product, and since we know that this algebra is an EJA, any
1407 left-multiplication operator's matrix will be symmetric because
1408 natural->EJA basis representation is an isometry and within the
1409 EJA the operator is self-adjoint by the Jordan axiom::
1411 sage: set_random_seed()
1412 sage: J = ConcreteEuclideanJordanAlgebra.random_instance()
1413 sage: x = J.random_element()
1414 sage: x.operator().is_self_adjoint()
1419 def _max_random_instance_size():
1421 Return an integer "size" that is an upper bound on the size of
1422 this algebra when it is used in a random test
1423 case. Unfortunately, the term "size" is ambiguous -- when
1424 dealing with `R^n` under either the Hadamard or Jordan spin
1425 product, the "size" refers to the dimension `n`. When dealing
1426 with a matrix algebra (real symmetric or complex/quaternion
1427 Hermitian), it refers to the size of the matrix, which is far
1428 less than the dimension of the underlying vector space.
1430 This method must be implemented in each subclass.
1432 raise NotImplementedError
1435 def random_instance(cls
, *args
, **kwargs
):
1437 Return a random instance of this type of algebra.
1439 This method should be implemented in each subclass.
1441 from sage
.misc
.prandom
import choice
1442 eja_class
= choice(cls
.__subclasses
__())
1444 # These all bubble up to the RationalBasisEuclideanJordanAlgebra
1445 # superclass constructor, so any (kw)args valid there are also
1447 return eja_class
.random_instance(*args
, **kwargs
)
1450 class MatrixEuclideanJordanAlgebra
:
1454 Embed the matrix ``M`` into a space of real matrices.
1456 The matrix ``M`` can have entries in any field at the moment:
1457 the real numbers, complex numbers, or quaternions. And although
1458 they are not a field, we can probably support octonions at some
1459 point, too. This function returns a real matrix that "acts like"
1460 the original with respect to matrix multiplication; i.e.
1462 real_embed(M*N) = real_embed(M)*real_embed(N)
1465 raise NotImplementedError
1469 def real_unembed(M
):
1471 The inverse of :meth:`real_embed`.
1473 raise NotImplementedError
1476 def jordan_product(X
,Y
):
1477 return (X
*Y
+ Y
*X
)/2
1480 def trace_inner_product(cls
,X
,Y
):
1481 Xu
= cls
.real_unembed(X
)
1482 Yu
= cls
.real_unembed(Y
)
1483 tr
= (Xu
*Yu
).trace()
1486 # Works in QQ, AA, RDF, et cetera.
1488 except AttributeError:
1489 # A quaternion doesn't have a real() method, but does
1490 # have coefficient_tuple() method that returns the
1491 # coefficients of 1, i, j, and k -- in that order.
1492 return tr
.coefficient_tuple()[0]
1495 class RealMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1499 The identity function, for embedding real matrices into real
1505 def real_unembed(M
):
1507 The identity function, for unembedding real matrices from real
1513 class RealSymmetricEJA(ConcreteEuclideanJordanAlgebra
,
1514 RealMatrixEuclideanJordanAlgebra
):
1516 The rank-n simple EJA consisting of real symmetric n-by-n
1517 matrices, the usual symmetric Jordan product, and the trace inner
1518 product. It has dimension `(n^2 + n)/2` over the reals.
1522 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1526 sage: J = RealSymmetricEJA(2)
1527 sage: e0, e1, e2 = J.gens()
1535 In theory, our "field" can be any subfield of the reals::
1537 sage: RealSymmetricEJA(2, RDF)
1538 Euclidean Jordan algebra of dimension 3 over Real Double Field
1539 sage: RealSymmetricEJA(2, RR)
1540 Euclidean Jordan algebra of dimension 3 over Real Field with
1541 53 bits of precision
1545 The dimension of this algebra is `(n^2 + n) / 2`::
1547 sage: set_random_seed()
1548 sage: n_max = RealSymmetricEJA._max_random_instance_size()
1549 sage: n = ZZ.random_element(1, n_max)
1550 sage: J = RealSymmetricEJA(n)
1551 sage: J.dimension() == (n^2 + n)/2
1554 The Jordan multiplication is what we think it is::
1556 sage: set_random_seed()
1557 sage: J = RealSymmetricEJA.random_instance()
1558 sage: x,y = J.random_elements(2)
1559 sage: actual = (x*y).to_matrix()
1560 sage: X = x.to_matrix()
1561 sage: Y = y.to_matrix()
1562 sage: expected = (X*Y + Y*X)/2
1563 sage: actual == expected
1565 sage: J(expected) == x*y
1568 We can change the generator prefix::
1570 sage: RealSymmetricEJA(3, prefix='q').gens()
1571 (q0, q1, q2, q3, q4, q5)
1573 We can construct the (trivial) algebra of rank zero::
1575 sage: RealSymmetricEJA(0)
1576 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1580 def _denormalized_basis(cls
, n
, field
):
1582 Return a basis for the space of real symmetric n-by-n matrices.
1586 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1590 sage: set_random_seed()
1591 sage: n = ZZ.random_element(1,5)
1592 sage: B = RealSymmetricEJA._denormalized_basis(n,QQ)
1593 sage: all( M.is_symmetric() for M in B)
1597 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1601 for j
in range(i
+1):
1602 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
1606 Sij
= Eij
+ Eij
.transpose()
1612 def _max_random_instance_size():
1613 return 4 # Dimension 10
1616 def random_instance(cls
, field
=AA
, **kwargs
):
1618 Return a random instance of this type of algebra.
1620 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1621 return cls(n
, field
, **kwargs
)
1623 def __init__(self
, n
, field
=AA
, **kwargs
):
1624 basis
= self
._denormalized
_basis
(n
, field
)
1625 super(RealSymmetricEJA
, self
).__init
__(field
,
1627 self
.jordan_product
,
1628 self
.trace_inner_product
,
1630 self
.rank
.set_cache(n
)
1631 self
.one
.set_cache(self(matrix
.identity(field
,n
)))
1634 class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1638 Embed the n-by-n complex matrix ``M`` into the space of real
1639 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1640 bi` to the block matrix ``[[a,b],[-b,a]]``.
1644 sage: from mjo.eja.eja_algebra import \
1645 ....: ComplexMatrixEuclideanJordanAlgebra
1649 sage: F = QuadraticField(-1, 'I')
1650 sage: x1 = F(4 - 2*i)
1651 sage: x2 = F(1 + 2*i)
1654 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1655 sage: ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1664 Embedding is a homomorphism (isomorphism, in fact)::
1666 sage: set_random_seed()
1667 sage: n = ZZ.random_element(3)
1668 sage: F = QuadraticField(-1, 'I')
1669 sage: X = random_matrix(F, n)
1670 sage: Y = random_matrix(F, n)
1671 sage: Xe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X)
1672 sage: Ye = ComplexMatrixEuclideanJordanAlgebra.real_embed(Y)
1673 sage: XYe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1680 raise ValueError("the matrix 'M' must be square")
1682 # We don't need any adjoined elements...
1683 field
= M
.base_ring().base_ring()
1687 a
= z
.list()[0] # real part, I guess
1688 b
= z
.list()[1] # imag part, I guess
1689 blocks
.append(matrix(field
, 2, [[a
,b
],[-b
,a
]]))
1691 return matrix
.block(field
, n
, blocks
)
1695 def real_unembed(M
):
1697 The inverse of _embed_complex_matrix().
1701 sage: from mjo.eja.eja_algebra import \
1702 ....: ComplexMatrixEuclideanJordanAlgebra
1706 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1707 ....: [-2, 1, -4, 3],
1708 ....: [ 9, 10, 11, 12],
1709 ....: [-10, 9, -12, 11] ])
1710 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(A)
1712 [ 10*I + 9 12*I + 11]
1716 Unembedding is the inverse of embedding::
1718 sage: set_random_seed()
1719 sage: F = QuadraticField(-1, 'I')
1720 sage: M = random_matrix(F, 3)
1721 sage: Me = ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1722 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1728 raise ValueError("the matrix 'M' must be square")
1729 if not n
.mod(2).is_zero():
1730 raise ValueError("the matrix 'M' must be a complex embedding")
1732 # If "M" was normalized, its base ring might have roots
1733 # adjoined and they can stick around after unembedding.
1734 field
= M
.base_ring()
1735 R
= PolynomialRing(field
, 'z')
1738 # Sage doesn't know how to embed AA into QQbar, i.e. how
1739 # to adjoin sqrt(-1) to AA.
1742 F
= field
.extension(z
**2 + 1, 'I', embedding
=CLF(-1).sqrt())
1745 # Go top-left to bottom-right (reading order), converting every
1746 # 2-by-2 block we see to a single complex element.
1748 for k
in range(n
/2):
1749 for j
in range(n
/2):
1750 submat
= M
[2*k
:2*k
+2,2*j
:2*j
+2]
1751 if submat
[0,0] != submat
[1,1]:
1752 raise ValueError('bad on-diagonal submatrix')
1753 if submat
[0,1] != -submat
[1,0]:
1754 raise ValueError('bad off-diagonal submatrix')
1755 z
= submat
[0,0] + submat
[0,1]*i
1758 return matrix(F
, n
/2, elements
)
1762 def trace_inner_product(cls
,X
,Y
):
1764 Compute a matrix inner product in this algebra directly from
1769 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1773 This gives the same answer as the slow, default method implemented
1774 in :class:`MatrixEuclideanJordanAlgebra`::
1776 sage: set_random_seed()
1777 sage: J = ComplexHermitianEJA.random_instance()
1778 sage: x,y = J.random_elements(2)
1779 sage: Xe = x.to_matrix()
1780 sage: Ye = y.to_matrix()
1781 sage: X = ComplexHermitianEJA.real_unembed(Xe)
1782 sage: Y = ComplexHermitianEJA.real_unembed(Ye)
1783 sage: expected = (X*Y).trace().real()
1784 sage: actual = ComplexHermitianEJA.trace_inner_product(Xe,Ye)
1785 sage: actual == expected
1789 return RealMatrixEuclideanJordanAlgebra
.trace_inner_product(X
,Y
)/2
1792 class ComplexHermitianEJA(ConcreteEuclideanJordanAlgebra
,
1793 ComplexMatrixEuclideanJordanAlgebra
):
1795 The rank-n simple EJA consisting of complex Hermitian n-by-n
1796 matrices over the real numbers, the usual symmetric Jordan product,
1797 and the real-part-of-trace inner product. It has dimension `n^2` over
1802 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1806 In theory, our "field" can be any subfield of the reals::
1808 sage: ComplexHermitianEJA(2, RDF)
1809 Euclidean Jordan algebra of dimension 4 over Real Double Field
1810 sage: ComplexHermitianEJA(2, RR)
1811 Euclidean Jordan algebra of dimension 4 over Real Field with
1812 53 bits of precision
1816 The dimension of this algebra is `n^2`::
1818 sage: set_random_seed()
1819 sage: n_max = ComplexHermitianEJA._max_random_instance_size()
1820 sage: n = ZZ.random_element(1, n_max)
1821 sage: J = ComplexHermitianEJA(n)
1822 sage: J.dimension() == n^2
1825 The Jordan multiplication is what we think it is::
1827 sage: set_random_seed()
1828 sage: J = ComplexHermitianEJA.random_instance()
1829 sage: x,y = J.random_elements(2)
1830 sage: actual = (x*y).to_matrix()
1831 sage: X = x.to_matrix()
1832 sage: Y = y.to_matrix()
1833 sage: expected = (X*Y + Y*X)/2
1834 sage: actual == expected
1836 sage: J(expected) == x*y
1839 We can change the generator prefix::
1841 sage: ComplexHermitianEJA(2, prefix='z').gens()
1844 We can construct the (trivial) algebra of rank zero::
1846 sage: ComplexHermitianEJA(0)
1847 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1852 def _denormalized_basis(cls
, n
, field
):
1854 Returns a basis for the space of complex Hermitian n-by-n matrices.
1856 Why do we embed these? Basically, because all of numerical linear
1857 algebra assumes that you're working with vectors consisting of `n`
1858 entries from a field and scalars from the same field. There's no way
1859 to tell SageMath that (for example) the vectors contain complex
1860 numbers, while the scalar field is real.
1864 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1868 sage: set_random_seed()
1869 sage: n = ZZ.random_element(1,5)
1870 sage: field = QuadraticField(2, 'sqrt2')
1871 sage: B = ComplexHermitianEJA._denormalized_basis(n, field)
1872 sage: all( M.is_symmetric() for M in B)
1876 R
= PolynomialRing(field
, 'z')
1878 F
= field
.extension(z
**2 + 1, 'I')
1881 # This is like the symmetric case, but we need to be careful:
1883 # * We want conjugate-symmetry, not just symmetry.
1884 # * The diagonal will (as a result) be real.
1888 for j
in range(i
+1):
1889 Eij
= matrix(F
, n
, lambda k
,l
: k
==i
and l
==j
)
1891 Sij
= cls
.real_embed(Eij
)
1894 # The second one has a minus because it's conjugated.
1895 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
1897 Sij_imag
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
1900 # Since we embedded these, we can drop back to the "field" that we
1901 # started with instead of the complex extension "F".
1902 return tuple( s
.change_ring(field
) for s
in S
)
1905 def __init__(self
, n
, field
=AA
, **kwargs
):
1906 basis
= self
._denormalized
_basis
(n
,field
)
1907 super(ComplexHermitianEJA
, self
).__init
__(field
,
1909 self
.jordan_product
,
1910 self
.trace_inner_product
,
1912 self
.rank
.set_cache(n
)
1913 # TODO: pre-cache the identity!
1916 def _max_random_instance_size():
1917 return 3 # Dimension 9
1920 def random_instance(cls
, field
=AA
, **kwargs
):
1922 Return a random instance of this type of algebra.
1924 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1925 return cls(n
, field
, **kwargs
)
1927 class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1931 Embed the n-by-n quaternion matrix ``M`` into the space of real
1932 matrices of size 4n-by-4n by first sending each quaternion entry `z
1933 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
1934 c+di],[-c + di, a-bi]]`, and then embedding those into a real
1939 sage: from mjo.eja.eja_algebra import \
1940 ....: QuaternionMatrixEuclideanJordanAlgebra
1944 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1945 sage: i,j,k = Q.gens()
1946 sage: x = 1 + 2*i + 3*j + 4*k
1947 sage: M = matrix(Q, 1, [[x]])
1948 sage: QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1954 Embedding is a homomorphism (isomorphism, in fact)::
1956 sage: set_random_seed()
1957 sage: n = ZZ.random_element(2)
1958 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1959 sage: X = random_matrix(Q, n)
1960 sage: Y = random_matrix(Q, n)
1961 sage: Xe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X)
1962 sage: Ye = QuaternionMatrixEuclideanJordanAlgebra.real_embed(Y)
1963 sage: XYe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1968 quaternions
= M
.base_ring()
1971 raise ValueError("the matrix 'M' must be square")
1973 F
= QuadraticField(-1, 'I')
1978 t
= z
.coefficient_tuple()
1983 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
1984 [-c
+ d
*i
, a
- b
*i
]])
1985 realM
= ComplexMatrixEuclideanJordanAlgebra
.real_embed(cplxM
)
1986 blocks
.append(realM
)
1988 # We should have real entries by now, so use the realest field
1989 # we've got for the return value.
1990 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
1995 def real_unembed(M
):
1997 The inverse of _embed_quaternion_matrix().
2001 sage: from mjo.eja.eja_algebra import \
2002 ....: QuaternionMatrixEuclideanJordanAlgebra
2006 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
2007 ....: [-2, 1, -4, 3],
2008 ....: [-3, 4, 1, -2],
2009 ....: [-4, -3, 2, 1]])
2010 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(M)
2011 [1 + 2*i + 3*j + 4*k]
2015 Unembedding is the inverse of embedding::
2017 sage: set_random_seed()
2018 sage: Q = QuaternionAlgebra(QQ, -1, -1)
2019 sage: M = random_matrix(Q, 3)
2020 sage: Me = QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
2021 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
2027 raise ValueError("the matrix 'M' must be square")
2028 if not n
.mod(4).is_zero():
2029 raise ValueError("the matrix 'M' must be a quaternion embedding")
2031 # Use the base ring of the matrix to ensure that its entries can be
2032 # multiplied by elements of the quaternion algebra.
2033 field
= M
.base_ring()
2034 Q
= QuaternionAlgebra(field
,-1,-1)
2037 # Go top-left to bottom-right (reading order), converting every
2038 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
2041 for l
in range(n
/4):
2042 for m
in range(n
/4):
2043 submat
= ComplexMatrixEuclideanJordanAlgebra
.real_unembed(
2044 M
[4*l
:4*l
+4,4*m
:4*m
+4] )
2045 if submat
[0,0] != submat
[1,1].conjugate():
2046 raise ValueError('bad on-diagonal submatrix')
2047 if submat
[0,1] != -submat
[1,0].conjugate():
2048 raise ValueError('bad off-diagonal submatrix')
2049 z
= submat
[0,0].real()
2050 z
+= submat
[0,0].imag()*i
2051 z
+= submat
[0,1].real()*j
2052 z
+= submat
[0,1].imag()*k
2055 return matrix(Q
, n
/4, elements
)
2059 def trace_inner_product(cls
,X
,Y
):
2061 Compute a matrix inner product in this algebra directly from
2066 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2070 This gives the same answer as the slow, default method implemented
2071 in :class:`MatrixEuclideanJordanAlgebra`::
2073 sage: set_random_seed()
2074 sage: J = QuaternionHermitianEJA.random_instance()
2075 sage: x,y = J.random_elements(2)
2076 sage: Xe = x.to_matrix()
2077 sage: Ye = y.to_matrix()
2078 sage: X = QuaternionHermitianEJA.real_unembed(Xe)
2079 sage: Y = QuaternionHermitianEJA.real_unembed(Ye)
2080 sage: expected = (X*Y).trace().coefficient_tuple()[0]
2081 sage: actual = QuaternionHermitianEJA.trace_inner_product(Xe,Ye)
2082 sage: actual == expected
2086 return RealMatrixEuclideanJordanAlgebra
.trace_inner_product(X
,Y
)/4
2089 class QuaternionHermitianEJA(ConcreteEuclideanJordanAlgebra
,
2090 QuaternionMatrixEuclideanJordanAlgebra
):
2092 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
2093 matrices, the usual symmetric Jordan product, and the
2094 real-part-of-trace inner product. It has dimension `2n^2 - n` over
2099 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2103 In theory, our "field" can be any subfield of the reals::
2105 sage: QuaternionHermitianEJA(2, RDF)
2106 Euclidean Jordan algebra of dimension 6 over Real Double Field
2107 sage: QuaternionHermitianEJA(2, RR)
2108 Euclidean Jordan algebra of dimension 6 over Real Field with
2109 53 bits of precision
2113 The dimension of this algebra is `2*n^2 - n`::
2115 sage: set_random_seed()
2116 sage: n_max = QuaternionHermitianEJA._max_random_instance_size()
2117 sage: n = ZZ.random_element(1, n_max)
2118 sage: J = QuaternionHermitianEJA(n)
2119 sage: J.dimension() == 2*(n^2) - n
2122 The Jordan multiplication is what we think it is::
2124 sage: set_random_seed()
2125 sage: J = QuaternionHermitianEJA.random_instance()
2126 sage: x,y = J.random_elements(2)
2127 sage: actual = (x*y).to_matrix()
2128 sage: X = x.to_matrix()
2129 sage: Y = y.to_matrix()
2130 sage: expected = (X*Y + Y*X)/2
2131 sage: actual == expected
2133 sage: J(expected) == x*y
2136 We can change the generator prefix::
2138 sage: QuaternionHermitianEJA(2, prefix='a').gens()
2139 (a0, a1, a2, a3, a4, a5)
2141 We can construct the (trivial) algebra of rank zero::
2143 sage: QuaternionHermitianEJA(0)
2144 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2148 def _denormalized_basis(cls
, n
, field
):
2150 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
2152 Why do we embed these? Basically, because all of numerical
2153 linear algebra assumes that you're working with vectors consisting
2154 of `n` entries from a field and scalars from the same field. There's
2155 no way to tell SageMath that (for example) the vectors contain
2156 complex numbers, while the scalar field is real.
2160 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2164 sage: set_random_seed()
2165 sage: n = ZZ.random_element(1,5)
2166 sage: B = QuaternionHermitianEJA._denormalized_basis(n,QQ)
2167 sage: all( M.is_symmetric() for M in B )
2171 Q
= QuaternionAlgebra(QQ
,-1,-1)
2174 # This is like the symmetric case, but we need to be careful:
2176 # * We want conjugate-symmetry, not just symmetry.
2177 # * The diagonal will (as a result) be real.
2181 for j
in range(i
+1):
2182 Eij
= matrix(Q
, n
, lambda k
,l
: k
==i
and l
==j
)
2184 Sij
= cls
.real_embed(Eij
)
2187 # The second, third, and fourth ones have a minus
2188 # because they're conjugated.
2189 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
2191 Sij_I
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
2193 Sij_J
= cls
.real_embed(J
*Eij
- J
*Eij
.transpose())
2195 Sij_K
= cls
.real_embed(K
*Eij
- K
*Eij
.transpose())
2198 # Since we embedded these, we can drop back to the "field" that we
2199 # started with instead of the quaternion algebra "Q".
2200 return tuple( s
.change_ring(field
) for s
in S
)
2203 def __init__(self
, n
, field
=AA
, **kwargs
):
2204 basis
= self
._denormalized
_basis
(n
,field
)
2205 super(QuaternionHermitianEJA
, self
).__init
__(field
,
2207 self
.jordan_product
,
2208 self
.trace_inner_product
,
2210 self
.rank
.set_cache(n
)
2211 # TODO: cache one()!
2214 def _max_random_instance_size():
2216 The maximum rank of a random QuaternionHermitianEJA.
2218 return 2 # Dimension 6
2221 def random_instance(cls
, field
=AA
, **kwargs
):
2223 Return a random instance of this type of algebra.
2225 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2226 return cls(n
, field
, **kwargs
)
2229 class HadamardEJA(ConcreteEuclideanJordanAlgebra
):
2231 Return the Euclidean Jordan Algebra corresponding to the set
2232 `R^n` under the Hadamard product.
2234 Note: this is nothing more than the Cartesian product of ``n``
2235 copies of the spin algebra. Once Cartesian product algebras
2236 are implemented, this can go.
2240 sage: from mjo.eja.eja_algebra import HadamardEJA
2244 This multiplication table can be verified by hand::
2246 sage: J = HadamardEJA(3)
2247 sage: e0,e1,e2 = J.gens()
2263 We can change the generator prefix::
2265 sage: HadamardEJA(3, prefix='r').gens()
2269 def __init__(self
, n
, field
=AA
, **kwargs
):
2270 V
= VectorSpace(field
, n
)
2273 def jordan_product(x
,y
):
2274 return V([ xi
*yi
for (xi
,yi
) in zip(x
,y
) ])
2275 def inner_product(x
,y
):
2276 return x
.inner_product(y
)
2278 # Don't orthonormalize because our basis is already orthonormal
2279 # with respect to our inner-product.
2280 super(HadamardEJA
, self
).__init
__(field
,
2284 orthonormalize
=False,
2288 self
.rank
.set_cache(n
)
2291 self
.one
.set_cache( self
.zero() )
2293 self
.one
.set_cache( sum(self
.gens()) )
2296 def _max_random_instance_size():
2298 The maximum dimension of a random HadamardEJA.
2303 def random_instance(cls
, field
=AA
, **kwargs
):
2305 Return a random instance of this type of algebra.
2307 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2308 return cls(n
, field
, **kwargs
)
2311 class BilinearFormEJA(ConcreteEuclideanJordanAlgebra
):
2313 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2314 with the half-trace inner product and jordan product ``x*y =
2315 (<Bx,y>,y_bar>, x0*y_bar + y0*x_bar)`` where `B = 1 \times B22` is
2316 a symmetric positive-definite "bilinear form" matrix. Its
2317 dimension is the size of `B`, and it has rank two in dimensions
2318 larger than two. It reduces to the ``JordanSpinEJA`` when `B` is
2319 the identity matrix of order ``n``.
2321 We insist that the one-by-one upper-left identity block of `B` be
2322 passed in as well so that we can be passed a matrix of size zero
2323 to construct a trivial algebra.
2327 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
2328 ....: JordanSpinEJA)
2332 When no bilinear form is specified, the identity matrix is used,
2333 and the resulting algebra is the Jordan spin algebra::
2335 sage: B = matrix.identity(AA,3)
2336 sage: J0 = BilinearFormEJA(B)
2337 sage: J1 = JordanSpinEJA(3)
2338 sage: J0.multiplication_table() == J0.multiplication_table()
2341 An error is raised if the matrix `B` does not correspond to a
2342 positive-definite bilinear form::
2344 sage: B = matrix.random(QQ,2,3)
2345 sage: J = BilinearFormEJA(B)
2346 Traceback (most recent call last):
2348 ValueError: bilinear form is not positive-definite
2349 sage: B = matrix.zero(QQ,3)
2350 sage: J = BilinearFormEJA(B)
2351 Traceback (most recent call last):
2353 ValueError: bilinear form is not positive-definite
2357 We can create a zero-dimensional algebra::
2359 sage: B = matrix.identity(AA,0)
2360 sage: J = BilinearFormEJA(B)
2364 We can check the multiplication condition given in the Jordan, von
2365 Neumann, and Wigner paper (and also discussed on my "On the
2366 symmetry..." paper). Note that this relies heavily on the standard
2367 choice of basis, as does anything utilizing the bilinear form
2368 matrix. We opt not to orthonormalize the basis, because if we
2369 did, we would have to normalize the `s_{i}` in a similar manner::
2371 sage: set_random_seed()
2372 sage: n = ZZ.random_element(5)
2373 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2374 sage: B11 = matrix.identity(QQ,1)
2375 sage: B22 = M.transpose()*M
2376 sage: B = block_matrix(2,2,[ [B11,0 ],
2378 sage: J = BilinearFormEJA(B, orthonormalize=False)
2379 sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
2380 sage: V = J.vector_space()
2381 sage: sis = [ J( V([0] + (M.inverse()*ei).list()).column() )
2382 ....: for ei in eis ]
2383 sage: actual = [ sis[i]*sis[j]
2384 ....: for i in range(n-1)
2385 ....: for j in range(n-1) ]
2386 sage: expected = [ J.one() if i == j else J.zero()
2387 ....: for i in range(n-1)
2388 ....: for j in range(n-1) ]
2389 sage: actual == expected
2392 def __init__(self
, B
, field
=AA
, **kwargs
):
2393 if not B
.is_positive_definite():
2394 raise ValueError("bilinear form is not positive-definite")
2397 V
= VectorSpace(field
, n
)
2399 def inner_product(x
,y
):
2400 return (B
*x
).inner_product(y
)
2402 def jordan_product(x
,y
):
2407 z0
= inner_product(x
,y
)
2408 zbar
= y0
*xbar
+ x0
*ybar
2409 return V([z0
] + zbar
.list())
2411 super(BilinearFormEJA
, self
).__init
__(field
,
2417 # The rank of this algebra is two, unless we're in a
2418 # one-dimensional ambient space (because the rank is bounded
2419 # by the ambient dimension).
2420 self
.rank
.set_cache(min(n
,2))
2423 self
.one
.set_cache( self
.zero() )
2425 self
.one
.set_cache( self
.monomial(0) )
2428 def _max_random_instance_size():
2430 The maximum dimension of a random BilinearFormEJA.
2435 def random_instance(cls
, field
=AA
, **kwargs
):
2437 Return a random instance of this algebra.
2439 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2441 B
= matrix
.identity(field
, n
)
2442 return cls(B
, field
, **kwargs
)
2444 B11
= matrix
.identity(field
,1)
2445 M
= matrix
.random(field
, n
-1)
2446 I
= matrix
.identity(field
, n
-1)
2447 alpha
= field
.zero()
2448 while alpha
.is_zero():
2449 alpha
= field
.random_element().abs()
2450 B22
= M
.transpose()*M
+ alpha
*I
2452 from sage
.matrix
.special
import block_matrix
2453 B
= block_matrix(2,2, [ [B11
, ZZ(0) ],
2456 return cls(B
, field
, **kwargs
)
2459 class JordanSpinEJA(BilinearFormEJA
):
2461 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2462 with the usual inner product and jordan product ``x*y =
2463 (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2468 sage: from mjo.eja.eja_algebra import JordanSpinEJA
2472 This multiplication table can be verified by hand::
2474 sage: J = JordanSpinEJA(4)
2475 sage: e0,e1,e2,e3 = J.gens()
2491 We can change the generator prefix::
2493 sage: JordanSpinEJA(2, prefix='B').gens()
2498 Ensure that we have the usual inner product on `R^n`::
2500 sage: set_random_seed()
2501 sage: J = JordanSpinEJA.random_instance()
2502 sage: x,y = J.random_elements(2)
2503 sage: actual = x.inner_product(y)
2504 sage: expected = x.to_vector().inner_product(y.to_vector())
2505 sage: actual == expected
2509 def __init__(self
, n
, field
=AA
, **kwargs
):
2510 # This is a special case of the BilinearFormEJA with the
2511 # identity matrix as its bilinear form.
2512 B
= matrix
.identity(field
, n
)
2514 # Don't orthonormalize because our basis is already
2515 # orthonormal with respect to our inner-product.
2516 super(JordanSpinEJA
, self
).__init
__(B
,
2518 orthonormalize
=False,
2524 def _max_random_instance_size():
2526 The maximum dimension of a random JordanSpinEJA.
2531 def random_instance(cls
, field
=AA
, **kwargs
):
2533 Return a random instance of this type of algebra.
2535 Needed here to override the implementation for ``BilinearFormEJA``.
2537 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2538 return cls(n
, field
, **kwargs
)
2541 class TrivialEJA(ConcreteEuclideanJordanAlgebra
):
2543 The trivial Euclidean Jordan algebra consisting of only a zero element.
2547 sage: from mjo.eja.eja_algebra import TrivialEJA
2551 sage: J = TrivialEJA()
2558 sage: 7*J.one()*12*J.one()
2560 sage: J.one().inner_product(J.one())
2562 sage: J.one().norm()
2564 sage: J.one().subalgebra_generated_by()
2565 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2570 def __init__(self
, field
=AA
, **kwargs
):
2571 jordan_product
= lambda x
,y
: x
2572 inner_product
= lambda x
,y
: field(0)
2574 super(TrivialEJA
, self
).__init
__(field
,
2579 # The rank is zero using my definition, namely the dimension of the
2580 # largest subalgebra generated by any element.
2581 self
.rank
.set_cache(0)
2582 self
.one
.set_cache( self
.zero() )
2585 def random_instance(cls
, field
=AA
, **kwargs
):
2586 # We don't take a "size" argument so the superclass method is
2587 # inappropriate for us.
2588 return cls(field
, **kwargs
)
2590 class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra
):
2592 The external (orthogonal) direct sum of two other Euclidean Jordan
2593 algebras. Essentially the Cartesian product of its two factors.
2594 Every Euclidean Jordan algebra decomposes into an orthogonal
2595 direct sum of simple Euclidean Jordan algebras, so no generality
2596 is lost by providing only this construction.
2600 sage: from mjo.eja.eja_algebra import (random_eja,
2602 ....: RealSymmetricEJA,
2607 sage: J1 = HadamardEJA(2)
2608 sage: J2 = RealSymmetricEJA(3)
2609 sage: J = DirectSumEJA(J1,J2)
2617 The external direct sum construction is only valid when the two factors
2618 have the same base ring; an error is raised otherwise::
2620 sage: set_random_seed()
2621 sage: J1 = random_eja(AA)
2622 sage: J2 = random_eja(QQ,orthonormalize=False)
2623 sage: J = DirectSumEJA(J1,J2)
2624 Traceback (most recent call last):
2626 ValueError: algebras must share the same base field
2629 def __init__(self
, J1
, J2
, **kwargs
):
2630 if J1
.base_ring() != J2
.base_ring():
2631 raise ValueError("algebras must share the same base field")
2632 field
= J1
.base_ring()
2634 self
._factors
= (J1
, J2
)
2638 V
= VectorSpace(field
, n
)
2639 mult_table
= [ [ V
.zero() for j
in range(i
+1) ]
2642 for j
in range(i
+1):
2643 p
= (J1
.monomial(i
)*J1
.monomial(j
)).to_vector()
2644 mult_table
[i
][j
] = V(p
.list() + [field
.zero()]*n2
)
2647 for j
in range(i
+1):
2648 p
= (J2
.monomial(i
)*J2
.monomial(j
)).to_vector()
2649 mult_table
[n1
+i
][n1
+j
] = V([field
.zero()]*n1
+ p
.list())
2651 # TODO: build the IP table here from the two constituent IP
2652 # matrices (it'll be block diagonal, I think).
2653 ip_table
= [ [ field
.zero() for j
in range(i
+1) ]
2655 super(DirectSumEJA
, self
).__init
__(field
,
2660 self
.rank
.set_cache(J1
.rank() + J2
.rank())
2665 Return the pair of this algebra's factors.
2669 sage: from mjo.eja.eja_algebra import (HadamardEJA,
2670 ....: JordanSpinEJA,
2675 sage: J1 = HadamardEJA(2,QQ)
2676 sage: J2 = JordanSpinEJA(3,QQ)
2677 sage: J = DirectSumEJA(J1,J2)
2679 (Euclidean Jordan algebra of dimension 2 over Rational Field,
2680 Euclidean Jordan algebra of dimension 3 over Rational Field)
2683 return self
._factors
2685 def projections(self
):
2687 Return a pair of projections onto this algebra's factors.
2691 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
2692 ....: ComplexHermitianEJA,
2697 sage: J1 = JordanSpinEJA(2)
2698 sage: J2 = ComplexHermitianEJA(2)
2699 sage: J = DirectSumEJA(J1,J2)
2700 sage: (pi_left, pi_right) = J.projections()
2701 sage: J.one().to_vector()
2703 sage: pi_left(J.one()).to_vector()
2705 sage: pi_right(J.one()).to_vector()
2709 (J1
,J2
) = self
.factors()
2712 V_basis
= self
.vector_space().basis()
2713 # Need to specify the dimensions explicitly so that we don't
2714 # wind up with a zero-by-zero matrix when we want e.g. a
2715 # zero-by-two matrix (important for composing things).
2716 P1
= matrix(self
.base_ring(), m
, m
+n
, V_basis
[:m
])
2717 P2
= matrix(self
.base_ring(), n
, m
+n
, V_basis
[m
:])
2718 pi_left
= FiniteDimensionalEuclideanJordanAlgebraOperator(self
,J1
,P1
)
2719 pi_right
= FiniteDimensionalEuclideanJordanAlgebraOperator(self
,J2
,P2
)
2720 return (pi_left
, pi_right
)
2722 def inclusions(self
):
2724 Return the pair of inclusion maps from our factors into us.
2728 sage: from mjo.eja.eja_algebra import (random_eja,
2729 ....: JordanSpinEJA,
2730 ....: RealSymmetricEJA,
2735 sage: J1 = JordanSpinEJA(3)
2736 sage: J2 = RealSymmetricEJA(2)
2737 sage: J = DirectSumEJA(J1,J2)
2738 sage: (iota_left, iota_right) = J.inclusions()
2739 sage: iota_left(J1.zero()) == J.zero()
2741 sage: iota_right(J2.zero()) == J.zero()
2743 sage: J1.one().to_vector()
2745 sage: iota_left(J1.one()).to_vector()
2747 sage: J2.one().to_vector()
2749 sage: iota_right(J2.one()).to_vector()
2751 sage: J.one().to_vector()
2756 Composing a projection with the corresponding inclusion should
2757 produce the identity map, and mismatching them should produce
2760 sage: set_random_seed()
2761 sage: J1 = random_eja()
2762 sage: J2 = random_eja()
2763 sage: J = DirectSumEJA(J1,J2)
2764 sage: (iota_left, iota_right) = J.inclusions()
2765 sage: (pi_left, pi_right) = J.projections()
2766 sage: pi_left*iota_left == J1.one().operator()
2768 sage: pi_right*iota_right == J2.one().operator()
2770 sage: (pi_left*iota_right).is_zero()
2772 sage: (pi_right*iota_left).is_zero()
2776 (J1
,J2
) = self
.factors()
2779 V_basis
= self
.vector_space().basis()
2780 # Need to specify the dimensions explicitly so that we don't
2781 # wind up with a zero-by-zero matrix when we want e.g. a
2782 # two-by-zero matrix (important for composing things).
2783 I1
= matrix
.column(self
.base_ring(), m
, m
+n
, V_basis
[:m
])
2784 I2
= matrix
.column(self
.base_ring(), n
, m
+n
, V_basis
[m
:])
2785 iota_left
= FiniteDimensionalEuclideanJordanAlgebraOperator(J1
,self
,I1
)
2786 iota_right
= FiniteDimensionalEuclideanJordanAlgebraOperator(J2
,self
,I2
)
2787 return (iota_left
, iota_right
)
2789 def inner_product(self
, x
, y
):
2791 The standard Cartesian inner-product.
2793 We project ``x`` and ``y`` onto our factors, and add up the
2794 inner-products from the subalgebras.
2799 sage: from mjo.eja.eja_algebra import (HadamardEJA,
2800 ....: QuaternionHermitianEJA,
2805 sage: J1 = HadamardEJA(3,QQ)
2806 sage: J2 = QuaternionHermitianEJA(2,QQ,orthonormalize=False)
2807 sage: J = DirectSumEJA(J1,J2)
2812 sage: x1.inner_product(x2)
2814 sage: y1.inner_product(y2)
2816 sage: J.one().inner_product(J.one())
2820 (pi_left
, pi_right
) = self
.projections()
2826 return (x1
.inner_product(y1
) + x2
.inner_product(y2
))
2830 random_eja
= ConcreteEuclideanJordanAlgebra
.random_instance