2 Euclidean Jordan Algebras. These are formally-real Jordan Algebras;
3 specifically those where u^2 + v^2 = 0 implies that u = v = 0. They
4 are used in optimization, and have some additional nice methods beyond
5 what can be supported in a general Jordan Algebra.
10 sage: from mjo.eja.eja_algebra import random_eja
15 Euclidean Jordan algebra of dimension...
19 from itertools
import repeat
21 from sage
.algebras
.quatalg
.quaternion_algebra
import QuaternionAlgebra
22 from sage
.categories
.magmatic_algebras
import MagmaticAlgebras
23 from sage
.combinat
.free_module
import CombinatorialFreeModule
24 from sage
.matrix
.constructor
import matrix
25 from sage
.matrix
.matrix_space
import MatrixSpace
26 from sage
.misc
.cachefunc
import cached_method
27 from sage
.misc
.table
import table
28 from sage
.modules
.free_module
import FreeModule
, VectorSpace
29 from sage
.rings
.all
import (ZZ
, QQ
, AA
, QQbar
, RR
, RLF
, CLF
,
32 from mjo
.eja
.eja_element
import FiniteDimensionalEuclideanJordanAlgebraElement
33 from mjo
.eja
.eja_operator
import FiniteDimensionalEuclideanJordanAlgebraOperator
34 from mjo
.eja
.eja_utils
import _mat2vec
36 class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule
):
38 The lowest-level class for representing a Euclidean Jordan algebra.
40 def _coerce_map_from_base_ring(self
):
42 Disable the map from the base ring into the algebra.
44 Performing a nonsense conversion like this automatically
45 is counterpedagogical. The fallback is to try the usual
46 element constructor, which should also fail.
50 sage: from mjo.eja.eja_algebra import random_eja
54 sage: set_random_seed()
55 sage: J = random_eja()
57 Traceback (most recent call last):
59 ValueError: not an element of this algebra
76 * field -- the scalar field for this algebra (must be real)
78 * multiplication_table -- the multiplication table for this
79 algebra's implicit basis. Only the lower-triangular portion
80 of the table is used, since the multiplication is assumed
85 sage: from mjo.eja.eja_algebra import (
86 ....: FiniteDimensionalEuclideanJordanAlgebra,
92 By definition, Jordan multiplication commutes::
94 sage: set_random_seed()
95 sage: J = random_eja()
96 sage: x,y = J.random_elements(2)
100 An error is raised if the Jordan product is not commutative::
102 sage: JP = ((1,2),(0,0))
103 sage: IP = ((1,0),(0,1))
104 sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,JP,IP)
105 Traceback (most recent call last):
107 ValueError: Jordan product is not commutative
109 An error is raised if the inner-product is not commutative::
111 sage: JP = ((1,0),(0,1))
112 sage: IP = ((1,2),(0,0))
113 sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,JP,IP)
114 Traceback (most recent call last):
116 ValueError: inner-product is not commutative
120 The ``field`` we're given must be real with ``check_field=True``::
122 sage: JordanSpinEJA(2,QQbar)
123 Traceback (most recent call last):
125 ValueError: scalar field is not real
127 The multiplication table must be square with ``check_axioms=True``::
129 sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,((),()),((1,),))
130 Traceback (most recent call last):
132 ValueError: multiplication table is not square
134 The multiplication and inner-product tables must be the same
135 size (and in particular, the inner-product table must also be
136 square) with ``check_axioms=True``::
138 sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,((1,),),(()))
139 Traceback (most recent call last):
141 ValueError: multiplication and inner-product tables are
143 sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,((1,),),((1,2),))
144 Traceback (most recent call last):
146 ValueError: multiplication and inner-product tables are
151 if not field
.is_subring(RR
):
152 # Note: this does return true for the real algebraic
153 # field, the rationals, and any quadratic field where
154 # we've specified a real embedding.
155 raise ValueError("scalar field is not real")
158 # The multiplication and inner-product tables should be square
159 # if the user wants us to verify them. And we verify them as
160 # soon as possible, because we want to exploit their symmetry.
161 n
= len(multiplication_table
)
163 if not all( len(l
) == n
for l
in multiplication_table
):
164 raise ValueError("multiplication table is not square")
166 # If the multiplication table is square, we can check if
167 # the inner-product table is square by comparing it to the
168 # multiplication table's dimensions.
169 msg
= "multiplication and inner-product tables are different sizes"
170 if not len(inner_product_table
) == n
:
171 raise ValueError(msg
)
173 if not all( len(l
) == n
for l
in inner_product_table
):
174 raise ValueError(msg
)
176 if not all( multiplication_table
[j
][i
]
177 == multiplication_table
[i
][j
]
179 for j
in range(i
+1) ):
180 raise ValueError("Jordan product is not commutative")
181 if not all( inner_product_table
[j
][i
]
182 == inner_product_table
[i
][j
]
184 for j
in range(i
+1) ):
185 raise ValueError("inner-product is not commutative")
186 self
._matrix
_basis
= matrix_basis
189 category
= MagmaticAlgebras(field
).FiniteDimensional()
190 category
= category
.WithBasis().Unital()
192 fda
= super(FiniteDimensionalEuclideanJordanAlgebra
, self
)
197 self
.print_options(bracket
='')
199 # The multiplication table we're given is necessarily in terms
200 # of vectors, because we don't have an algebra yet for
201 # anything to be an element of. However, it's faster in the
202 # long run to have the multiplication table be in terms of
203 # algebra elements. We do this after calling the superclass
204 # constructor so that from_vector() knows what to do.
206 # Note: we take advantage of symmetry here, and only store
207 # the lower-triangular portion of the table.
208 self
._multiplication
_table
= [ [ self
.vector_space().zero()
209 for j
in range(i
+1) ]
214 elt
= self
.from_vector(multiplication_table
[i
][j
])
215 self
._multiplication
_table
[i
][j
] = elt
217 self
._multiplication
_table
= tuple(map(tuple, self
._multiplication
_table
))
219 # Save our inner product as a matrix, since the efficiency of
220 # matrix multiplication will usually outweigh the fact that we
221 # have to store a redundant upper- or lower-triangular part.
222 # Pre-cache the fact that these are Hermitian (real symmetric,
223 # in fact) in case some e.g. matrix multiplication routine can
224 # take advantage of it.
225 self
._inner
_product
_matrix
= matrix(field
, inner_product_table
)
226 self
._inner
_product
_matrix
._cache
= {'hermitian': True}
227 self
._inner
_product
_matrix
.set_immutable()
230 if not self
._is
_jordanian
():
231 raise ValueError("Jordan identity does not hold")
232 if not self
._inner
_product
_is
_associative
():
233 raise ValueError("inner product is not associative")
235 def _element_constructor_(self
, elt
):
237 Construct an element of this algebra from its vector or matrix
240 This gets called only after the parent element _call_ method
241 fails to find a coercion for the argument.
245 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
247 ....: RealSymmetricEJA)
251 The identity in `S^n` is converted to the identity in the EJA::
253 sage: J = RealSymmetricEJA(3)
254 sage: I = matrix.identity(QQ,3)
255 sage: J(I) == J.one()
258 This skew-symmetric matrix can't be represented in the EJA::
260 sage: J = RealSymmetricEJA(3)
261 sage: A = matrix(QQ,3, lambda i,j: i-j)
263 Traceback (most recent call last):
265 ValueError: not an element of this algebra
269 Ensure that we can convert any element of the two non-matrix
270 simple algebras (whose matrix representations are columns)
271 back and forth faithfully::
273 sage: set_random_seed()
274 sage: J = HadamardEJA.random_instance()
275 sage: x = J.random_element()
276 sage: J(x.to_vector().column()) == x
278 sage: J = JordanSpinEJA.random_instance()
279 sage: x = J.random_element()
280 sage: J(x.to_vector().column()) == x
283 msg
= "not an element of this algebra"
285 # The superclass implementation of random_element()
286 # needs to be able to coerce "0" into the algebra.
288 elif elt
in self
.base_ring():
289 # Ensure that no base ring -> algebra coercion is performed
290 # by this method. There's some stupidity in sage that would
291 # otherwise propagate to this method; for example, sage thinks
292 # that the integer 3 belongs to the space of 2-by-2 matrices.
293 raise ValueError(msg
)
295 if elt
not in self
.matrix_space():
296 raise ValueError(msg
)
298 # Thanks for nothing! Matrix spaces aren't vector spaces in
299 # Sage, so we have to figure out its matrix-basis coordinates
300 # ourselves. We use the basis space's ring instead of the
301 # element's ring because the basis space might be an algebraic
302 # closure whereas the base ring of the 3-by-3 identity matrix
303 # could be QQ instead of QQbar.
304 V
= VectorSpace(self
.base_ring(), elt
.nrows()*elt
.ncols())
305 W
= V
.span_of_basis( _mat2vec(s
) for s
in self
.matrix_basis() )
308 coords
= W
.coordinate_vector(_mat2vec(elt
))
309 except ArithmeticError: # vector is not in free module
310 raise ValueError(msg
)
312 return self
.from_vector(coords
)
316 Return a string representation of ``self``.
320 sage: from mjo.eja.eja_algebra import JordanSpinEJA
324 Ensure that it says what we think it says::
326 sage: JordanSpinEJA(2, field=AA)
327 Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
328 sage: JordanSpinEJA(3, field=RDF)
329 Euclidean Jordan algebra of dimension 3 over Real Double Field
332 fmt
= "Euclidean Jordan algebra of dimension {} over {}"
333 return fmt
.format(self
.dimension(), self
.base_ring())
335 def product_on_basis(self
, i
, j
):
336 # We only stored the lower-triangular portion of the
337 # multiplication table.
339 return self
._multiplication
_table
[i
][j
]
341 return self
._multiplication
_table
[j
][i
]
343 def _is_commutative(self
):
345 Whether or not this algebra's multiplication table is commutative.
347 This method should of course always return ``True``, unless
348 this algebra was constructed with ``check_axioms=False`` and
349 passed an invalid multiplication table.
351 return all( self
.product_on_basis(i
,j
) == self
.product_on_basis(i
,j
)
352 for i
in range(self
.dimension())
353 for j
in range(self
.dimension()) )
355 def _is_jordanian(self
):
357 Whether or not this algebra's multiplication table respects the
358 Jordan identity `(x^{2})(xy) = x(x^{2}y)`.
360 We only check one arrangement of `x` and `y`, so for a
361 ``True`` result to be truly true, you should also check
362 :meth:`_is_commutative`. This method should of course always
363 return ``True``, unless this algebra was constructed with
364 ``check_axioms=False`` and passed an invalid multiplication table.
366 return all( (self
.monomial(i
)**2)*(self
.monomial(i
)*self
.monomial(j
))
368 (self
.monomial(i
))*((self
.monomial(i
)**2)*self
.monomial(j
))
369 for i
in range(self
.dimension())
370 for j
in range(self
.dimension()) )
372 def _inner_product_is_associative(self
):
374 Return whether or not this algebra's inner product `B` is
375 associative; that is, whether or not `B(xy,z) = B(x,yz)`.
377 This method should of course always return ``True``, unless
378 this algebra was constructed with ``check_axioms=False`` and
379 passed an invalid multiplication table.
382 # Used to check whether or not something is zero in an inexact
383 # ring. This number is sufficient to allow the construction of
384 # QuaternionHermitianEJA(2, RDF) with check_axioms=True.
387 for i
in range(self
.dimension()):
388 for j
in range(self
.dimension()):
389 for k
in range(self
.dimension()):
393 diff
= (x
*y
).inner_product(z
) - x
.inner_product(y
*z
)
395 if self
.base_ring().is_exact():
399 if diff
.abs() > epsilon
:
405 def characteristic_polynomial_of(self
):
407 Return the algebra's "characteristic polynomial of" function,
408 which is itself a multivariate polynomial that, when evaluated
409 at the coordinates of some algebra element, returns that
410 element's characteristic polynomial.
412 The resulting polynomial has `n+1` variables, where `n` is the
413 dimension of this algebra. The first `n` variables correspond to
414 the coordinates of an algebra element: when evaluated at the
415 coordinates of an algebra element with respect to a certain
416 basis, the result is a univariate polynomial (in the one
417 remaining variable ``t``), namely the characteristic polynomial
422 sage: from mjo.eja.eja_algebra import JordanSpinEJA, TrivialEJA
426 The characteristic polynomial in the spin algebra is given in
427 Alizadeh, Example 11.11::
429 sage: J = JordanSpinEJA(3)
430 sage: p = J.characteristic_polynomial_of(); p
431 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
432 sage: xvec = J.one().to_vector()
436 By definition, the characteristic polynomial is a monic
437 degree-zero polynomial in a rank-zero algebra. Note that
438 Cayley-Hamilton is indeed satisfied since the polynomial
439 ``1`` evaluates to the identity element of the algebra on
442 sage: J = TrivialEJA()
443 sage: J.characteristic_polynomial_of()
450 # The list of coefficient polynomials a_0, a_1, a_2, ..., a_(r-1).
451 a
= self
._charpoly
_coefficients
()
453 # We go to a bit of trouble here to reorder the
454 # indeterminates, so that it's easier to evaluate the
455 # characteristic polynomial at x's coordinates and get back
456 # something in terms of t, which is what we want.
457 S
= PolynomialRing(self
.base_ring(),'t')
461 S
= PolynomialRing(S
, R
.variable_names())
464 return (t
**r
+ sum( a
[k
]*(t
**k
) for k
in range(r
) ))
466 def coordinate_polynomial_ring(self
):
468 The multivariate polynomial ring in which this algebra's
469 :meth:`characteristic_polynomial_of` lives.
473 sage: from mjo.eja.eja_algebra import (HadamardEJA,
474 ....: RealSymmetricEJA)
478 sage: J = HadamardEJA(2)
479 sage: J.coordinate_polynomial_ring()
480 Multivariate Polynomial Ring in X1, X2...
481 sage: J = RealSymmetricEJA(3,QQ,orthonormalize=False)
482 sage: J.coordinate_polynomial_ring()
483 Multivariate Polynomial Ring in X1, X2, X3, X4, X5, X6...
486 var_names
= tuple( "X%d" % z
for z
in range(1, self
.dimension()+1) )
487 return PolynomialRing(self
.base_ring(), var_names
)
489 def inner_product(self
, x
, y
):
491 The inner product associated with this Euclidean Jordan algebra.
493 Defaults to the trace inner product, but can be overridden by
494 subclasses if they are sure that the necessary properties are
499 sage: from mjo.eja.eja_algebra import (random_eja,
501 ....: BilinearFormEJA)
505 Our inner product is "associative," which means the following for
506 a symmetric bilinear form::
508 sage: set_random_seed()
509 sage: J = random_eja()
510 sage: x,y,z = J.random_elements(3)
511 sage: (x*y).inner_product(z) == y.inner_product(x*z)
516 Ensure that this is the usual inner product for the algebras
519 sage: set_random_seed()
520 sage: J = HadamardEJA.random_instance()
521 sage: x,y = J.random_elements(2)
522 sage: actual = x.inner_product(y)
523 sage: expected = x.to_vector().inner_product(y.to_vector())
524 sage: actual == expected
527 Ensure that this is one-half of the trace inner-product in a
528 BilinearFormEJA that isn't just the reals (when ``n`` isn't
529 one). This is in Faraut and Koranyi, and also my "On the
532 sage: set_random_seed()
533 sage: J = BilinearFormEJA.random_instance()
534 sage: n = J.dimension()
535 sage: x = J.random_element()
536 sage: y = J.random_element()
537 sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
540 B
= self
._inner
_product
_matrix
541 return (B
*x
.to_vector()).inner_product(y
.to_vector())
544 def is_trivial(self
):
546 Return whether or not this algebra is trivial.
548 A trivial algebra contains only the zero element.
552 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
557 sage: J = ComplexHermitianEJA(3)
563 sage: J = TrivialEJA()
568 return self
.dimension() == 0
571 def multiplication_table(self
):
573 Return a visual representation of this algebra's multiplication
574 table (on basis elements).
578 sage: from mjo.eja.eja_algebra import JordanSpinEJA
582 sage: J = JordanSpinEJA(4)
583 sage: J.multiplication_table()
584 +----++----+----+----+----+
585 | * || e0 | e1 | e2 | e3 |
586 +====++====+====+====+====+
587 | e0 || e0 | e1 | e2 | e3 |
588 +----++----+----+----+----+
589 | e1 || e1 | e0 | 0 | 0 |
590 +----++----+----+----+----+
591 | e2 || e2 | 0 | e0 | 0 |
592 +----++----+----+----+----+
593 | e3 || e3 | 0 | 0 | e0 |
594 +----++----+----+----+----+
598 M
= [ [ self
.zero() for j
in range(n
) ]
602 M
[i
][j
] = self
._multiplication
_table
[i
][j
]
606 # Prepend the left "header" column entry Can't do this in
607 # the loop because it messes up the symmetry.
608 M
[i
] = [self
.monomial(i
)] + M
[i
]
610 # Prepend the header row.
611 M
= [["*"] + list(self
.gens())] + M
612 return table(M
, header_row
=True, header_column
=True, frame
=True)
615 def matrix_basis(self
):
617 Return an (often more natural) representation of this algebras
618 basis as an ordered tuple of matrices.
620 Every finite-dimensional Euclidean Jordan Algebra is a, up to
621 Jordan isomorphism, a direct sum of five simple
622 algebras---four of which comprise Hermitian matrices. And the
623 last type of algebra can of course be thought of as `n`-by-`1`
624 column matrices (ambiguusly called column vectors) to avoid
625 special cases. As a result, matrices (and column vectors) are
626 a natural representation format for Euclidean Jordan algebra
629 But, when we construct an algebra from a basis of matrices,
630 those matrix representations are lost in favor of coordinate
631 vectors *with respect to* that basis. We could eventually
632 convert back if we tried hard enough, but having the original
633 representations handy is valuable enough that we simply store
634 them and return them from this method.
636 Why implement this for non-matrix algebras? Avoiding special
637 cases for the :class:`BilinearFormEJA` pays with simplicity in
638 its own right. But mainly, we would like to be able to assume
639 that elements of a :class:`DirectSumEJA` can be displayed
640 nicely, without having to have special classes for direct sums
641 one of whose components was a matrix algebra.
645 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
646 ....: RealSymmetricEJA)
650 sage: J = RealSymmetricEJA(2)
652 Finite family {0: e0, 1: e1, 2: e2}
653 sage: J.matrix_basis()
655 [1 0] [ 0 0.7071067811865475?] [0 0]
656 [0 0], [0.7071067811865475? 0], [0 1]
661 sage: J = JordanSpinEJA(2)
663 Finite family {0: e0, 1: e1}
664 sage: J.matrix_basis()
670 if self
._matrix
_basis
is None:
671 M
= self
.matrix_space()
672 return tuple( M(b
.to_vector()) for b
in self
.basis() )
674 return self
._matrix
_basis
677 def matrix_space(self
):
679 Return the matrix space in which this algebra's elements live, if
680 we think of them as matrices (including column vectors of the
683 Generally this will be an `n`-by-`1` column-vector space,
684 except when the algebra is trivial. There it's `n`-by-`n`
685 (where `n` is zero), to ensure that two elements of the matrix
686 space (empty matrices) can be multiplied.
688 Matrix algebras override this with something more useful.
690 if self
.is_trivial():
691 return MatrixSpace(self
.base_ring(), 0)
692 elif self
._matrix
_basis
is None or len(self
._matrix
_basis
) == 0:
693 return MatrixSpace(self
.base_ring(), self
.dimension(), 1)
695 return self
._matrix
_basis
[0].matrix_space()
701 Return the unit element of this algebra.
705 sage: from mjo.eja.eja_algebra import (HadamardEJA,
710 sage: J = HadamardEJA(5)
712 e0 + e1 + e2 + e3 + e4
716 The identity element acts like the identity::
718 sage: set_random_seed()
719 sage: J = random_eja()
720 sage: x = J.random_element()
721 sage: J.one()*x == x and x*J.one() == x
724 The matrix of the unit element's operator is the identity::
726 sage: set_random_seed()
727 sage: J = random_eja()
728 sage: actual = J.one().operator().matrix()
729 sage: expected = matrix.identity(J.base_ring(), J.dimension())
730 sage: actual == expected
733 Ensure that the cached unit element (often precomputed by
734 hand) agrees with the computed one::
736 sage: set_random_seed()
737 sage: J = random_eja()
738 sage: cached = J.one()
739 sage: J.one.clear_cache()
740 sage: J.one() == cached
744 # We can brute-force compute the matrices of the operators
745 # that correspond to the basis elements of this algebra.
746 # If some linear combination of those basis elements is the
747 # algebra identity, then the same linear combination of
748 # their matrices has to be the identity matrix.
750 # Of course, matrices aren't vectors in sage, so we have to
751 # appeal to the "long vectors" isometry.
752 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
754 # Now we use basic linear algebra to find the coefficients,
755 # of the matrices-as-vectors-linear-combination, which should
756 # work for the original algebra basis too.
757 A
= matrix(self
.base_ring(), oper_vecs
)
759 # We used the isometry on the left-hand side already, but we
760 # still need to do it for the right-hand side. Recall that we
761 # wanted something that summed to the identity matrix.
762 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
764 # Now if there's an identity element in the algebra, this
765 # should work. We solve on the left to avoid having to
766 # transpose the matrix "A".
767 return self
.from_vector(A
.solve_left(b
))
770 def peirce_decomposition(self
, c
):
772 The Peirce decomposition of this algebra relative to the
775 In the future, this can be extended to a complete system of
776 orthogonal idempotents.
780 - ``c`` -- an idempotent of this algebra.
784 A triple (J0, J5, J1) containing two subalgebras and one subspace
787 - ``J0`` -- the algebra on the eigenspace of ``c.operator()``
788 corresponding to the eigenvalue zero.
790 - ``J5`` -- the eigenspace (NOT a subalgebra) of ``c.operator()``
791 corresponding to the eigenvalue one-half.
793 - ``J1`` -- the algebra on the eigenspace of ``c.operator()``
794 corresponding to the eigenvalue one.
796 These are the only possible eigenspaces for that operator, and this
797 algebra is a direct sum of them. The spaces ``J0`` and ``J1`` are
798 orthogonal, and are subalgebras of this algebra with the appropriate
803 sage: from mjo.eja.eja_algebra import random_eja, RealSymmetricEJA
807 The canonical example comes from the symmetric matrices, which
808 decompose into diagonal and off-diagonal parts::
810 sage: J = RealSymmetricEJA(3)
811 sage: C = matrix(QQ, [ [1,0,0],
815 sage: J0,J5,J1 = J.peirce_decomposition(c)
817 Euclidean Jordan algebra of dimension 1...
819 Vector space of degree 6 and dimension 2...
821 Euclidean Jordan algebra of dimension 3...
822 sage: J0.one().to_matrix()
826 sage: orig_df = AA.options.display_format
827 sage: AA.options.display_format = 'radical'
828 sage: J.from_vector(J5.basis()[0]).to_matrix()
832 sage: J.from_vector(J5.basis()[1]).to_matrix()
836 sage: AA.options.display_format = orig_df
837 sage: J1.one().to_matrix()
844 Every algebra decomposes trivially with respect to its identity
847 sage: set_random_seed()
848 sage: J = random_eja()
849 sage: J0,J5,J1 = J.peirce_decomposition(J.one())
850 sage: J0.dimension() == 0 and J5.dimension() == 0
852 sage: J1.superalgebra() == J and J1.dimension() == J.dimension()
855 The decomposition is into eigenspaces, and its components are
856 therefore necessarily orthogonal. Moreover, the identity
857 elements in the two subalgebras are the projections onto their
858 respective subspaces of the superalgebra's identity element::
860 sage: set_random_seed()
861 sage: J = random_eja()
862 sage: x = J.random_element()
863 sage: if not J.is_trivial():
864 ....: while x.is_nilpotent():
865 ....: x = J.random_element()
866 sage: c = x.subalgebra_idempotent()
867 sage: J0,J5,J1 = J.peirce_decomposition(c)
869 sage: for (w,y,z) in zip(J0.basis(), J5.basis(), J1.basis()):
870 ....: w = w.superalgebra_element()
871 ....: y = J.from_vector(y)
872 ....: z = z.superalgebra_element()
873 ....: ipsum += w.inner_product(y).abs()
874 ....: ipsum += w.inner_product(z).abs()
875 ....: ipsum += y.inner_product(z).abs()
878 sage: J1(c) == J1.one()
880 sage: J0(J.one() - c) == J0.one()
884 if not c
.is_idempotent():
885 raise ValueError("element is not idempotent: %s" % c
)
887 from mjo
.eja
.eja_subalgebra
import FiniteDimensionalEuclideanJordanSubalgebra
889 # Default these to what they should be if they turn out to be
890 # trivial, because eigenspaces_left() won't return eigenvalues
891 # corresponding to trivial spaces (e.g. it returns only the
892 # eigenspace corresponding to lambda=1 if you take the
893 # decomposition relative to the identity element).
894 trivial
= FiniteDimensionalEuclideanJordanSubalgebra(self
, ())
895 J0
= trivial
# eigenvalue zero
896 J5
= VectorSpace(self
.base_ring(), 0) # eigenvalue one-half
897 J1
= trivial
# eigenvalue one
899 for (eigval
, eigspace
) in c
.operator().matrix().right_eigenspaces():
900 if eigval
== ~
(self
.base_ring()(2)):
903 gens
= tuple( self
.from_vector(b
) for b
in eigspace
.basis() )
904 subalg
= FiniteDimensionalEuclideanJordanSubalgebra(self
,
912 raise ValueError("unexpected eigenvalue: %s" % eigval
)
917 def random_element(self
, thorough
=False):
919 Return a random element of this algebra.
921 Our algebra superclass method only returns a linear
922 combination of at most two basis elements. We instead
923 want the vector space "random element" method that
924 returns a more diverse selection.
928 - ``thorough`` -- (boolean; default False) whether or not we
929 should generate irrational coefficients for the random
930 element when our base ring is irrational; this slows the
931 algebra operations to a crawl, but any truly random method
935 # For a general base ring... maybe we can trust this to do the
936 # right thing? Unlikely, but.
937 V
= self
.vector_space()
938 v
= V
.random_element()
940 if self
.base_ring() is AA
:
941 # The "random element" method of the algebraic reals is
942 # stupid at the moment, and only returns integers between
943 # -2 and 2, inclusive:
945 # https://trac.sagemath.org/ticket/30875
947 # Instead, we implement our own "random vector" method,
948 # and then coerce that into the algebra. We use the vector
949 # space degree here instead of the dimension because a
950 # subalgebra could (for example) be spanned by only two
951 # vectors, each with five coordinates. We need to
952 # generate all five coordinates.
954 v
*= QQbar
.random_element().real()
956 v
*= QQ
.random_element()
958 return self
.from_vector(V
.coordinate_vector(v
))
960 def random_elements(self
, count
, thorough
=False):
962 Return ``count`` random elements as a tuple.
966 - ``thorough`` -- (boolean; default False) whether or not we
967 should generate irrational coefficients for the random
968 elements when our base ring is irrational; this slows the
969 algebra operations to a crawl, but any truly random method
974 sage: from mjo.eja.eja_algebra import JordanSpinEJA
978 sage: J = JordanSpinEJA(3)
979 sage: x,y,z = J.random_elements(3)
980 sage: all( [ x in J, y in J, z in J ])
982 sage: len( J.random_elements(10) ) == 10
986 return tuple( self
.random_element(thorough
)
987 for idx
in range(count
) )
991 def _charpoly_coefficients(self
):
993 The `r` polynomial coefficients of the "characteristic polynomial
997 R
= self
.coordinate_polynomial_ring()
999 F
= R
.fraction_field()
1002 # From a result in my book, these are the entries of the
1003 # basis representation of L_x.
1004 return sum( vars[k
]*self
.monomial(k
).operator().matrix()[i
,j
]
1007 L_x
= matrix(F
, n
, n
, L_x_i_j
)
1010 if self
.rank
.is_in_cache():
1012 # There's no need to pad the system with redundant
1013 # columns if we *know* they'll be redundant.
1016 # Compute an extra power in case the rank is equal to
1017 # the dimension (otherwise, we would stop at x^(r-1)).
1018 x_powers
= [ (L_x
**k
)*self
.one().to_vector()
1019 for k
in range(n
+1) ]
1020 A
= matrix
.column(F
, x_powers
[:n
])
1021 AE
= A
.extended_echelon_form()
1028 # The theory says that only the first "r" coefficients are
1029 # nonzero, and they actually live in the original polynomial
1030 # ring and not the fraction field. We negate them because
1031 # in the actual characteristic polynomial, they get moved
1032 # to the other side where x^r lives.
1033 return -A_rref
.solve_right(E
*b
).change_ring(R
)[:r
]
1038 Return the rank of this EJA.
1040 This is a cached method because we know the rank a priori for
1041 all of the algebras we can construct. Thus we can avoid the
1042 expensive ``_charpoly_coefficients()`` call unless we truly
1043 need to compute the whole characteristic polynomial.
1047 sage: from mjo.eja.eja_algebra import (HadamardEJA,
1048 ....: JordanSpinEJA,
1049 ....: RealSymmetricEJA,
1050 ....: ComplexHermitianEJA,
1051 ....: QuaternionHermitianEJA,
1056 The rank of the Jordan spin algebra is always two::
1058 sage: JordanSpinEJA(2).rank()
1060 sage: JordanSpinEJA(3).rank()
1062 sage: JordanSpinEJA(4).rank()
1065 The rank of the `n`-by-`n` Hermitian real, complex, or
1066 quaternion matrices is `n`::
1068 sage: RealSymmetricEJA(4).rank()
1070 sage: ComplexHermitianEJA(3).rank()
1072 sage: QuaternionHermitianEJA(2).rank()
1077 Ensure that every EJA that we know how to construct has a
1078 positive integer rank, unless the algebra is trivial in
1079 which case its rank will be zero::
1081 sage: set_random_seed()
1082 sage: J = random_eja()
1086 sage: r > 0 or (r == 0 and J.is_trivial())
1089 Ensure that computing the rank actually works, since the ranks
1090 of all simple algebras are known and will be cached by default::
1092 sage: set_random_seed() # long time
1093 sage: J = random_eja() # long time
1094 sage: caches = J.rank() # long time
1095 sage: J.rank.clear_cache() # long time
1096 sage: J.rank() == cached # long time
1100 return len(self
._charpoly
_coefficients
())
1103 def vector_space(self
):
1105 Return the vector space that underlies this algebra.
1109 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1113 sage: J = RealSymmetricEJA(2)
1114 sage: J.vector_space()
1115 Vector space of dimension 3 over...
1118 return self
.zero().to_vector().parent().ambient_vector_space()
1121 Element
= FiniteDimensionalEuclideanJordanAlgebraElement
1123 class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra
):
1125 New class for algebras whose supplied basis elements have all rational entries.
1129 sage: from mjo.eja.eja_algebra import BilinearFormEJA
1133 The supplied basis is orthonormalized by default::
1135 sage: B = matrix(QQ, [[1, 0, 0], [0, 25, -32], [0, -32, 41]])
1136 sage: J = BilinearFormEJA(B)
1137 sage: J.matrix_basis()
1150 orthonormalize
=True,
1157 vector_basis
= basis
1159 from sage
.structure
.element
import is_Matrix
1160 basis_is_matrices
= False
1164 if is_Matrix(basis
[0]):
1165 basis_is_matrices
= True
1166 from mjo
.eja
.eja_utils
import _vec2mat
1167 vector_basis
= tuple( map(_mat2vec
,basis
) )
1168 degree
= basis
[0].nrows()**2
1170 degree
= basis
[0].degree()
1172 V
= VectorSpace(field
, degree
)
1174 # If we were asked to orthonormalize, and if the orthonormal
1175 # basis is different from the given one, then we also want to
1176 # compute multiplication and inner-product tables for the
1177 # deorthonormalized basis. These can be used later to
1178 # construct a deorthonormalized copy of this algebra over QQ
1179 # in which several operations are much faster.
1180 self
._deortho
_multiplication
_table
= None
1181 self
._deortho
_inner
_product
_table
= None
1184 # Compute the deorthonormalized tables before we orthonormalize
1186 W
= V
.span_of_basis( vector_basis
)
1188 # TODO: use symmetry
1189 self
._deortho
_multiplication
_table
= [ [0 for j
in range(n
)]
1191 self
._deortho
_inner
_product
_table
= [ [0 for j
in range(n
)]
1194 # Note: the Jordan and inner-products are defined in terms
1195 # of the ambient basis. It's important that their arguments
1196 # are in ambient coordinates as well.
1198 for j
in range(i
+1):
1199 # given basis w.r.t. ambient coords
1200 q_i
= vector_basis
[i
]
1201 q_j
= vector_basis
[j
]
1203 if basis_is_matrices
:
1207 elt
= jordan_product(q_i
, q_j
)
1208 ip
= inner_product(q_i
, q_j
)
1210 if basis_is_matrices
:
1211 # do another mat2vec because the multiplication
1212 # table is in terms of vectors
1215 # TODO: use symmetry
1216 elt
= W
.coordinate_vector(elt
)
1217 self
._deortho
_multiplication
_table
[i
][j
] = elt
1218 self
._deortho
_multiplication
_table
[j
][i
] = elt
1219 self
._deortho
_inner
_product
_table
[i
][j
] = ip
1220 self
._deortho
_inner
_product
_table
[j
][i
] = ip
1222 if self
._deortho
_multiplication
_table
is not None:
1223 self
._deortho
_multiplication
_table
= tuple(map(tuple, self
._deortho
_multiplication
_table
))
1224 if self
._deortho
_inner
_product
_table
is not None:
1225 self
._deortho
_inner
_product
_table
= tuple(map(tuple, self
._deortho
_inner
_product
_table
))
1227 # We overwrite the name "vector_basis" in a second, but never modify it
1228 # in place, to this effectively makes a copy of it.
1229 deortho_vector_basis
= vector_basis
1230 self
._deortho
_matrix
= None
1233 from mjo
.eja
.eja_utils
import gram_schmidt
1234 if basis_is_matrices
:
1235 vector_ip
= lambda x
,y
: inner_product(_vec2mat(x
), _vec2mat(y
))
1236 vector_basis
= gram_schmidt(vector_basis
, vector_ip
)
1238 vector_basis
= gram_schmidt(vector_basis
, inner_product
)
1240 W
= V
.span_of_basis( vector_basis
)
1242 # Normalize the "matrix" basis, too!
1243 basis
= vector_basis
1245 if basis_is_matrices
:
1246 basis
= tuple( map(_vec2mat
,basis
) )
1248 W
= V
.span_of_basis( vector_basis
)
1250 # Now "W" is the vector space of our algebra coordinates. The
1251 # variables "X1", "X2",... refer to the entries of vectors in
1252 # W. Thus to convert back and forth between the orthonormal
1253 # coordinates and the given ones, we need to stick the original
1255 U
= V
.span_of_basis( deortho_vector_basis
)
1256 self
._deortho
_matrix
= matrix( U
.coordinate_vector(q
)
1257 for q
in vector_basis
)
1259 # TODO: use symmetry
1260 mult_table
= [ [0 for j
in range(n
)] for i
in range(n
) ]
1261 ip_table
= [ [0 for j
in range(n
)] for i
in range(n
) ]
1263 # Note: the Jordan and inner-products are defined in terms
1264 # of the ambient basis. It's important that their arguments
1265 # are in ambient coordinates as well.
1267 for j
in range(i
+1):
1268 # ortho basis w.r.t. ambient coords
1269 q_i
= vector_basis
[i
]
1270 q_j
= vector_basis
[j
]
1272 if basis_is_matrices
:
1276 elt
= jordan_product(q_i
, q_j
)
1277 ip
= inner_product(q_i
, q_j
)
1279 if basis_is_matrices
:
1280 # do another mat2vec because the multiplication
1281 # table is in terms of vectors
1284 # TODO: use symmetry
1285 elt
= W
.coordinate_vector(elt
)
1286 mult_table
[i
][j
] = elt
1287 mult_table
[j
][i
] = elt
1291 if basis_is_matrices
:
1295 basis
= tuple( x
.column() for x
in basis
)
1297 super().__init
__(field
,
1302 basis
, # matrix basis
1307 def _charpoly_coefficients(self
):
1311 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
1312 ....: JordanSpinEJA)
1316 The base ring of the resulting polynomial coefficients is what
1317 it should be, and not the rationals (unless the algebra was
1318 already over the rationals)::
1320 sage: J = JordanSpinEJA(3)
1321 sage: J._charpoly_coefficients()
1322 (X1^2 - X2^2 - X3^2, -2*X1)
1323 sage: a0 = J._charpoly_coefficients()[0]
1325 Algebraic Real Field
1326 sage: a0.base_ring()
1327 Algebraic Real Field
1330 if self
.base_ring() is QQ
:
1331 # There's no need to construct *another* algebra over the
1332 # rationals if this one is already over the rationals.
1333 superclass
= super(RationalBasisEuclideanJordanAlgebra
, self
)
1334 return superclass
._charpoly
_coefficients
()
1336 # Do the computation over the rationals. The answer will be
1337 # the same, because all we've done is a change of basis.
1338 J
= FiniteDimensionalEuclideanJordanAlgebra(QQ
,
1339 self
._deortho
_multiplication
_table
,
1340 self
._deortho
_inner
_product
_table
)
1342 # Change back from QQ to our real base ring
1343 a
= ( a_i
.change_ring(self
.base_ring())
1344 for a_i
in J
._charpoly
_coefficients
() )
1346 # Now convert the coordinate variables back to the
1347 # deorthonormalized ones.
1348 R
= self
.coordinate_polynomial_ring()
1349 from sage
.modules
.free_module_element
import vector
1350 X
= vector(R
, R
.gens())
1351 BX
= self
._deortho
_matrix
*X
1353 subs_dict
= { X[i]: BX[i] for i in range(len(X)) }
1354 return tuple( a_i
.subs(subs_dict
) for a_i
in a
)
1356 class ConcreteEuclideanJordanAlgebra(RationalBasisEuclideanJordanAlgebra
):
1358 A class for the Euclidean Jordan algebras that we know by name.
1360 These are the Jordan algebras whose basis, multiplication table,
1361 rank, and so on are known a priori. More to the point, they are
1362 the Euclidean Jordan algebras for which we are able to conjure up
1363 a "random instance."
1367 sage: from mjo.eja.eja_algebra import ConcreteEuclideanJordanAlgebra
1371 Our basis is normalized with respect to the algebra's inner
1372 product, unless we specify otherwise::
1374 sage: set_random_seed()
1375 sage: J = ConcreteEuclideanJordanAlgebra.random_instance()
1376 sage: all( b.norm() == 1 for b in J.gens() )
1379 Since our basis is orthonormal with respect to the algebra's inner
1380 product, and since we know that this algebra is an EJA, any
1381 left-multiplication operator's matrix will be symmetric because
1382 natural->EJA basis representation is an isometry and within the
1383 EJA the operator is self-adjoint by the Jordan axiom::
1385 sage: set_random_seed()
1386 sage: J = ConcreteEuclideanJordanAlgebra.random_instance()
1387 sage: x = J.random_element()
1388 sage: x.operator().is_self_adjoint()
1393 def _max_random_instance_size():
1395 Return an integer "size" that is an upper bound on the size of
1396 this algebra when it is used in a random test
1397 case. Unfortunately, the term "size" is ambiguous -- when
1398 dealing with `R^n` under either the Hadamard or Jordan spin
1399 product, the "size" refers to the dimension `n`. When dealing
1400 with a matrix algebra (real symmetric or complex/quaternion
1401 Hermitian), it refers to the size of the matrix, which is far
1402 less than the dimension of the underlying vector space.
1404 This method must be implemented in each subclass.
1406 raise NotImplementedError
1409 def random_instance(cls
, *args
, **kwargs
):
1411 Return a random instance of this type of algebra.
1413 This method should be implemented in each subclass.
1415 from sage
.misc
.prandom
import choice
1416 eja_class
= choice(cls
.__subclasses
__())
1418 # These all bubble up to the RationalBasisEuclideanJordanAlgebra
1419 # superclass constructor, so any (kw)args valid there are also
1421 return eja_class
.random_instance(*args
, **kwargs
)
1424 class MatrixEuclideanJordanAlgebra
:
1428 Embed the matrix ``M`` into a space of real matrices.
1430 The matrix ``M`` can have entries in any field at the moment:
1431 the real numbers, complex numbers, or quaternions. And although
1432 they are not a field, we can probably support octonions at some
1433 point, too. This function returns a real matrix that "acts like"
1434 the original with respect to matrix multiplication; i.e.
1436 real_embed(M*N) = real_embed(M)*real_embed(N)
1439 raise NotImplementedError
1443 def real_unembed(M
):
1445 The inverse of :meth:`real_embed`.
1447 raise NotImplementedError
1450 def jordan_product(X
,Y
):
1451 return (X
*Y
+ Y
*X
)/2
1454 def trace_inner_product(cls
,X
,Y
):
1455 Xu
= cls
.real_unembed(X
)
1456 Yu
= cls
.real_unembed(Y
)
1457 tr
= (Xu
*Yu
).trace()
1460 # Works in QQ, AA, RDF, et cetera.
1462 except AttributeError:
1463 # A quaternion doesn't have a real() method, but does
1464 # have coefficient_tuple() method that returns the
1465 # coefficients of 1, i, j, and k -- in that order.
1466 return tr
.coefficient_tuple()[0]
1469 class RealMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1473 The identity function, for embedding real matrices into real
1479 def real_unembed(M
):
1481 The identity function, for unembedding real matrices from real
1487 class RealSymmetricEJA(ConcreteEuclideanJordanAlgebra
,
1488 RealMatrixEuclideanJordanAlgebra
):
1490 The rank-n simple EJA consisting of real symmetric n-by-n
1491 matrices, the usual symmetric Jordan product, and the trace inner
1492 product. It has dimension `(n^2 + n)/2` over the reals.
1496 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1500 sage: J = RealSymmetricEJA(2)
1501 sage: e0, e1, e2 = J.gens()
1509 In theory, our "field" can be any subfield of the reals::
1511 sage: RealSymmetricEJA(2, RDF)
1512 Euclidean Jordan algebra of dimension 3 over Real Double Field
1513 sage: RealSymmetricEJA(2, RR)
1514 Euclidean Jordan algebra of dimension 3 over Real Field with
1515 53 bits of precision
1519 The dimension of this algebra is `(n^2 + n) / 2`::
1521 sage: set_random_seed()
1522 sage: n_max = RealSymmetricEJA._max_random_instance_size()
1523 sage: n = ZZ.random_element(1, n_max)
1524 sage: J = RealSymmetricEJA(n)
1525 sage: J.dimension() == (n^2 + n)/2
1528 The Jordan multiplication is what we think it is::
1530 sage: set_random_seed()
1531 sage: J = RealSymmetricEJA.random_instance()
1532 sage: x,y = J.random_elements(2)
1533 sage: actual = (x*y).to_matrix()
1534 sage: X = x.to_matrix()
1535 sage: Y = y.to_matrix()
1536 sage: expected = (X*Y + Y*X)/2
1537 sage: actual == expected
1539 sage: J(expected) == x*y
1542 We can change the generator prefix::
1544 sage: RealSymmetricEJA(3, prefix='q').gens()
1545 (q0, q1, q2, q3, q4, q5)
1547 We can construct the (trivial) algebra of rank zero::
1549 sage: RealSymmetricEJA(0)
1550 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1554 def _denormalized_basis(cls
, n
, field
):
1556 Return a basis for the space of real symmetric n-by-n matrices.
1560 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1564 sage: set_random_seed()
1565 sage: n = ZZ.random_element(1,5)
1566 sage: B = RealSymmetricEJA._denormalized_basis(n,QQ)
1567 sage: all( M.is_symmetric() for M in B)
1571 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1575 for j
in range(i
+1):
1576 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
1580 Sij
= Eij
+ Eij
.transpose()
1586 def _max_random_instance_size():
1587 return 4 # Dimension 10
1590 def random_instance(cls
, field
=AA
, **kwargs
):
1592 Return a random instance of this type of algebra.
1594 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1595 return cls(n
, field
, **kwargs
)
1597 def __init__(self
, n
, field
=AA
, **kwargs
):
1598 basis
= self
._denormalized
_basis
(n
, field
)
1599 super(RealSymmetricEJA
, self
).__init
__(field
,
1601 self
.jordan_product
,
1602 self
.trace_inner_product
,
1604 self
.rank
.set_cache(n
)
1605 self
.one
.set_cache(self(matrix
.identity(field
,n
)))
1608 class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1612 Embed the n-by-n complex matrix ``M`` into the space of real
1613 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1614 bi` to the block matrix ``[[a,b],[-b,a]]``.
1618 sage: from mjo.eja.eja_algebra import \
1619 ....: ComplexMatrixEuclideanJordanAlgebra
1623 sage: F = QuadraticField(-1, 'I')
1624 sage: x1 = F(4 - 2*i)
1625 sage: x2 = F(1 + 2*i)
1628 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1629 sage: ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1638 Embedding is a homomorphism (isomorphism, in fact)::
1640 sage: set_random_seed()
1641 sage: n = ZZ.random_element(3)
1642 sage: F = QuadraticField(-1, 'I')
1643 sage: X = random_matrix(F, n)
1644 sage: Y = random_matrix(F, n)
1645 sage: Xe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X)
1646 sage: Ye = ComplexMatrixEuclideanJordanAlgebra.real_embed(Y)
1647 sage: XYe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1654 raise ValueError("the matrix 'M' must be square")
1656 # We don't need any adjoined elements...
1657 field
= M
.base_ring().base_ring()
1661 a
= z
.list()[0] # real part, I guess
1662 b
= z
.list()[1] # imag part, I guess
1663 blocks
.append(matrix(field
, 2, [[a
,b
],[-b
,a
]]))
1665 return matrix
.block(field
, n
, blocks
)
1669 def real_unembed(M
):
1671 The inverse of _embed_complex_matrix().
1675 sage: from mjo.eja.eja_algebra import \
1676 ....: ComplexMatrixEuclideanJordanAlgebra
1680 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1681 ....: [-2, 1, -4, 3],
1682 ....: [ 9, 10, 11, 12],
1683 ....: [-10, 9, -12, 11] ])
1684 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(A)
1686 [ 10*I + 9 12*I + 11]
1690 Unembedding is the inverse of embedding::
1692 sage: set_random_seed()
1693 sage: F = QuadraticField(-1, 'I')
1694 sage: M = random_matrix(F, 3)
1695 sage: Me = ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1696 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1702 raise ValueError("the matrix 'M' must be square")
1703 if not n
.mod(2).is_zero():
1704 raise ValueError("the matrix 'M' must be a complex embedding")
1706 # If "M" was normalized, its base ring might have roots
1707 # adjoined and they can stick around after unembedding.
1708 field
= M
.base_ring()
1709 R
= PolynomialRing(field
, 'z')
1712 # Sage doesn't know how to embed AA into QQbar, i.e. how
1713 # to adjoin sqrt(-1) to AA.
1716 F
= field
.extension(z
**2 + 1, 'I', embedding
=CLF(-1).sqrt())
1719 # Go top-left to bottom-right (reading order), converting every
1720 # 2-by-2 block we see to a single complex element.
1722 for k
in range(n
/2):
1723 for j
in range(n
/2):
1724 submat
= M
[2*k
:2*k
+2,2*j
:2*j
+2]
1725 if submat
[0,0] != submat
[1,1]:
1726 raise ValueError('bad on-diagonal submatrix')
1727 if submat
[0,1] != -submat
[1,0]:
1728 raise ValueError('bad off-diagonal submatrix')
1729 z
= submat
[0,0] + submat
[0,1]*i
1732 return matrix(F
, n
/2, elements
)
1736 def trace_inner_product(cls
,X
,Y
):
1738 Compute a matrix inner product in this algebra directly from
1743 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1747 This gives the same answer as the slow, default method implemented
1748 in :class:`MatrixEuclideanJordanAlgebra`::
1750 sage: set_random_seed()
1751 sage: J = ComplexHermitianEJA.random_instance()
1752 sage: x,y = J.random_elements(2)
1753 sage: Xe = x.to_matrix()
1754 sage: Ye = y.to_matrix()
1755 sage: X = ComplexHermitianEJA.real_unembed(Xe)
1756 sage: Y = ComplexHermitianEJA.real_unembed(Ye)
1757 sage: expected = (X*Y).trace().real()
1758 sage: actual = ComplexHermitianEJA.trace_inner_product(Xe,Ye)
1759 sage: actual == expected
1763 return RealMatrixEuclideanJordanAlgebra
.trace_inner_product(X
,Y
)/2
1766 class ComplexHermitianEJA(ConcreteEuclideanJordanAlgebra
,
1767 ComplexMatrixEuclideanJordanAlgebra
):
1769 The rank-n simple EJA consisting of complex Hermitian n-by-n
1770 matrices over the real numbers, the usual symmetric Jordan product,
1771 and the real-part-of-trace inner product. It has dimension `n^2` over
1776 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1780 In theory, our "field" can be any subfield of the reals::
1782 sage: ComplexHermitianEJA(2, RDF)
1783 Euclidean Jordan algebra of dimension 4 over Real Double Field
1784 sage: ComplexHermitianEJA(2, RR)
1785 Euclidean Jordan algebra of dimension 4 over Real Field with
1786 53 bits of precision
1790 The dimension of this algebra is `n^2`::
1792 sage: set_random_seed()
1793 sage: n_max = ComplexHermitianEJA._max_random_instance_size()
1794 sage: n = ZZ.random_element(1, n_max)
1795 sage: J = ComplexHermitianEJA(n)
1796 sage: J.dimension() == n^2
1799 The Jordan multiplication is what we think it is::
1801 sage: set_random_seed()
1802 sage: J = ComplexHermitianEJA.random_instance()
1803 sage: x,y = J.random_elements(2)
1804 sage: actual = (x*y).to_matrix()
1805 sage: X = x.to_matrix()
1806 sage: Y = y.to_matrix()
1807 sage: expected = (X*Y + Y*X)/2
1808 sage: actual == expected
1810 sage: J(expected) == x*y
1813 We can change the generator prefix::
1815 sage: ComplexHermitianEJA(2, prefix='z').gens()
1818 We can construct the (trivial) algebra of rank zero::
1820 sage: ComplexHermitianEJA(0)
1821 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1826 def _denormalized_basis(cls
, n
, field
):
1828 Returns a basis for the space of complex Hermitian n-by-n matrices.
1830 Why do we embed these? Basically, because all of numerical linear
1831 algebra assumes that you're working with vectors consisting of `n`
1832 entries from a field and scalars from the same field. There's no way
1833 to tell SageMath that (for example) the vectors contain complex
1834 numbers, while the scalar field is real.
1838 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1842 sage: set_random_seed()
1843 sage: n = ZZ.random_element(1,5)
1844 sage: field = QuadraticField(2, 'sqrt2')
1845 sage: B = ComplexHermitianEJA._denormalized_basis(n, field)
1846 sage: all( M.is_symmetric() for M in B)
1850 R
= PolynomialRing(field
, 'z')
1852 F
= field
.extension(z
**2 + 1, 'I')
1855 # This is like the symmetric case, but we need to be careful:
1857 # * We want conjugate-symmetry, not just symmetry.
1858 # * The diagonal will (as a result) be real.
1862 for j
in range(i
+1):
1863 Eij
= matrix(F
, n
, lambda k
,l
: k
==i
and l
==j
)
1865 Sij
= cls
.real_embed(Eij
)
1868 # The second one has a minus because it's conjugated.
1869 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
1871 Sij_imag
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
1874 # Since we embedded these, we can drop back to the "field" that we
1875 # started with instead of the complex extension "F".
1876 return tuple( s
.change_ring(field
) for s
in S
)
1879 def __init__(self
, n
, field
=AA
, **kwargs
):
1880 basis
= self
._denormalized
_basis
(n
,field
)
1881 super(ComplexHermitianEJA
, self
).__init
__(field
,
1883 self
.jordan_product
,
1884 self
.trace_inner_product
,
1886 self
.rank
.set_cache(n
)
1887 # TODO: pre-cache the identity!
1890 def _max_random_instance_size():
1891 return 3 # Dimension 9
1894 def random_instance(cls
, field
=AA
, **kwargs
):
1896 Return a random instance of this type of algebra.
1898 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1899 return cls(n
, field
, **kwargs
)
1901 class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1905 Embed the n-by-n quaternion matrix ``M`` into the space of real
1906 matrices of size 4n-by-4n by first sending each quaternion entry `z
1907 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
1908 c+di],[-c + di, a-bi]]`, and then embedding those into a real
1913 sage: from mjo.eja.eja_algebra import \
1914 ....: QuaternionMatrixEuclideanJordanAlgebra
1918 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1919 sage: i,j,k = Q.gens()
1920 sage: x = 1 + 2*i + 3*j + 4*k
1921 sage: M = matrix(Q, 1, [[x]])
1922 sage: QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1928 Embedding is a homomorphism (isomorphism, in fact)::
1930 sage: set_random_seed()
1931 sage: n = ZZ.random_element(2)
1932 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1933 sage: X = random_matrix(Q, n)
1934 sage: Y = random_matrix(Q, n)
1935 sage: Xe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X)
1936 sage: Ye = QuaternionMatrixEuclideanJordanAlgebra.real_embed(Y)
1937 sage: XYe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1942 quaternions
= M
.base_ring()
1945 raise ValueError("the matrix 'M' must be square")
1947 F
= QuadraticField(-1, 'I')
1952 t
= z
.coefficient_tuple()
1957 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
1958 [-c
+ d
*i
, a
- b
*i
]])
1959 realM
= ComplexMatrixEuclideanJordanAlgebra
.real_embed(cplxM
)
1960 blocks
.append(realM
)
1962 # We should have real entries by now, so use the realest field
1963 # we've got for the return value.
1964 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
1969 def real_unembed(M
):
1971 The inverse of _embed_quaternion_matrix().
1975 sage: from mjo.eja.eja_algebra import \
1976 ....: QuaternionMatrixEuclideanJordanAlgebra
1980 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
1981 ....: [-2, 1, -4, 3],
1982 ....: [-3, 4, 1, -2],
1983 ....: [-4, -3, 2, 1]])
1984 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(M)
1985 [1 + 2*i + 3*j + 4*k]
1989 Unembedding is the inverse of embedding::
1991 sage: set_random_seed()
1992 sage: Q = QuaternionAlgebra(QQ, -1, -1)
1993 sage: M = random_matrix(Q, 3)
1994 sage: Me = QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1995 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
2001 raise ValueError("the matrix 'M' must be square")
2002 if not n
.mod(4).is_zero():
2003 raise ValueError("the matrix 'M' must be a quaternion embedding")
2005 # Use the base ring of the matrix to ensure that its entries can be
2006 # multiplied by elements of the quaternion algebra.
2007 field
= M
.base_ring()
2008 Q
= QuaternionAlgebra(field
,-1,-1)
2011 # Go top-left to bottom-right (reading order), converting every
2012 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
2015 for l
in range(n
/4):
2016 for m
in range(n
/4):
2017 submat
= ComplexMatrixEuclideanJordanAlgebra
.real_unembed(
2018 M
[4*l
:4*l
+4,4*m
:4*m
+4] )
2019 if submat
[0,0] != submat
[1,1].conjugate():
2020 raise ValueError('bad on-diagonal submatrix')
2021 if submat
[0,1] != -submat
[1,0].conjugate():
2022 raise ValueError('bad off-diagonal submatrix')
2023 z
= submat
[0,0].real()
2024 z
+= submat
[0,0].imag()*i
2025 z
+= submat
[0,1].real()*j
2026 z
+= submat
[0,1].imag()*k
2029 return matrix(Q
, n
/4, elements
)
2033 def trace_inner_product(cls
,X
,Y
):
2035 Compute a matrix inner product in this algebra directly from
2040 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2044 This gives the same answer as the slow, default method implemented
2045 in :class:`MatrixEuclideanJordanAlgebra`::
2047 sage: set_random_seed()
2048 sage: J = QuaternionHermitianEJA.random_instance()
2049 sage: x,y = J.random_elements(2)
2050 sage: Xe = x.to_matrix()
2051 sage: Ye = y.to_matrix()
2052 sage: X = QuaternionHermitianEJA.real_unembed(Xe)
2053 sage: Y = QuaternionHermitianEJA.real_unembed(Ye)
2054 sage: expected = (X*Y).trace().coefficient_tuple()[0]
2055 sage: actual = QuaternionHermitianEJA.trace_inner_product(Xe,Ye)
2056 sage: actual == expected
2060 return RealMatrixEuclideanJordanAlgebra
.trace_inner_product(X
,Y
)/4
2063 class QuaternionHermitianEJA(ConcreteEuclideanJordanAlgebra
,
2064 QuaternionMatrixEuclideanJordanAlgebra
):
2066 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
2067 matrices, the usual symmetric Jordan product, and the
2068 real-part-of-trace inner product. It has dimension `2n^2 - n` over
2073 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2077 In theory, our "field" can be any subfield of the reals::
2079 sage: QuaternionHermitianEJA(2, RDF)
2080 Euclidean Jordan algebra of dimension 6 over Real Double Field
2081 sage: QuaternionHermitianEJA(2, RR)
2082 Euclidean Jordan algebra of dimension 6 over Real Field with
2083 53 bits of precision
2087 The dimension of this algebra is `2*n^2 - n`::
2089 sage: set_random_seed()
2090 sage: n_max = QuaternionHermitianEJA._max_random_instance_size()
2091 sage: n = ZZ.random_element(1, n_max)
2092 sage: J = QuaternionHermitianEJA(n)
2093 sage: J.dimension() == 2*(n^2) - n
2096 The Jordan multiplication is what we think it is::
2098 sage: set_random_seed()
2099 sage: J = QuaternionHermitianEJA.random_instance()
2100 sage: x,y = J.random_elements(2)
2101 sage: actual = (x*y).to_matrix()
2102 sage: X = x.to_matrix()
2103 sage: Y = y.to_matrix()
2104 sage: expected = (X*Y + Y*X)/2
2105 sage: actual == expected
2107 sage: J(expected) == x*y
2110 We can change the generator prefix::
2112 sage: QuaternionHermitianEJA(2, prefix='a').gens()
2113 (a0, a1, a2, a3, a4, a5)
2115 We can construct the (trivial) algebra of rank zero::
2117 sage: QuaternionHermitianEJA(0)
2118 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2122 def _denormalized_basis(cls
, n
, field
):
2124 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
2126 Why do we embed these? Basically, because all of numerical
2127 linear algebra assumes that you're working with vectors consisting
2128 of `n` entries from a field and scalars from the same field. There's
2129 no way to tell SageMath that (for example) the vectors contain
2130 complex numbers, while the scalar field is real.
2134 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2138 sage: set_random_seed()
2139 sage: n = ZZ.random_element(1,5)
2140 sage: B = QuaternionHermitianEJA._denormalized_basis(n,QQ)
2141 sage: all( M.is_symmetric() for M in B )
2145 Q
= QuaternionAlgebra(QQ
,-1,-1)
2148 # This is like the symmetric case, but we need to be careful:
2150 # * We want conjugate-symmetry, not just symmetry.
2151 # * The diagonal will (as a result) be real.
2155 for j
in range(i
+1):
2156 Eij
= matrix(Q
, n
, lambda k
,l
: k
==i
and l
==j
)
2158 Sij
= cls
.real_embed(Eij
)
2161 # The second, third, and fourth ones have a minus
2162 # because they're conjugated.
2163 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
2165 Sij_I
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
2167 Sij_J
= cls
.real_embed(J
*Eij
- J
*Eij
.transpose())
2169 Sij_K
= cls
.real_embed(K
*Eij
- K
*Eij
.transpose())
2172 # Since we embedded these, we can drop back to the "field" that we
2173 # started with instead of the quaternion algebra "Q".
2174 return tuple( s
.change_ring(field
) for s
in S
)
2177 def __init__(self
, n
, field
=AA
, **kwargs
):
2178 basis
= self
._denormalized
_basis
(n
,field
)
2179 super(QuaternionHermitianEJA
, self
).__init
__(field
,
2181 self
.jordan_product
,
2182 self
.trace_inner_product
,
2184 self
.rank
.set_cache(n
)
2185 # TODO: cache one()!
2188 def _max_random_instance_size():
2190 The maximum rank of a random QuaternionHermitianEJA.
2192 return 2 # Dimension 6
2195 def random_instance(cls
, field
=AA
, **kwargs
):
2197 Return a random instance of this type of algebra.
2199 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2200 return cls(n
, field
, **kwargs
)
2203 class HadamardEJA(ConcreteEuclideanJordanAlgebra
):
2205 Return the Euclidean Jordan Algebra corresponding to the set
2206 `R^n` under the Hadamard product.
2208 Note: this is nothing more than the Cartesian product of ``n``
2209 copies of the spin algebra. Once Cartesian product algebras
2210 are implemented, this can go.
2214 sage: from mjo.eja.eja_algebra import HadamardEJA
2218 This multiplication table can be verified by hand::
2220 sage: J = HadamardEJA(3)
2221 sage: e0,e1,e2 = J.gens()
2237 We can change the generator prefix::
2239 sage: HadamardEJA(3, prefix='r').gens()
2243 def __init__(self
, n
, field
=AA
, **kwargs
):
2244 V
= VectorSpace(field
, n
)
2247 def jordan_product(x
,y
):
2248 return V([ xi
*yi
for (xi
,yi
) in zip(x
,y
) ])
2249 def inner_product(x
,y
):
2250 return x
.inner_product(y
)
2252 super(HadamardEJA
, self
).__init
__(field
,
2257 self
.rank
.set_cache(n
)
2260 self
.one
.set_cache( self
.zero() )
2262 self
.one
.set_cache( sum(self
.gens()) )
2265 def _max_random_instance_size():
2267 The maximum dimension of a random HadamardEJA.
2272 def random_instance(cls
, field
=AA
, **kwargs
):
2274 Return a random instance of this type of algebra.
2276 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2277 return cls(n
, field
, **kwargs
)
2280 class BilinearFormEJA(ConcreteEuclideanJordanAlgebra
):
2282 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2283 with the half-trace inner product and jordan product ``x*y =
2284 (<Bx,y>,y_bar>, x0*y_bar + y0*x_bar)`` where `B = 1 \times B22` is
2285 a symmetric positive-definite "bilinear form" matrix. Its
2286 dimension is the size of `B`, and it has rank two in dimensions
2287 larger than two. It reduces to the ``JordanSpinEJA`` when `B` is
2288 the identity matrix of order ``n``.
2290 We insist that the one-by-one upper-left identity block of `B` be
2291 passed in as well so that we can be passed a matrix of size zero
2292 to construct a trivial algebra.
2296 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
2297 ....: JordanSpinEJA)
2301 When no bilinear form is specified, the identity matrix is used,
2302 and the resulting algebra is the Jordan spin algebra::
2304 sage: B = matrix.identity(AA,3)
2305 sage: J0 = BilinearFormEJA(B)
2306 sage: J1 = JordanSpinEJA(3)
2307 sage: J0.multiplication_table() == J0.multiplication_table()
2310 An error is raised if the matrix `B` does not correspond to a
2311 positive-definite bilinear form::
2313 sage: B = matrix.random(QQ,2,3)
2314 sage: J = BilinearFormEJA(B)
2315 Traceback (most recent call last):
2317 ValueError: bilinear form is not positive-definite
2318 sage: B = matrix.zero(QQ,3)
2319 sage: J = BilinearFormEJA(B)
2320 Traceback (most recent call last):
2322 ValueError: bilinear form is not positive-definite
2326 We can create a zero-dimensional algebra::
2328 sage: B = matrix.identity(AA,0)
2329 sage: J = BilinearFormEJA(B)
2333 We can check the multiplication condition given in the Jordan, von
2334 Neumann, and Wigner paper (and also discussed on my "On the
2335 symmetry..." paper). Note that this relies heavily on the standard
2336 choice of basis, as does anything utilizing the bilinear form
2337 matrix. We opt not to orthonormalize the basis, because if we
2338 did, we would have to normalize the `s_{i}` in a similar manner::
2340 sage: set_random_seed()
2341 sage: n = ZZ.random_element(5)
2342 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2343 sage: B11 = matrix.identity(QQ,1)
2344 sage: B22 = M.transpose()*M
2345 sage: B = block_matrix(2,2,[ [B11,0 ],
2347 sage: J = BilinearFormEJA(B, orthonormalize=False)
2348 sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
2349 sage: V = J.vector_space()
2350 sage: sis = [ J( V([0] + (M.inverse()*ei).list()).column() )
2351 ....: for ei in eis ]
2352 sage: actual = [ sis[i]*sis[j]
2353 ....: for i in range(n-1)
2354 ....: for j in range(n-1) ]
2355 sage: expected = [ J.one() if i == j else J.zero()
2356 ....: for i in range(n-1)
2357 ....: for j in range(n-1) ]
2358 sage: actual == expected
2361 def __init__(self
, B
, field
=AA
, **kwargs
):
2362 if not B
.is_positive_definite():
2363 raise ValueError("bilinear form is not positive-definite")
2366 V
= VectorSpace(field
, n
)
2368 def inner_product(x
,y
):
2369 return (B
*x
).inner_product(y
)
2371 def jordan_product(x
,y
):
2376 z0
= inner_product(x
,y
)
2377 zbar
= y0
*xbar
+ x0
*ybar
2378 return V([z0
] + zbar
.list())
2380 super(BilinearFormEJA
, self
).__init
__(field
,
2386 # The rank of this algebra is two, unless we're in a
2387 # one-dimensional ambient space (because the rank is bounded
2388 # by the ambient dimension).
2389 self
.rank
.set_cache(min(n
,2))
2392 self
.one
.set_cache( self
.zero() )
2394 self
.one
.set_cache( self
.monomial(0) )
2397 def _max_random_instance_size():
2399 The maximum dimension of a random BilinearFormEJA.
2404 def random_instance(cls
, field
=AA
, **kwargs
):
2406 Return a random instance of this algebra.
2408 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2410 B
= matrix
.identity(field
, n
)
2411 return cls(B
, field
, **kwargs
)
2413 B11
= matrix
.identity(field
,1)
2414 M
= matrix
.random(field
, n
-1)
2415 I
= matrix
.identity(field
, n
-1)
2416 alpha
= field
.zero()
2417 while alpha
.is_zero():
2418 alpha
= field
.random_element().abs()
2419 B22
= M
.transpose()*M
+ alpha
*I
2421 from sage
.matrix
.special
import block_matrix
2422 B
= block_matrix(2,2, [ [B11
, ZZ(0) ],
2425 return cls(B
, field
, **kwargs
)
2428 class JordanSpinEJA(BilinearFormEJA
):
2430 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2431 with the usual inner product and jordan product ``x*y =
2432 (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2437 sage: from mjo.eja.eja_algebra import JordanSpinEJA
2441 This multiplication table can be verified by hand::
2443 sage: J = JordanSpinEJA(4)
2444 sage: e0,e1,e2,e3 = J.gens()
2460 We can change the generator prefix::
2462 sage: JordanSpinEJA(2, prefix='B').gens()
2467 Ensure that we have the usual inner product on `R^n`::
2469 sage: set_random_seed()
2470 sage: J = JordanSpinEJA.random_instance()
2471 sage: x,y = J.random_elements(2)
2472 sage: actual = x.inner_product(y)
2473 sage: expected = x.to_vector().inner_product(y.to_vector())
2474 sage: actual == expected
2478 def __init__(self
, n
, field
=AA
, **kwargs
):
2479 # This is a special case of the BilinearFormEJA with the identity
2480 # matrix as its bilinear form.
2481 B
= matrix
.identity(field
, n
)
2482 super(JordanSpinEJA
, self
).__init
__(B
, field
, **kwargs
)
2485 def _max_random_instance_size():
2487 The maximum dimension of a random JordanSpinEJA.
2492 def random_instance(cls
, field
=AA
, **kwargs
):
2494 Return a random instance of this type of algebra.
2496 Needed here to override the implementation for ``BilinearFormEJA``.
2498 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2499 return cls(n
, field
, **kwargs
)
2502 class TrivialEJA(ConcreteEuclideanJordanAlgebra
):
2504 The trivial Euclidean Jordan algebra consisting of only a zero element.
2508 sage: from mjo.eja.eja_algebra import TrivialEJA
2512 sage: J = TrivialEJA()
2519 sage: 7*J.one()*12*J.one()
2521 sage: J.one().inner_product(J.one())
2523 sage: J.one().norm()
2525 sage: J.one().subalgebra_generated_by()
2526 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2531 def __init__(self
, field
=AA
, **kwargs
):
2532 jordan_product
= lambda x
,y
: x
2533 inner_product
= lambda x
,y
: field(0)
2535 super(TrivialEJA
, self
).__init
__(field
,
2540 # The rank is zero using my definition, namely the dimension of the
2541 # largest subalgebra generated by any element.
2542 self
.rank
.set_cache(0)
2543 self
.one
.set_cache( self
.zero() )
2546 def random_instance(cls
, field
=AA
, **kwargs
):
2547 # We don't take a "size" argument so the superclass method is
2548 # inappropriate for us.
2549 return cls(field
, **kwargs
)
2551 class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra
):
2553 The external (orthogonal) direct sum of two other Euclidean Jordan
2554 algebras. Essentially the Cartesian product of its two factors.
2555 Every Euclidean Jordan algebra decomposes into an orthogonal
2556 direct sum of simple Euclidean Jordan algebras, so no generality
2557 is lost by providing only this construction.
2561 sage: from mjo.eja.eja_algebra import (random_eja,
2563 ....: RealSymmetricEJA,
2568 sage: J1 = HadamardEJA(2)
2569 sage: J2 = RealSymmetricEJA(3)
2570 sage: J = DirectSumEJA(J1,J2)
2578 The external direct sum construction is only valid when the two factors
2579 have the same base ring; an error is raised otherwise::
2581 sage: set_random_seed()
2582 sage: J1 = random_eja(AA)
2583 sage: J2 = random_eja(QQ)
2584 sage: J = DirectSumEJA(J1,J2)
2585 Traceback (most recent call last):
2587 ValueError: algebras must share the same base field
2590 def __init__(self
, J1
, J2
, **kwargs
):
2591 if J1
.base_ring() != J2
.base_ring():
2592 raise ValueError("algebras must share the same base field")
2593 field
= J1
.base_ring()
2595 self
._factors
= (J1
, J2
)
2599 V
= VectorSpace(field
, n
)
2600 mult_table
= [ [ V
.zero() for j
in range(n
) ]
2604 p
= (J1
.monomial(i
)*J1
.monomial(j
)).to_vector()
2605 mult_table
[i
][j
] = V(p
.list() + [field
.zero()]*n2
)
2609 p
= (J2
.monomial(i
)*J2
.monomial(j
)).to_vector()
2610 mult_table
[n1
+i
][n1
+j
] = V([field
.zero()]*n1
+ p
.list())
2612 # TODO: build the IP table here from the two constituent IP
2613 # matrices (it'll be block diagonal, I think).
2615 super(DirectSumEJA
, self
).__init
__(field
,
2620 self
.rank
.set_cache(J1
.rank() + J2
.rank())
2625 Return the pair of this algebra's factors.
2629 sage: from mjo.eja.eja_algebra import (HadamardEJA,
2630 ....: JordanSpinEJA,
2635 sage: J1 = HadamardEJA(2,QQ)
2636 sage: J2 = JordanSpinEJA(3,QQ)
2637 sage: J = DirectSumEJA(J1,J2)
2639 (Euclidean Jordan algebra of dimension 2 over Rational Field,
2640 Euclidean Jordan algebra of dimension 3 over Rational Field)
2643 return self
._factors
2645 def projections(self
):
2647 Return a pair of projections onto this algebra's factors.
2651 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
2652 ....: ComplexHermitianEJA,
2657 sage: J1 = JordanSpinEJA(2)
2658 sage: J2 = ComplexHermitianEJA(2)
2659 sage: J = DirectSumEJA(J1,J2)
2660 sage: (pi_left, pi_right) = J.projections()
2661 sage: J.one().to_vector()
2663 sage: pi_left(J.one()).to_vector()
2665 sage: pi_right(J.one()).to_vector()
2669 (J1
,J2
) = self
.factors()
2672 V_basis
= self
.vector_space().basis()
2673 # Need to specify the dimensions explicitly so that we don't
2674 # wind up with a zero-by-zero matrix when we want e.g. a
2675 # zero-by-two matrix (important for composing things).
2676 P1
= matrix(self
.base_ring(), m
, m
+n
, V_basis
[:m
])
2677 P2
= matrix(self
.base_ring(), n
, m
+n
, V_basis
[m
:])
2678 pi_left
= FiniteDimensionalEuclideanJordanAlgebraOperator(self
,J1
,P1
)
2679 pi_right
= FiniteDimensionalEuclideanJordanAlgebraOperator(self
,J2
,P2
)
2680 return (pi_left
, pi_right
)
2682 def inclusions(self
):
2684 Return the pair of inclusion maps from our factors into us.
2688 sage: from mjo.eja.eja_algebra import (random_eja,
2689 ....: JordanSpinEJA,
2690 ....: RealSymmetricEJA,
2695 sage: J1 = JordanSpinEJA(3)
2696 sage: J2 = RealSymmetricEJA(2)
2697 sage: J = DirectSumEJA(J1,J2)
2698 sage: (iota_left, iota_right) = J.inclusions()
2699 sage: iota_left(J1.zero()) == J.zero()
2701 sage: iota_right(J2.zero()) == J.zero()
2703 sage: J1.one().to_vector()
2705 sage: iota_left(J1.one()).to_vector()
2707 sage: J2.one().to_vector()
2709 sage: iota_right(J2.one()).to_vector()
2711 sage: J.one().to_vector()
2716 Composing a projection with the corresponding inclusion should
2717 produce the identity map, and mismatching them should produce
2720 sage: set_random_seed()
2721 sage: J1 = random_eja()
2722 sage: J2 = random_eja()
2723 sage: J = DirectSumEJA(J1,J2)
2724 sage: (iota_left, iota_right) = J.inclusions()
2725 sage: (pi_left, pi_right) = J.projections()
2726 sage: pi_left*iota_left == J1.one().operator()
2728 sage: pi_right*iota_right == J2.one().operator()
2730 sage: (pi_left*iota_right).is_zero()
2732 sage: (pi_right*iota_left).is_zero()
2736 (J1
,J2
) = self
.factors()
2739 V_basis
= self
.vector_space().basis()
2740 # Need to specify the dimensions explicitly so that we don't
2741 # wind up with a zero-by-zero matrix when we want e.g. a
2742 # two-by-zero matrix (important for composing things).
2743 I1
= matrix
.column(self
.base_ring(), m
, m
+n
, V_basis
[:m
])
2744 I2
= matrix
.column(self
.base_ring(), n
, m
+n
, V_basis
[m
:])
2745 iota_left
= FiniteDimensionalEuclideanJordanAlgebraOperator(J1
,self
,I1
)
2746 iota_right
= FiniteDimensionalEuclideanJordanAlgebraOperator(J2
,self
,I2
)
2747 return (iota_left
, iota_right
)
2749 def inner_product(self
, x
, y
):
2751 The standard Cartesian inner-product.
2753 We project ``x`` and ``y`` onto our factors, and add up the
2754 inner-products from the subalgebras.
2759 sage: from mjo.eja.eja_algebra import (HadamardEJA,
2760 ....: QuaternionHermitianEJA,
2765 sage: J1 = HadamardEJA(3,QQ)
2766 sage: J2 = QuaternionHermitianEJA(2,QQ,orthonormalize=False)
2767 sage: J = DirectSumEJA(J1,J2)
2772 sage: x1.inner_product(x2)
2774 sage: y1.inner_product(y2)
2776 sage: J.one().inner_product(J.one())
2780 (pi_left
, pi_right
) = self
.projections()
2786 return (x1
.inner_product(y1
) + x2
.inner_product(y2
))
2790 random_eja
= ConcreteEuclideanJordanAlgebra
.random_instance