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
:
1359 # There's no need to construct *another* algebra over the
1360 # rationals if this one is already over the rationals.
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
:
1452 Embed the matrix ``M`` into a space of real matrices.
1454 The matrix ``M`` can have entries in any field at the moment:
1455 the real numbers, complex numbers, or quaternions. And although
1456 they are not a field, we can probably support octonions at some
1457 point, too. This function returns a real matrix that "acts like"
1458 the original with respect to matrix multiplication; i.e.
1460 real_embed(M*N) = real_embed(M)*real_embed(N)
1463 raise NotImplementedError
1467 def real_unembed(M
):
1469 The inverse of :meth:`real_embed`.
1471 raise NotImplementedError
1474 def jordan_product(X
,Y
):
1475 return (X
*Y
+ Y
*X
)/2
1478 def trace_inner_product(cls
,X
,Y
):
1479 Xu
= cls
.real_unembed(X
)
1480 Yu
= cls
.real_unembed(Y
)
1481 tr
= (Xu
*Yu
).trace()
1484 # Works in QQ, AA, RDF, et cetera.
1486 except AttributeError:
1487 # A quaternion doesn't have a real() method, but does
1488 # have coefficient_tuple() method that returns the
1489 # coefficients of 1, i, j, and k -- in that order.
1490 return tr
.coefficient_tuple()[0]
1493 class RealMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1497 The identity function, for embedding real matrices into real
1503 def real_unembed(M
):
1505 The identity function, for unembedding real matrices from real
1511 class RealSymmetricEJA(ConcreteEuclideanJordanAlgebra
,
1512 RealMatrixEuclideanJordanAlgebra
):
1514 The rank-n simple EJA consisting of real symmetric n-by-n
1515 matrices, the usual symmetric Jordan product, and the trace inner
1516 product. It has dimension `(n^2 + n)/2` over the reals.
1520 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1524 sage: J = RealSymmetricEJA(2)
1525 sage: e0, e1, e2 = J.gens()
1533 In theory, our "field" can be any subfield of the reals::
1535 sage: RealSymmetricEJA(2, RDF)
1536 Euclidean Jordan algebra of dimension 3 over Real Double Field
1537 sage: RealSymmetricEJA(2, RR)
1538 Euclidean Jordan algebra of dimension 3 over Real Field with
1539 53 bits of precision
1543 The dimension of this algebra is `(n^2 + n) / 2`::
1545 sage: set_random_seed()
1546 sage: n_max = RealSymmetricEJA._max_random_instance_size()
1547 sage: n = ZZ.random_element(1, n_max)
1548 sage: J = RealSymmetricEJA(n)
1549 sage: J.dimension() == (n^2 + n)/2
1552 The Jordan multiplication is what we think it is::
1554 sage: set_random_seed()
1555 sage: J = RealSymmetricEJA.random_instance()
1556 sage: x,y = J.random_elements(2)
1557 sage: actual = (x*y).to_matrix()
1558 sage: X = x.to_matrix()
1559 sage: Y = y.to_matrix()
1560 sage: expected = (X*Y + Y*X)/2
1561 sage: actual == expected
1563 sage: J(expected) == x*y
1566 We can change the generator prefix::
1568 sage: RealSymmetricEJA(3, prefix='q').gens()
1569 (q0, q1, q2, q3, q4, q5)
1571 We can construct the (trivial) algebra of rank zero::
1573 sage: RealSymmetricEJA(0)
1574 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1578 def _denormalized_basis(cls
, n
, field
):
1580 Return a basis for the space of real symmetric n-by-n matrices.
1584 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1588 sage: set_random_seed()
1589 sage: n = ZZ.random_element(1,5)
1590 sage: B = RealSymmetricEJA._denormalized_basis(n,QQ)
1591 sage: all( M.is_symmetric() for M in B)
1595 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1599 for j
in range(i
+1):
1600 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
1604 Sij
= Eij
+ Eij
.transpose()
1610 def _max_random_instance_size():
1611 return 4 # Dimension 10
1614 def random_instance(cls
, field
=AA
, **kwargs
):
1616 Return a random instance of this type of algebra.
1618 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1619 return cls(n
, field
, **kwargs
)
1621 def __init__(self
, n
, field
=AA
, **kwargs
):
1622 basis
= self
._denormalized
_basis
(n
, field
)
1623 super(RealSymmetricEJA
, self
).__init
__(field
,
1625 self
.jordan_product
,
1626 self
.trace_inner_product
,
1628 self
.rank
.set_cache(n
)
1629 self
.one
.set_cache(self(matrix
.identity(field
,n
)))
1632 class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1636 Embed the n-by-n complex matrix ``M`` into the space of real
1637 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1638 bi` to the block matrix ``[[a,b],[-b,a]]``.
1642 sage: from mjo.eja.eja_algebra import \
1643 ....: ComplexMatrixEuclideanJordanAlgebra
1647 sage: F = QuadraticField(-1, 'I')
1648 sage: x1 = F(4 - 2*i)
1649 sage: x2 = F(1 + 2*i)
1652 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1653 sage: ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1662 Embedding is a homomorphism (isomorphism, in fact)::
1664 sage: set_random_seed()
1665 sage: n = ZZ.random_element(3)
1666 sage: F = QuadraticField(-1, 'I')
1667 sage: X = random_matrix(F, n)
1668 sage: Y = random_matrix(F, n)
1669 sage: Xe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X)
1670 sage: Ye = ComplexMatrixEuclideanJordanAlgebra.real_embed(Y)
1671 sage: XYe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1678 raise ValueError("the matrix 'M' must be square")
1680 # We don't need any adjoined elements...
1681 field
= M
.base_ring().base_ring()
1685 a
= z
.list()[0] # real part, I guess
1686 b
= z
.list()[1] # imag part, I guess
1687 blocks
.append(matrix(field
, 2, [[a
,b
],[-b
,a
]]))
1689 return matrix
.block(field
, n
, blocks
)
1693 def real_unembed(M
):
1695 The inverse of _embed_complex_matrix().
1699 sage: from mjo.eja.eja_algebra import \
1700 ....: ComplexMatrixEuclideanJordanAlgebra
1704 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1705 ....: [-2, 1, -4, 3],
1706 ....: [ 9, 10, 11, 12],
1707 ....: [-10, 9, -12, 11] ])
1708 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(A)
1710 [ 10*I + 9 12*I + 11]
1714 Unembedding is the inverse of embedding::
1716 sage: set_random_seed()
1717 sage: F = QuadraticField(-1, 'I')
1718 sage: M = random_matrix(F, 3)
1719 sage: Me = ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1720 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1726 raise ValueError("the matrix 'M' must be square")
1727 if not n
.mod(2).is_zero():
1728 raise ValueError("the matrix 'M' must be a complex embedding")
1730 # If "M" was normalized, its base ring might have roots
1731 # adjoined and they can stick around after unembedding.
1732 field
= M
.base_ring()
1733 R
= PolynomialRing(field
, 'z')
1736 # Sage doesn't know how to embed AA into QQbar, i.e. how
1737 # to adjoin sqrt(-1) to AA.
1740 F
= field
.extension(z
**2 + 1, 'I', embedding
=CLF(-1).sqrt())
1743 # Go top-left to bottom-right (reading order), converting every
1744 # 2-by-2 block we see to a single complex element.
1746 for k
in range(n
/2):
1747 for j
in range(n
/2):
1748 submat
= M
[2*k
:2*k
+2,2*j
:2*j
+2]
1749 if submat
[0,0] != submat
[1,1]:
1750 raise ValueError('bad on-diagonal submatrix')
1751 if submat
[0,1] != -submat
[1,0]:
1752 raise ValueError('bad off-diagonal submatrix')
1753 z
= submat
[0,0] + submat
[0,1]*i
1756 return matrix(F
, n
/2, elements
)
1760 def trace_inner_product(cls
,X
,Y
):
1762 Compute a matrix inner product in this algebra directly from
1767 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1771 This gives the same answer as the slow, default method implemented
1772 in :class:`MatrixEuclideanJordanAlgebra`::
1774 sage: set_random_seed()
1775 sage: J = ComplexHermitianEJA.random_instance()
1776 sage: x,y = J.random_elements(2)
1777 sage: Xe = x.to_matrix()
1778 sage: Ye = y.to_matrix()
1779 sage: X = ComplexHermitianEJA.real_unembed(Xe)
1780 sage: Y = ComplexHermitianEJA.real_unembed(Ye)
1781 sage: expected = (X*Y).trace().real()
1782 sage: actual = ComplexHermitianEJA.trace_inner_product(Xe,Ye)
1783 sage: actual == expected
1787 return RealMatrixEuclideanJordanAlgebra
.trace_inner_product(X
,Y
)/2
1790 class ComplexHermitianEJA(ConcreteEuclideanJordanAlgebra
,
1791 ComplexMatrixEuclideanJordanAlgebra
):
1793 The rank-n simple EJA consisting of complex Hermitian n-by-n
1794 matrices over the real numbers, the usual symmetric Jordan product,
1795 and the real-part-of-trace inner product. It has dimension `n^2` over
1800 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1804 In theory, our "field" can be any subfield of the reals::
1806 sage: ComplexHermitianEJA(2, RDF)
1807 Euclidean Jordan algebra of dimension 4 over Real Double Field
1808 sage: ComplexHermitianEJA(2, RR)
1809 Euclidean Jordan algebra of dimension 4 over Real Field with
1810 53 bits of precision
1814 The dimension of this algebra is `n^2`::
1816 sage: set_random_seed()
1817 sage: n_max = ComplexHermitianEJA._max_random_instance_size()
1818 sage: n = ZZ.random_element(1, n_max)
1819 sage: J = ComplexHermitianEJA(n)
1820 sage: J.dimension() == n^2
1823 The Jordan multiplication is what we think it is::
1825 sage: set_random_seed()
1826 sage: J = ComplexHermitianEJA.random_instance()
1827 sage: x,y = J.random_elements(2)
1828 sage: actual = (x*y).to_matrix()
1829 sage: X = x.to_matrix()
1830 sage: Y = y.to_matrix()
1831 sage: expected = (X*Y + Y*X)/2
1832 sage: actual == expected
1834 sage: J(expected) == x*y
1837 We can change the generator prefix::
1839 sage: ComplexHermitianEJA(2, prefix='z').gens()
1842 We can construct the (trivial) algebra of rank zero::
1844 sage: ComplexHermitianEJA(0)
1845 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1850 def _denormalized_basis(cls
, n
, field
):
1852 Returns a basis for the space of complex Hermitian n-by-n matrices.
1854 Why do we embed these? Basically, because all of numerical linear
1855 algebra assumes that you're working with vectors consisting of `n`
1856 entries from a field and scalars from the same field. There's no way
1857 to tell SageMath that (for example) the vectors contain complex
1858 numbers, while the scalar field is real.
1862 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1866 sage: set_random_seed()
1867 sage: n = ZZ.random_element(1,5)
1868 sage: field = QuadraticField(2, 'sqrt2')
1869 sage: B = ComplexHermitianEJA._denormalized_basis(n, field)
1870 sage: all( M.is_symmetric() for M in B)
1874 R
= PolynomialRing(field
, 'z')
1876 F
= field
.extension(z
**2 + 1, 'I')
1879 # This is like the symmetric case, but we need to be careful:
1881 # * We want conjugate-symmetry, not just symmetry.
1882 # * The diagonal will (as a result) be real.
1886 for j
in range(i
+1):
1887 Eij
= matrix(F
, n
, lambda k
,l
: k
==i
and l
==j
)
1889 Sij
= cls
.real_embed(Eij
)
1892 # The second one has a minus because it's conjugated.
1893 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
1895 Sij_imag
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
1898 # Since we embedded these, we can drop back to the "field" that we
1899 # started with instead of the complex extension "F".
1900 return tuple( s
.change_ring(field
) for s
in S
)
1903 def __init__(self
, n
, field
=AA
, **kwargs
):
1904 basis
= self
._denormalized
_basis
(n
,field
)
1905 super(ComplexHermitianEJA
, self
).__init
__(field
,
1907 self
.jordan_product
,
1908 self
.trace_inner_product
,
1910 self
.rank
.set_cache(n
)
1911 # TODO: pre-cache the identity!
1914 def _max_random_instance_size():
1915 return 3 # Dimension 9
1918 def random_instance(cls
, field
=AA
, **kwargs
):
1920 Return a random instance of this type of algebra.
1922 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1923 return cls(n
, field
, **kwargs
)
1925 class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1929 Embed the n-by-n quaternion matrix ``M`` into the space of real
1930 matrices of size 4n-by-4n by first sending each quaternion entry `z
1931 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
1932 c+di],[-c + di, a-bi]]`, and then embedding those into a real
1937 sage: from mjo.eja.eja_algebra import \
1938 ....: QuaternionMatrixEuclideanJordanAlgebra
1942 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1943 sage: i,j,k = Q.gens()
1944 sage: x = 1 + 2*i + 3*j + 4*k
1945 sage: M = matrix(Q, 1, [[x]])
1946 sage: QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1952 Embedding is a homomorphism (isomorphism, in fact)::
1954 sage: set_random_seed()
1955 sage: n = ZZ.random_element(2)
1956 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1957 sage: X = random_matrix(Q, n)
1958 sage: Y = random_matrix(Q, n)
1959 sage: Xe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X)
1960 sage: Ye = QuaternionMatrixEuclideanJordanAlgebra.real_embed(Y)
1961 sage: XYe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1966 quaternions
= M
.base_ring()
1969 raise ValueError("the matrix 'M' must be square")
1971 F
= QuadraticField(-1, 'I')
1976 t
= z
.coefficient_tuple()
1981 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
1982 [-c
+ d
*i
, a
- b
*i
]])
1983 realM
= ComplexMatrixEuclideanJordanAlgebra
.real_embed(cplxM
)
1984 blocks
.append(realM
)
1986 # We should have real entries by now, so use the realest field
1987 # we've got for the return value.
1988 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
1993 def real_unembed(M
):
1995 The inverse of _embed_quaternion_matrix().
1999 sage: from mjo.eja.eja_algebra import \
2000 ....: QuaternionMatrixEuclideanJordanAlgebra
2004 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
2005 ....: [-2, 1, -4, 3],
2006 ....: [-3, 4, 1, -2],
2007 ....: [-4, -3, 2, 1]])
2008 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(M)
2009 [1 + 2*i + 3*j + 4*k]
2013 Unembedding is the inverse of embedding::
2015 sage: set_random_seed()
2016 sage: Q = QuaternionAlgebra(QQ, -1, -1)
2017 sage: M = random_matrix(Q, 3)
2018 sage: Me = QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
2019 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
2025 raise ValueError("the matrix 'M' must be square")
2026 if not n
.mod(4).is_zero():
2027 raise ValueError("the matrix 'M' must be a quaternion embedding")
2029 # Use the base ring of the matrix to ensure that its entries can be
2030 # multiplied by elements of the quaternion algebra.
2031 field
= M
.base_ring()
2032 Q
= QuaternionAlgebra(field
,-1,-1)
2035 # Go top-left to bottom-right (reading order), converting every
2036 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
2039 for l
in range(n
/4):
2040 for m
in range(n
/4):
2041 submat
= ComplexMatrixEuclideanJordanAlgebra
.real_unembed(
2042 M
[4*l
:4*l
+4,4*m
:4*m
+4] )
2043 if submat
[0,0] != submat
[1,1].conjugate():
2044 raise ValueError('bad on-diagonal submatrix')
2045 if submat
[0,1] != -submat
[1,0].conjugate():
2046 raise ValueError('bad off-diagonal submatrix')
2047 z
= submat
[0,0].real()
2048 z
+= submat
[0,0].imag()*i
2049 z
+= submat
[0,1].real()*j
2050 z
+= submat
[0,1].imag()*k
2053 return matrix(Q
, n
/4, elements
)
2057 def trace_inner_product(cls
,X
,Y
):
2059 Compute a matrix inner product in this algebra directly from
2064 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2068 This gives the same answer as the slow, default method implemented
2069 in :class:`MatrixEuclideanJordanAlgebra`::
2071 sage: set_random_seed()
2072 sage: J = QuaternionHermitianEJA.random_instance()
2073 sage: x,y = J.random_elements(2)
2074 sage: Xe = x.to_matrix()
2075 sage: Ye = y.to_matrix()
2076 sage: X = QuaternionHermitianEJA.real_unembed(Xe)
2077 sage: Y = QuaternionHermitianEJA.real_unembed(Ye)
2078 sage: expected = (X*Y).trace().coefficient_tuple()[0]
2079 sage: actual = QuaternionHermitianEJA.trace_inner_product(Xe,Ye)
2080 sage: actual == expected
2084 return RealMatrixEuclideanJordanAlgebra
.trace_inner_product(X
,Y
)/4
2087 class QuaternionHermitianEJA(ConcreteEuclideanJordanAlgebra
,
2088 QuaternionMatrixEuclideanJordanAlgebra
):
2090 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
2091 matrices, the usual symmetric Jordan product, and the
2092 real-part-of-trace inner product. It has dimension `2n^2 - n` over
2097 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2101 In theory, our "field" can be any subfield of the reals::
2103 sage: QuaternionHermitianEJA(2, RDF)
2104 Euclidean Jordan algebra of dimension 6 over Real Double Field
2105 sage: QuaternionHermitianEJA(2, RR)
2106 Euclidean Jordan algebra of dimension 6 over Real Field with
2107 53 bits of precision
2111 The dimension of this algebra is `2*n^2 - n`::
2113 sage: set_random_seed()
2114 sage: n_max = QuaternionHermitianEJA._max_random_instance_size()
2115 sage: n = ZZ.random_element(1, n_max)
2116 sage: J = QuaternionHermitianEJA(n)
2117 sage: J.dimension() == 2*(n^2) - n
2120 The Jordan multiplication is what we think it is::
2122 sage: set_random_seed()
2123 sage: J = QuaternionHermitianEJA.random_instance()
2124 sage: x,y = J.random_elements(2)
2125 sage: actual = (x*y).to_matrix()
2126 sage: X = x.to_matrix()
2127 sage: Y = y.to_matrix()
2128 sage: expected = (X*Y + Y*X)/2
2129 sage: actual == expected
2131 sage: J(expected) == x*y
2134 We can change the generator prefix::
2136 sage: QuaternionHermitianEJA(2, prefix='a').gens()
2137 (a0, a1, a2, a3, a4, a5)
2139 We can construct the (trivial) algebra of rank zero::
2141 sage: QuaternionHermitianEJA(0)
2142 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2146 def _denormalized_basis(cls
, n
, field
):
2148 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
2150 Why do we embed these? Basically, because all of numerical
2151 linear algebra assumes that you're working with vectors consisting
2152 of `n` entries from a field and scalars from the same field. There's
2153 no way to tell SageMath that (for example) the vectors contain
2154 complex numbers, while the scalar field is real.
2158 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2162 sage: set_random_seed()
2163 sage: n = ZZ.random_element(1,5)
2164 sage: B = QuaternionHermitianEJA._denormalized_basis(n,QQ)
2165 sage: all( M.is_symmetric() for M in B )
2169 Q
= QuaternionAlgebra(QQ
,-1,-1)
2172 # This is like the symmetric case, but we need to be careful:
2174 # * We want conjugate-symmetry, not just symmetry.
2175 # * The diagonal will (as a result) be real.
2179 for j
in range(i
+1):
2180 Eij
= matrix(Q
, n
, lambda k
,l
: k
==i
and l
==j
)
2182 Sij
= cls
.real_embed(Eij
)
2185 # The second, third, and fourth ones have a minus
2186 # because they're conjugated.
2187 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
2189 Sij_I
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
2191 Sij_J
= cls
.real_embed(J
*Eij
- J
*Eij
.transpose())
2193 Sij_K
= cls
.real_embed(K
*Eij
- K
*Eij
.transpose())
2196 # Since we embedded these, we can drop back to the "field" that we
2197 # started with instead of the quaternion algebra "Q".
2198 return tuple( s
.change_ring(field
) for s
in S
)
2201 def __init__(self
, n
, field
=AA
, **kwargs
):
2202 basis
= self
._denormalized
_basis
(n
,field
)
2203 super(QuaternionHermitianEJA
, self
).__init
__(field
,
2205 self
.jordan_product
,
2206 self
.trace_inner_product
,
2208 self
.rank
.set_cache(n
)
2209 # TODO: cache one()!
2212 def _max_random_instance_size():
2214 The maximum rank of a random QuaternionHermitianEJA.
2216 return 2 # Dimension 6
2219 def random_instance(cls
, field
=AA
, **kwargs
):
2221 Return a random instance of this type of algebra.
2223 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2224 return cls(n
, field
, **kwargs
)
2227 class HadamardEJA(ConcreteEuclideanJordanAlgebra
):
2229 Return the Euclidean Jordan Algebra corresponding to the set
2230 `R^n` under the Hadamard product.
2232 Note: this is nothing more than the Cartesian product of ``n``
2233 copies of the spin algebra. Once Cartesian product algebras
2234 are implemented, this can go.
2238 sage: from mjo.eja.eja_algebra import HadamardEJA
2242 This multiplication table can be verified by hand::
2244 sage: J = HadamardEJA(3)
2245 sage: e0,e1,e2 = J.gens()
2261 We can change the generator prefix::
2263 sage: HadamardEJA(3, prefix='r').gens()
2267 def __init__(self
, n
, field
=AA
, **kwargs
):
2268 V
= VectorSpace(field
, n
)
2271 def jordan_product(x
,y
):
2272 return V([ xi
*yi
for (xi
,yi
) in zip(x
,y
) ])
2273 def inner_product(x
,y
):
2274 return x
.inner_product(y
)
2276 super(HadamardEJA
, self
).__init
__(field
,
2281 self
.rank
.set_cache(n
)
2284 self
.one
.set_cache( self
.zero() )
2286 self
.one
.set_cache( sum(self
.gens()) )
2289 def _max_random_instance_size():
2291 The maximum dimension of a random HadamardEJA.
2296 def random_instance(cls
, field
=AA
, **kwargs
):
2298 Return a random instance of this type of algebra.
2300 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2301 return cls(n
, field
, **kwargs
)
2304 class BilinearFormEJA(ConcreteEuclideanJordanAlgebra
):
2306 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2307 with the half-trace inner product and jordan product ``x*y =
2308 (<Bx,y>,y_bar>, x0*y_bar + y0*x_bar)`` where `B = 1 \times B22` is
2309 a symmetric positive-definite "bilinear form" matrix. Its
2310 dimension is the size of `B`, and it has rank two in dimensions
2311 larger than two. It reduces to the ``JordanSpinEJA`` when `B` is
2312 the identity matrix of order ``n``.
2314 We insist that the one-by-one upper-left identity block of `B` be
2315 passed in as well so that we can be passed a matrix of size zero
2316 to construct a trivial algebra.
2320 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
2321 ....: JordanSpinEJA)
2325 When no bilinear form is specified, the identity matrix is used,
2326 and the resulting algebra is the Jordan spin algebra::
2328 sage: B = matrix.identity(AA,3)
2329 sage: J0 = BilinearFormEJA(B)
2330 sage: J1 = JordanSpinEJA(3)
2331 sage: J0.multiplication_table() == J0.multiplication_table()
2334 An error is raised if the matrix `B` does not correspond to a
2335 positive-definite bilinear form::
2337 sage: B = matrix.random(QQ,2,3)
2338 sage: J = BilinearFormEJA(B)
2339 Traceback (most recent call last):
2341 ValueError: bilinear form is not positive-definite
2342 sage: B = matrix.zero(QQ,3)
2343 sage: J = BilinearFormEJA(B)
2344 Traceback (most recent call last):
2346 ValueError: bilinear form is not positive-definite
2350 We can create a zero-dimensional algebra::
2352 sage: B = matrix.identity(AA,0)
2353 sage: J = BilinearFormEJA(B)
2357 We can check the multiplication condition given in the Jordan, von
2358 Neumann, and Wigner paper (and also discussed on my "On the
2359 symmetry..." paper). Note that this relies heavily on the standard
2360 choice of basis, as does anything utilizing the bilinear form
2361 matrix. We opt not to orthonormalize the basis, because if we
2362 did, we would have to normalize the `s_{i}` in a similar manner::
2364 sage: set_random_seed()
2365 sage: n = ZZ.random_element(5)
2366 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2367 sage: B11 = matrix.identity(QQ,1)
2368 sage: B22 = M.transpose()*M
2369 sage: B = block_matrix(2,2,[ [B11,0 ],
2371 sage: J = BilinearFormEJA(B, orthonormalize=False)
2372 sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
2373 sage: V = J.vector_space()
2374 sage: sis = [ J( V([0] + (M.inverse()*ei).list()).column() )
2375 ....: for ei in eis ]
2376 sage: actual = [ sis[i]*sis[j]
2377 ....: for i in range(n-1)
2378 ....: for j in range(n-1) ]
2379 sage: expected = [ J.one() if i == j else J.zero()
2380 ....: for i in range(n-1)
2381 ....: for j in range(n-1) ]
2382 sage: actual == expected
2385 def __init__(self
, B
, field
=AA
, **kwargs
):
2386 if not B
.is_positive_definite():
2387 raise ValueError("bilinear form is not positive-definite")
2390 V
= VectorSpace(field
, n
)
2392 def inner_product(x
,y
):
2393 return (B
*x
).inner_product(y
)
2395 def jordan_product(x
,y
):
2400 z0
= inner_product(x
,y
)
2401 zbar
= y0
*xbar
+ x0
*ybar
2402 return V([z0
] + zbar
.list())
2404 super(BilinearFormEJA
, self
).__init
__(field
,
2410 # The rank of this algebra is two, unless we're in a
2411 # one-dimensional ambient space (because the rank is bounded
2412 # by the ambient dimension).
2413 self
.rank
.set_cache(min(n
,2))
2416 self
.one
.set_cache( self
.zero() )
2418 self
.one
.set_cache( self
.monomial(0) )
2421 def _max_random_instance_size():
2423 The maximum dimension of a random BilinearFormEJA.
2428 def random_instance(cls
, field
=AA
, **kwargs
):
2430 Return a random instance of this algebra.
2432 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2434 B
= matrix
.identity(field
, n
)
2435 return cls(B
, field
, **kwargs
)
2437 B11
= matrix
.identity(field
,1)
2438 M
= matrix
.random(field
, n
-1)
2439 I
= matrix
.identity(field
, n
-1)
2440 alpha
= field
.zero()
2441 while alpha
.is_zero():
2442 alpha
= field
.random_element().abs()
2443 B22
= M
.transpose()*M
+ alpha
*I
2445 from sage
.matrix
.special
import block_matrix
2446 B
= block_matrix(2,2, [ [B11
, ZZ(0) ],
2449 return cls(B
, field
, **kwargs
)
2452 class JordanSpinEJA(BilinearFormEJA
):
2454 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2455 with the usual inner product and jordan product ``x*y =
2456 (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2461 sage: from mjo.eja.eja_algebra import JordanSpinEJA
2465 This multiplication table can be verified by hand::
2467 sage: J = JordanSpinEJA(4)
2468 sage: e0,e1,e2,e3 = J.gens()
2484 We can change the generator prefix::
2486 sage: JordanSpinEJA(2, prefix='B').gens()
2491 Ensure that we have the usual inner product on `R^n`::
2493 sage: set_random_seed()
2494 sage: J = JordanSpinEJA.random_instance()
2495 sage: x,y = J.random_elements(2)
2496 sage: actual = x.inner_product(y)
2497 sage: expected = x.to_vector().inner_product(y.to_vector())
2498 sage: actual == expected
2502 def __init__(self
, n
, field
=AA
, **kwargs
):
2503 # This is a special case of the BilinearFormEJA with the identity
2504 # matrix as its bilinear form.
2505 B
= matrix
.identity(field
, n
)
2506 super(JordanSpinEJA
, self
).__init
__(B
, field
, **kwargs
)
2509 def _max_random_instance_size():
2511 The maximum dimension of a random JordanSpinEJA.
2516 def random_instance(cls
, field
=AA
, **kwargs
):
2518 Return a random instance of this type of algebra.
2520 Needed here to override the implementation for ``BilinearFormEJA``.
2522 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2523 return cls(n
, field
, **kwargs
)
2526 class TrivialEJA(ConcreteEuclideanJordanAlgebra
):
2528 The trivial Euclidean Jordan algebra consisting of only a zero element.
2532 sage: from mjo.eja.eja_algebra import TrivialEJA
2536 sage: J = TrivialEJA()
2543 sage: 7*J.one()*12*J.one()
2545 sage: J.one().inner_product(J.one())
2547 sage: J.one().norm()
2549 sage: J.one().subalgebra_generated_by()
2550 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2555 def __init__(self
, field
=AA
, **kwargs
):
2556 jordan_product
= lambda x
,y
: x
2557 inner_product
= lambda x
,y
: field(0)
2559 super(TrivialEJA
, self
).__init
__(field
,
2564 # The rank is zero using my definition, namely the dimension of the
2565 # largest subalgebra generated by any element.
2566 self
.rank
.set_cache(0)
2567 self
.one
.set_cache( self
.zero() )
2570 def random_instance(cls
, field
=AA
, **kwargs
):
2571 # We don't take a "size" argument so the superclass method is
2572 # inappropriate for us.
2573 return cls(field
, **kwargs
)
2575 class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra
):
2577 The external (orthogonal) direct sum of two other Euclidean Jordan
2578 algebras. Essentially the Cartesian product of its two factors.
2579 Every Euclidean Jordan algebra decomposes into an orthogonal
2580 direct sum of simple Euclidean Jordan algebras, so no generality
2581 is lost by providing only this construction.
2585 sage: from mjo.eja.eja_algebra import (random_eja,
2587 ....: RealSymmetricEJA,
2592 sage: J1 = HadamardEJA(2)
2593 sage: J2 = RealSymmetricEJA(3)
2594 sage: J = DirectSumEJA(J1,J2)
2602 The external direct sum construction is only valid when the two factors
2603 have the same base ring; an error is raised otherwise::
2605 sage: set_random_seed()
2606 sage: J1 = random_eja(AA)
2607 sage: J2 = random_eja(QQ,orthonormalize=False)
2608 sage: J = DirectSumEJA(J1,J2)
2609 Traceback (most recent call last):
2611 ValueError: algebras must share the same base field
2614 def __init__(self
, J1
, J2
, **kwargs
):
2615 if J1
.base_ring() != J2
.base_ring():
2616 raise ValueError("algebras must share the same base field")
2617 field
= J1
.base_ring()
2619 self
._factors
= (J1
, J2
)
2623 V
= VectorSpace(field
, n
)
2624 mult_table
= [ [ V
.zero() for j
in range(i
+1) ]
2627 for j
in range(i
+1):
2628 p
= (J1
.monomial(i
)*J1
.monomial(j
)).to_vector()
2629 mult_table
[i
][j
] = V(p
.list() + [field
.zero()]*n2
)
2632 for j
in range(i
+1):
2633 p
= (J2
.monomial(i
)*J2
.monomial(j
)).to_vector()
2634 mult_table
[n1
+i
][n1
+j
] = V([field
.zero()]*n1
+ p
.list())
2636 # TODO: build the IP table here from the two constituent IP
2637 # matrices (it'll be block diagonal, I think).
2638 ip_table
= [ [ field
.zero() for j
in range(i
+1) ]
2640 super(DirectSumEJA
, self
).__init
__(field
,
2645 self
.rank
.set_cache(J1
.rank() + J2
.rank())
2650 Return the pair of this algebra's factors.
2654 sage: from mjo.eja.eja_algebra import (HadamardEJA,
2655 ....: JordanSpinEJA,
2660 sage: J1 = HadamardEJA(2,QQ)
2661 sage: J2 = JordanSpinEJA(3,QQ)
2662 sage: J = DirectSumEJA(J1,J2)
2664 (Euclidean Jordan algebra of dimension 2 over Rational Field,
2665 Euclidean Jordan algebra of dimension 3 over Rational Field)
2668 return self
._factors
2670 def projections(self
):
2672 Return a pair of projections onto this algebra's factors.
2676 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
2677 ....: ComplexHermitianEJA,
2682 sage: J1 = JordanSpinEJA(2)
2683 sage: J2 = ComplexHermitianEJA(2)
2684 sage: J = DirectSumEJA(J1,J2)
2685 sage: (pi_left, pi_right) = J.projections()
2686 sage: J.one().to_vector()
2688 sage: pi_left(J.one()).to_vector()
2690 sage: pi_right(J.one()).to_vector()
2694 (J1
,J2
) = self
.factors()
2697 V_basis
= self
.vector_space().basis()
2698 # Need to specify the dimensions explicitly so that we don't
2699 # wind up with a zero-by-zero matrix when we want e.g. a
2700 # zero-by-two matrix (important for composing things).
2701 P1
= matrix(self
.base_ring(), m
, m
+n
, V_basis
[:m
])
2702 P2
= matrix(self
.base_ring(), n
, m
+n
, V_basis
[m
:])
2703 pi_left
= FiniteDimensionalEuclideanJordanAlgebraOperator(self
,J1
,P1
)
2704 pi_right
= FiniteDimensionalEuclideanJordanAlgebraOperator(self
,J2
,P2
)
2705 return (pi_left
, pi_right
)
2707 def inclusions(self
):
2709 Return the pair of inclusion maps from our factors into us.
2713 sage: from mjo.eja.eja_algebra import (random_eja,
2714 ....: JordanSpinEJA,
2715 ....: RealSymmetricEJA,
2720 sage: J1 = JordanSpinEJA(3)
2721 sage: J2 = RealSymmetricEJA(2)
2722 sage: J = DirectSumEJA(J1,J2)
2723 sage: (iota_left, iota_right) = J.inclusions()
2724 sage: iota_left(J1.zero()) == J.zero()
2726 sage: iota_right(J2.zero()) == J.zero()
2728 sage: J1.one().to_vector()
2730 sage: iota_left(J1.one()).to_vector()
2732 sage: J2.one().to_vector()
2734 sage: iota_right(J2.one()).to_vector()
2736 sage: J.one().to_vector()
2741 Composing a projection with the corresponding inclusion should
2742 produce the identity map, and mismatching them should produce
2745 sage: set_random_seed()
2746 sage: J1 = random_eja()
2747 sage: J2 = random_eja()
2748 sage: J = DirectSumEJA(J1,J2)
2749 sage: (iota_left, iota_right) = J.inclusions()
2750 sage: (pi_left, pi_right) = J.projections()
2751 sage: pi_left*iota_left == J1.one().operator()
2753 sage: pi_right*iota_right == J2.one().operator()
2755 sage: (pi_left*iota_right).is_zero()
2757 sage: (pi_right*iota_left).is_zero()
2761 (J1
,J2
) = self
.factors()
2764 V_basis
= self
.vector_space().basis()
2765 # Need to specify the dimensions explicitly so that we don't
2766 # wind up with a zero-by-zero matrix when we want e.g. a
2767 # two-by-zero matrix (important for composing things).
2768 I1
= matrix
.column(self
.base_ring(), m
, m
+n
, V_basis
[:m
])
2769 I2
= matrix
.column(self
.base_ring(), n
, m
+n
, V_basis
[m
:])
2770 iota_left
= FiniteDimensionalEuclideanJordanAlgebraOperator(J1
,self
,I1
)
2771 iota_right
= FiniteDimensionalEuclideanJordanAlgebraOperator(J2
,self
,I2
)
2772 return (iota_left
, iota_right
)
2774 def inner_product(self
, x
, y
):
2776 The standard Cartesian inner-product.
2778 We project ``x`` and ``y`` onto our factors, and add up the
2779 inner-products from the subalgebras.
2784 sage: from mjo.eja.eja_algebra import (HadamardEJA,
2785 ....: QuaternionHermitianEJA,
2790 sage: J1 = HadamardEJA(3,QQ)
2791 sage: J2 = QuaternionHermitianEJA(2,QQ,orthonormalize=False)
2792 sage: J = DirectSumEJA(J1,J2)
2797 sage: x1.inner_product(x2)
2799 sage: y1.inner_product(y2)
2801 sage: J.one().inner_product(J.one())
2805 (pi_left
, pi_right
) = self
.projections()
2811 return (x1
.inner_product(y1
) + x2
.inner_product(y2
))
2815 random_eja
= ConcreteEuclideanJordanAlgebra
.random_instance