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 FiniteDimensionalEJAElement
33 from mjo
.eja
.eja_operator
import FiniteDimensionalEJAOperator
34 from mjo
.eja
.eja_utils
import _mat2vec
36 class FiniteDimensionalEJA(CombinatorialFreeModule
):
38 A finite-dimensional Euclidean Jordan algebra.
42 - basis -- a tuple of basis elements in their matrix form.
44 - jordan_product -- function of two elements (in matrix form)
45 that returns their jordan product in this algebra; this will
46 be applied to ``basis`` to compute a multiplication table for
49 - inner_product -- function of two elements (in matrix form) that
50 returns their inner product. This will be applied to ``basis`` to
51 compute an inner-product table (basically a matrix) for this algebra.
54 Element
= FiniteDimensionalEJAElement
68 if not field
.is_subring(RR
):
69 # Note: this does return true for the real algebraic
70 # field, the rationals, and any quadratic field where
71 # we've specified a real embedding.
72 raise ValueError("scalar field is not real")
74 # If the basis given to us wasn't over the field that it's
75 # supposed to be over, fix that. Or, you know, crash.
76 basis
= tuple( b
.change_ring(field
) for b
in basis
)
79 # Check commutativity of the Jordan and inner-products.
80 # This has to be done before we build the multiplication
81 # and inner-product tables/matrices, because we take
82 # advantage of symmetry in the process.
83 if not all( jordan_product(bi
,bj
) == jordan_product(bj
,bi
)
86 raise ValueError("Jordan product is not commutative")
88 if not all( inner_product(bi
,bj
) == inner_product(bj
,bi
)
91 raise ValueError("inner-product is not commutative")
94 category
= MagmaticAlgebras(field
).FiniteDimensional()
95 category
= category
.WithBasis().Unital()
97 # Element subalgebras can take advantage of this.
98 category
= category
.Associative()
100 # Call the superclass constructor so that we can use its from_vector()
101 # method to build our multiplication table.
103 super().__init
__(field
,
109 # Now comes all of the hard work. We'll be constructing an
110 # ambient vector space V that our (vectorized) basis lives in,
111 # as well as a subspace W of V spanned by those (vectorized)
112 # basis elements. The W-coordinates are the coefficients that
113 # we see in things like x = 1*e1 + 2*e2.
116 from sage
.structure
.element
import is_Matrix
117 basis_is_matrices
= False
121 if is_Matrix(basis
[0]):
122 if basis
[0].is_square():
123 # TODO: this ugly is_square() hack works around the problem
124 # of passing to_matrix()ed vectors in as the basis from a
125 # subalgebra. They aren't REALLY matrices, at least not of
126 # the type that we assume here... Ugh.
127 basis_is_matrices
= True
128 from mjo
.eja
.eja_utils
import _vec2mat
129 vector_basis
= tuple( map(_mat2vec
,basis
) )
130 degree
= basis
[0].nrows()**2
132 # convert from column matrices to vectors, yuck
133 basis
= tuple( map(_mat2vec
,basis
) )
135 degree
= basis
[0].degree()
137 degree
= basis
[0].degree()
139 # Build an ambient space that fits...
140 V
= VectorSpace(field
, degree
)
142 # We overwrite the name "vector_basis" in a second, but never modify it
143 # in place, to this effectively makes a copy of it.
144 deortho_vector_basis
= vector_basis
145 self
._deortho
_matrix
= None
148 from mjo
.eja
.eja_utils
import gram_schmidt
149 if basis_is_matrices
:
150 vector_ip
= lambda x
,y
: inner_product(_vec2mat(x
), _vec2mat(y
))
151 vector_basis
= gram_schmidt(vector_basis
, vector_ip
)
153 vector_basis
= gram_schmidt(vector_basis
, inner_product
)
155 # Normalize the "matrix" basis, too!
158 if basis_is_matrices
:
159 basis
= tuple( map(_vec2mat
,basis
) )
161 # Save the matrix "basis" for later... this is the last time we'll
162 # reference it in this constructor.
163 if basis_is_matrices
:
164 self
._matrix
_basis
= basis
166 MS
= MatrixSpace(self
.base_ring(), degree
, 1)
167 self
._matrix
_basis
= tuple( MS(b
) for b
in basis
)
169 # Now create the vector space for the algebra...
170 W
= V
.span_of_basis( vector_basis
, check
=check_axioms
)
173 # Now "W" is the vector space of our algebra coordinates. The
174 # variables "X1", "X2",... refer to the entries of vectors in
175 # W. Thus to convert back and forth between the orthonormal
176 # coordinates and the given ones, we need to stick the original
178 U
= V
.span_of_basis( deortho_vector_basis
, check
=check_axioms
)
179 self
._deortho
_matrix
= matrix( U
.coordinate_vector(q
)
180 for q
in vector_basis
)
183 # Now we actually compute the multiplication and inner-product
184 # tables/matrices using the possibly-orthonormalized basis.
185 self
._inner
_product
_matrix
= matrix
.zero(field
, n
)
186 self
._multiplication
_table
= [ [0 for j
in range(i
+1)] for i
in range(n
) ]
188 print("vector_basis:")
190 # Note: the Jordan and inner-products are defined in terms
191 # of the ambient basis. It's important that their arguments
192 # are in ambient coordinates as well.
195 # ortho basis w.r.t. ambient coords
196 q_i
= vector_basis
[i
]
197 q_j
= vector_basis
[j
]
199 if basis_is_matrices
:
203 elt
= jordan_product(q_i
, q_j
)
204 ip
= inner_product(q_i
, q_j
)
206 if basis_is_matrices
:
207 # do another mat2vec because the multiplication
208 # table is in terms of vectors
211 # TODO: the jordan product turns things back into
212 # matrices here even if they're supposed to be
213 # vectors. ugh. Can we get rid of vectors all together
215 elt
= W
.coordinate_vector(elt
)
216 self
._multiplication
_table
[i
][j
] = self
.from_vector(elt
)
217 self
._inner
_product
_matrix
[i
,j
] = ip
218 self
._inner
_product
_matrix
[j
,i
] = ip
220 self
._inner
_product
_matrix
._cache
= {'hermitian': True}
221 self
._inner
_product
_matrix
.set_immutable()
224 if not self
._is
_jordanian
():
225 raise ValueError("Jordan identity does not hold")
226 if not self
._inner
_product
_is
_associative
():
227 raise ValueError("inner product is not associative")
230 def _coerce_map_from_base_ring(self
):
232 Disable the map from the base ring into the algebra.
234 Performing a nonsense conversion like this automatically
235 is counterpedagogical. The fallback is to try the usual
236 element constructor, which should also fail.
240 sage: from mjo.eja.eja_algebra import random_eja
244 sage: set_random_seed()
245 sage: J = random_eja()
247 Traceback (most recent call last):
249 ValueError: not an element of this algebra
255 def product_on_basis(self
, i
, j
):
256 # We only stored the lower-triangular portion of the
257 # multiplication table.
259 return self
._multiplication
_table
[i
][j
]
261 return self
._multiplication
_table
[j
][i
]
263 def inner_product(self
, x
, y
):
265 The inner product associated with this Euclidean Jordan algebra.
267 Defaults to the trace inner product, but can be overridden by
268 subclasses if they are sure that the necessary properties are
273 sage: from mjo.eja.eja_algebra import (random_eja,
275 ....: BilinearFormEJA)
279 Our inner product is "associative," which means the following for
280 a symmetric bilinear form::
282 sage: set_random_seed()
283 sage: J = random_eja()
284 sage: x,y,z = J.random_elements(3)
285 sage: (x*y).inner_product(z) == y.inner_product(x*z)
290 Ensure that this is the usual inner product for the algebras
293 sage: set_random_seed()
294 sage: J = HadamardEJA.random_instance()
295 sage: x,y = J.random_elements(2)
296 sage: actual = x.inner_product(y)
297 sage: expected = x.to_vector().inner_product(y.to_vector())
298 sage: actual == expected
301 Ensure that this is one-half of the trace inner-product in a
302 BilinearFormEJA that isn't just the reals (when ``n`` isn't
303 one). This is in Faraut and Koranyi, and also my "On the
306 sage: set_random_seed()
307 sage: J = BilinearFormEJA.random_instance()
308 sage: n = J.dimension()
309 sage: x = J.random_element()
310 sage: y = J.random_element()
311 sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
314 B
= self
._inner
_product
_matrix
315 return (B
*x
.to_vector()).inner_product(y
.to_vector())
318 def _is_commutative(self
):
320 Whether or not this algebra's multiplication table is commutative.
322 This method should of course always return ``True``, unless
323 this algebra was constructed with ``check_axioms=False`` and
324 passed an invalid multiplication table.
326 return all( self
.product_on_basis(i
,j
) == self
.product_on_basis(i
,j
)
327 for i
in range(self
.dimension())
328 for j
in range(self
.dimension()) )
330 def _is_jordanian(self
):
332 Whether or not this algebra's multiplication table respects the
333 Jordan identity `(x^{2})(xy) = x(x^{2}y)`.
335 We only check one arrangement of `x` and `y`, so for a
336 ``True`` result to be truly true, you should also check
337 :meth:`_is_commutative`. This method should of course always
338 return ``True``, unless this algebra was constructed with
339 ``check_axioms=False`` and passed an invalid multiplication table.
341 return all( (self
.monomial(i
)**2)*(self
.monomial(i
)*self
.monomial(j
))
343 (self
.monomial(i
))*((self
.monomial(i
)**2)*self
.monomial(j
))
344 for i
in range(self
.dimension())
345 for j
in range(self
.dimension()) )
347 def _inner_product_is_associative(self
):
349 Return whether or not this algebra's inner product `B` is
350 associative; that is, whether or not `B(xy,z) = B(x,yz)`.
352 This method should of course always return ``True``, unless
353 this algebra was constructed with ``check_axioms=False`` and
354 passed an invalid multiplication table.
357 # Used to check whether or not something is zero in an inexact
358 # ring. This number is sufficient to allow the construction of
359 # QuaternionHermitianEJA(2, field=RDF) with check_axioms=True.
362 for i
in range(self
.dimension()):
363 for j
in range(self
.dimension()):
364 for k
in range(self
.dimension()):
368 diff
= (x
*y
).inner_product(z
) - x
.inner_product(y
*z
)
370 if self
.base_ring().is_exact():
374 if diff
.abs() > epsilon
:
379 def _element_constructor_(self
, elt
):
381 Construct an element of this algebra from its vector or matrix
384 This gets called only after the parent element _call_ method
385 fails to find a coercion for the argument.
389 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
391 ....: RealSymmetricEJA)
395 The identity in `S^n` is converted to the identity in the EJA::
397 sage: J = RealSymmetricEJA(3)
398 sage: I = matrix.identity(QQ,3)
399 sage: J(I) == J.one()
402 This skew-symmetric matrix can't be represented in the EJA::
404 sage: J = RealSymmetricEJA(3)
405 sage: A = matrix(QQ,3, lambda i,j: i-j)
407 Traceback (most recent call last):
409 ValueError: not an element of this algebra
413 Ensure that we can convert any element of the two non-matrix
414 simple algebras (whose matrix representations are columns)
415 back and forth faithfully::
417 sage: set_random_seed()
418 sage: J = HadamardEJA.random_instance()
419 sage: x = J.random_element()
420 sage: J(x.to_vector().column()) == x
422 sage: J = JordanSpinEJA.random_instance()
423 sage: x = J.random_element()
424 sage: J(x.to_vector().column()) == x
428 msg
= "not an element of this algebra"
430 # The superclass implementation of random_element()
431 # needs to be able to coerce "0" into the algebra.
433 elif elt
in self
.base_ring():
434 # Ensure that no base ring -> algebra coercion is performed
435 # by this method. There's some stupidity in sage that would
436 # otherwise propagate to this method; for example, sage thinks
437 # that the integer 3 belongs to the space of 2-by-2 matrices.
438 raise ValueError(msg
)
442 except (AttributeError, TypeError):
443 # Try to convert a vector into a column-matrix
446 if elt
not in self
.matrix_space():
447 raise ValueError(msg
)
449 # Thanks for nothing! Matrix spaces aren't vector spaces in
450 # Sage, so we have to figure out its matrix-basis coordinates
451 # ourselves. We use the basis space's ring instead of the
452 # element's ring because the basis space might be an algebraic
453 # closure whereas the base ring of the 3-by-3 identity matrix
454 # could be QQ instead of QQbar.
456 # We pass check=False because the matrix basis is "guaranteed"
457 # to be linearly independent... right? Ha ha.
458 V
= VectorSpace(self
.base_ring(), elt
.nrows()*elt
.ncols())
459 W
= V
.span_of_basis( (_mat2vec(s
) for s
in self
.matrix_basis()),
463 coords
= W
.coordinate_vector(_mat2vec(elt
))
464 except ArithmeticError: # vector is not in free module
465 raise ValueError(msg
)
467 return self
.from_vector(coords
)
471 Return a string representation of ``self``.
475 sage: from mjo.eja.eja_algebra import JordanSpinEJA
479 Ensure that it says what we think it says::
481 sage: JordanSpinEJA(2, field=AA)
482 Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
483 sage: JordanSpinEJA(3, field=RDF)
484 Euclidean Jordan algebra of dimension 3 over Real Double Field
487 fmt
= "Euclidean Jordan algebra of dimension {} over {}"
488 return fmt
.format(self
.dimension(), self
.base_ring())
492 def characteristic_polynomial_of(self
):
494 Return the algebra's "characteristic polynomial of" function,
495 which is itself a multivariate polynomial that, when evaluated
496 at the coordinates of some algebra element, returns that
497 element's characteristic polynomial.
499 The resulting polynomial has `n+1` variables, where `n` is the
500 dimension of this algebra. The first `n` variables correspond to
501 the coordinates of an algebra element: when evaluated at the
502 coordinates of an algebra element with respect to a certain
503 basis, the result is a univariate polynomial (in the one
504 remaining variable ``t``), namely the characteristic polynomial
509 sage: from mjo.eja.eja_algebra import JordanSpinEJA, TrivialEJA
513 The characteristic polynomial in the spin algebra is given in
514 Alizadeh, Example 11.11::
516 sage: J = JordanSpinEJA(3)
517 sage: p = J.characteristic_polynomial_of(); p
518 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
519 sage: xvec = J.one().to_vector()
523 By definition, the characteristic polynomial is a monic
524 degree-zero polynomial in a rank-zero algebra. Note that
525 Cayley-Hamilton is indeed satisfied since the polynomial
526 ``1`` evaluates to the identity element of the algebra on
529 sage: J = TrivialEJA()
530 sage: J.characteristic_polynomial_of()
537 # The list of coefficient polynomials a_0, a_1, a_2, ..., a_(r-1).
538 a
= self
._charpoly
_coefficients
()
540 # We go to a bit of trouble here to reorder the
541 # indeterminates, so that it's easier to evaluate the
542 # characteristic polynomial at x's coordinates and get back
543 # something in terms of t, which is what we want.
544 S
= PolynomialRing(self
.base_ring(),'t')
548 S
= PolynomialRing(S
, R
.variable_names())
551 return (t
**r
+ sum( a
[k
]*(t
**k
) for k
in range(r
) ))
553 def coordinate_polynomial_ring(self
):
555 The multivariate polynomial ring in which this algebra's
556 :meth:`characteristic_polynomial_of` lives.
560 sage: from mjo.eja.eja_algebra import (HadamardEJA,
561 ....: RealSymmetricEJA)
565 sage: J = HadamardEJA(2)
566 sage: J.coordinate_polynomial_ring()
567 Multivariate Polynomial Ring in X1, X2...
568 sage: J = RealSymmetricEJA(3,field=QQ,orthonormalize=False)
569 sage: J.coordinate_polynomial_ring()
570 Multivariate Polynomial Ring in X1, X2, X3, X4, X5, X6...
573 var_names
= tuple( "X%d" % z
for z
in range(1, self
.dimension()+1) )
574 return PolynomialRing(self
.base_ring(), var_names
)
576 def inner_product(self
, x
, y
):
578 The inner product associated with this Euclidean Jordan algebra.
580 Defaults to the trace inner product, but can be overridden by
581 subclasses if they are sure that the necessary properties are
586 sage: from mjo.eja.eja_algebra import (random_eja,
588 ....: BilinearFormEJA)
592 Our inner product is "associative," which means the following for
593 a symmetric bilinear form::
595 sage: set_random_seed()
596 sage: J = random_eja()
597 sage: x,y,z = J.random_elements(3)
598 sage: (x*y).inner_product(z) == y.inner_product(x*z)
603 Ensure that this is the usual inner product for the algebras
606 sage: set_random_seed()
607 sage: J = HadamardEJA.random_instance()
608 sage: x,y = J.random_elements(2)
609 sage: actual = x.inner_product(y)
610 sage: expected = x.to_vector().inner_product(y.to_vector())
611 sage: actual == expected
614 Ensure that this is one-half of the trace inner-product in a
615 BilinearFormEJA that isn't just the reals (when ``n`` isn't
616 one). This is in Faraut and Koranyi, and also my "On the
619 sage: set_random_seed()
620 sage: J = BilinearFormEJA.random_instance()
621 sage: n = J.dimension()
622 sage: x = J.random_element()
623 sage: y = J.random_element()
624 sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
627 B
= self
._inner
_product
_matrix
628 return (B
*x
.to_vector()).inner_product(y
.to_vector())
631 def is_trivial(self
):
633 Return whether or not this algebra is trivial.
635 A trivial algebra contains only the zero element.
639 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
644 sage: J = ComplexHermitianEJA(3)
650 sage: J = TrivialEJA()
655 return self
.dimension() == 0
658 def multiplication_table(self
):
660 Return a visual representation of this algebra's multiplication
661 table (on basis elements).
665 sage: from mjo.eja.eja_algebra import JordanSpinEJA
669 sage: J = JordanSpinEJA(4)
670 sage: J.multiplication_table()
671 +----++----+----+----+----+
672 | * || e0 | e1 | e2 | e3 |
673 +====++====+====+====+====+
674 | e0 || e0 | e1 | e2 | e3 |
675 +----++----+----+----+----+
676 | e1 || e1 | e0 | 0 | 0 |
677 +----++----+----+----+----+
678 | e2 || e2 | 0 | e0 | 0 |
679 +----++----+----+----+----+
680 | e3 || e3 | 0 | 0 | e0 |
681 +----++----+----+----+----+
685 # Prepend the header row.
686 M
= [["*"] + list(self
.gens())]
688 # And to each subsequent row, prepend an entry that belongs to
689 # the left-side "header column."
690 M
+= [ [self
.monomial(i
)] + [ self
.product_on_basis(i
,j
)
694 return table(M
, header_row
=True, header_column
=True, frame
=True)
697 def matrix_basis(self
):
699 Return an (often more natural) representation of this algebras
700 basis as an ordered tuple of matrices.
702 Every finite-dimensional Euclidean Jordan Algebra is a, up to
703 Jordan isomorphism, a direct sum of five simple
704 algebras---four of which comprise Hermitian matrices. And the
705 last type of algebra can of course be thought of as `n`-by-`1`
706 column matrices (ambiguusly called column vectors) to avoid
707 special cases. As a result, matrices (and column vectors) are
708 a natural representation format for Euclidean Jordan algebra
711 But, when we construct an algebra from a basis of matrices,
712 those matrix representations are lost in favor of coordinate
713 vectors *with respect to* that basis. We could eventually
714 convert back if we tried hard enough, but having the original
715 representations handy is valuable enough that we simply store
716 them and return them from this method.
718 Why implement this for non-matrix algebras? Avoiding special
719 cases for the :class:`BilinearFormEJA` pays with simplicity in
720 its own right. But mainly, we would like to be able to assume
721 that elements of a :class:`DirectSumEJA` can be displayed
722 nicely, without having to have special classes for direct sums
723 one of whose components was a matrix algebra.
727 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
728 ....: RealSymmetricEJA)
732 sage: J = RealSymmetricEJA(2)
734 Finite family {0: e0, 1: e1, 2: e2}
735 sage: J.matrix_basis()
737 [1 0] [ 0 0.7071067811865475?] [0 0]
738 [0 0], [0.7071067811865475? 0], [0 1]
743 sage: J = JordanSpinEJA(2)
745 Finite family {0: e0, 1: e1}
746 sage: J.matrix_basis()
752 return self
._matrix
_basis
755 def matrix_space(self
):
757 Return the matrix space in which this algebra's elements live, if
758 we think of them as matrices (including column vectors of the
761 Generally this will be an `n`-by-`1` column-vector space,
762 except when the algebra is trivial. There it's `n`-by-`n`
763 (where `n` is zero), to ensure that two elements of the matrix
764 space (empty matrices) can be multiplied.
766 Matrix algebras override this with something more useful.
768 if self
.is_trivial():
769 return MatrixSpace(self
.base_ring(), 0)
771 return self
._matrix
_basis
[0].matrix_space()
777 Return the unit element of this algebra.
781 sage: from mjo.eja.eja_algebra import (HadamardEJA,
786 sage: J = HadamardEJA(5)
788 e0 + e1 + e2 + e3 + e4
792 The identity element acts like the identity::
794 sage: set_random_seed()
795 sage: J = random_eja()
796 sage: x = J.random_element()
797 sage: J.one()*x == x and x*J.one() == x
800 The matrix of the unit element's operator is the identity::
802 sage: set_random_seed()
803 sage: J = random_eja()
804 sage: actual = J.one().operator().matrix()
805 sage: expected = matrix.identity(J.base_ring(), J.dimension())
806 sage: actual == expected
809 Ensure that the cached unit element (often precomputed by
810 hand) agrees with the computed one::
812 sage: set_random_seed()
813 sage: J = random_eja()
814 sage: cached = J.one()
815 sage: J.one.clear_cache()
816 sage: J.one() == cached
820 # We can brute-force compute the matrices of the operators
821 # that correspond to the basis elements of this algebra.
822 # If some linear combination of those basis elements is the
823 # algebra identity, then the same linear combination of
824 # their matrices has to be the identity matrix.
826 # Of course, matrices aren't vectors in sage, so we have to
827 # appeal to the "long vectors" isometry.
828 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
830 # Now we use basic linear algebra to find the coefficients,
831 # of the matrices-as-vectors-linear-combination, which should
832 # work for the original algebra basis too.
833 A
= matrix(self
.base_ring(), oper_vecs
)
835 # We used the isometry on the left-hand side already, but we
836 # still need to do it for the right-hand side. Recall that we
837 # wanted something that summed to the identity matrix.
838 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
840 # Now if there's an identity element in the algebra, this
841 # should work. We solve on the left to avoid having to
842 # transpose the matrix "A".
843 return self
.from_vector(A
.solve_left(b
))
846 def peirce_decomposition(self
, c
):
848 The Peirce decomposition of this algebra relative to the
851 In the future, this can be extended to a complete system of
852 orthogonal idempotents.
856 - ``c`` -- an idempotent of this algebra.
860 A triple (J0, J5, J1) containing two subalgebras and one subspace
863 - ``J0`` -- the algebra on the eigenspace of ``c.operator()``
864 corresponding to the eigenvalue zero.
866 - ``J5`` -- the eigenspace (NOT a subalgebra) of ``c.operator()``
867 corresponding to the eigenvalue one-half.
869 - ``J1`` -- the algebra on the eigenspace of ``c.operator()``
870 corresponding to the eigenvalue one.
872 These are the only possible eigenspaces for that operator, and this
873 algebra is a direct sum of them. The spaces ``J0`` and ``J1`` are
874 orthogonal, and are subalgebras of this algebra with the appropriate
879 sage: from mjo.eja.eja_algebra import random_eja, RealSymmetricEJA
883 The canonical example comes from the symmetric matrices, which
884 decompose into diagonal and off-diagonal parts::
886 sage: J = RealSymmetricEJA(3)
887 sage: C = matrix(QQ, [ [1,0,0],
891 sage: J0,J5,J1 = J.peirce_decomposition(c)
893 Euclidean Jordan algebra of dimension 1...
895 Vector space of degree 6 and dimension 2...
897 Euclidean Jordan algebra of dimension 3...
898 sage: J0.one().to_matrix()
902 sage: orig_df = AA.options.display_format
903 sage: AA.options.display_format = 'radical'
904 sage: J.from_vector(J5.basis()[0]).to_matrix()
908 sage: J.from_vector(J5.basis()[1]).to_matrix()
912 sage: AA.options.display_format = orig_df
913 sage: J1.one().to_matrix()
920 Every algebra decomposes trivially with respect to its identity
923 sage: set_random_seed()
924 sage: J = random_eja()
925 sage: J0,J5,J1 = J.peirce_decomposition(J.one())
926 sage: J0.dimension() == 0 and J5.dimension() == 0
928 sage: J1.superalgebra() == J and J1.dimension() == J.dimension()
931 The decomposition is into eigenspaces, and its components are
932 therefore necessarily orthogonal. Moreover, the identity
933 elements in the two subalgebras are the projections onto their
934 respective subspaces of the superalgebra's identity element::
936 sage: set_random_seed()
937 sage: J = random_eja()
938 sage: x = J.random_element()
939 sage: if not J.is_trivial():
940 ....: while x.is_nilpotent():
941 ....: x = J.random_element()
942 sage: c = x.subalgebra_idempotent()
943 sage: J0,J5,J1 = J.peirce_decomposition(c)
945 sage: for (w,y,z) in zip(J0.basis(), J5.basis(), J1.basis()):
946 ....: w = w.superalgebra_element()
947 ....: y = J.from_vector(y)
948 ....: z = z.superalgebra_element()
949 ....: ipsum += w.inner_product(y).abs()
950 ....: ipsum += w.inner_product(z).abs()
951 ....: ipsum += y.inner_product(z).abs()
954 sage: J1(c) == J1.one()
956 sage: J0(J.one() - c) == J0.one()
960 if not c
.is_idempotent():
961 raise ValueError("element is not idempotent: %s" % c
)
963 from mjo
.eja
.eja_subalgebra
import FiniteDimensionalEJASubalgebra
965 # Default these to what they should be if they turn out to be
966 # trivial, because eigenspaces_left() won't return eigenvalues
967 # corresponding to trivial spaces (e.g. it returns only the
968 # eigenspace corresponding to lambda=1 if you take the
969 # decomposition relative to the identity element).
970 trivial
= FiniteDimensionalEJASubalgebra(self
, ())
971 J0
= trivial
# eigenvalue zero
972 J5
= VectorSpace(self
.base_ring(), 0) # eigenvalue one-half
973 J1
= trivial
# eigenvalue one
975 for (eigval
, eigspace
) in c
.operator().matrix().right_eigenspaces():
976 if eigval
== ~
(self
.base_ring()(2)):
979 gens
= tuple( self
.from_vector(b
) for b
in eigspace
.basis() )
980 subalg
= FiniteDimensionalEJASubalgebra(self
,
988 raise ValueError("unexpected eigenvalue: %s" % eigval
)
993 def random_element(self
, thorough
=False):
995 Return a random element of this algebra.
997 Our algebra superclass method only returns a linear
998 combination of at most two basis elements. We instead
999 want the vector space "random element" method that
1000 returns a more diverse selection.
1004 - ``thorough`` -- (boolean; default False) whether or not we
1005 should generate irrational coefficients for the random
1006 element when our base ring is irrational; this slows the
1007 algebra operations to a crawl, but any truly random method
1011 # For a general base ring... maybe we can trust this to do the
1012 # right thing? Unlikely, but.
1013 V
= self
.vector_space()
1014 v
= V
.random_element()
1016 if self
.base_ring() is AA
:
1017 # The "random element" method of the algebraic reals is
1018 # stupid at the moment, and only returns integers between
1019 # -2 and 2, inclusive:
1021 # https://trac.sagemath.org/ticket/30875
1023 # Instead, we implement our own "random vector" method,
1024 # and then coerce that into the algebra. We use the vector
1025 # space degree here instead of the dimension because a
1026 # subalgebra could (for example) be spanned by only two
1027 # vectors, each with five coordinates. We need to
1028 # generate all five coordinates.
1030 v
*= QQbar
.random_element().real()
1032 v
*= QQ
.random_element()
1034 return self
.from_vector(V
.coordinate_vector(v
))
1036 def random_elements(self
, count
, thorough
=False):
1038 Return ``count`` random elements as a tuple.
1042 - ``thorough`` -- (boolean; default False) whether or not we
1043 should generate irrational coefficients for the random
1044 elements when our base ring is irrational; this slows the
1045 algebra operations to a crawl, but any truly random method
1050 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1054 sage: J = JordanSpinEJA(3)
1055 sage: x,y,z = J.random_elements(3)
1056 sage: all( [ x in J, y in J, z in J ])
1058 sage: len( J.random_elements(10) ) == 10
1062 return tuple( self
.random_element(thorough
)
1063 for idx
in range(count
) )
1067 def _charpoly_coefficients(self
):
1069 The `r` polynomial coefficients of the "characteristic polynomial
1072 n
= self
.dimension()
1073 R
= self
.coordinate_polynomial_ring()
1075 F
= R
.fraction_field()
1078 # From a result in my book, these are the entries of the
1079 # basis representation of L_x.
1080 return sum( vars[k
]*self
.monomial(k
).operator().matrix()[i
,j
]
1083 L_x
= matrix(F
, n
, n
, L_x_i_j
)
1086 if self
.rank
.is_in_cache():
1088 # There's no need to pad the system with redundant
1089 # columns if we *know* they'll be redundant.
1092 # Compute an extra power in case the rank is equal to
1093 # the dimension (otherwise, we would stop at x^(r-1)).
1094 x_powers
= [ (L_x
**k
)*self
.one().to_vector()
1095 for k
in range(n
+1) ]
1096 A
= matrix
.column(F
, x_powers
[:n
])
1097 AE
= A
.extended_echelon_form()
1104 # The theory says that only the first "r" coefficients are
1105 # nonzero, and they actually live in the original polynomial
1106 # ring and not the fraction field. We negate them because
1107 # in the actual characteristic polynomial, they get moved
1108 # to the other side where x^r lives.
1109 return -A_rref
.solve_right(E
*b
).change_ring(R
)[:r
]
1114 Return the rank of this EJA.
1116 This is a cached method because we know the rank a priori for
1117 all of the algebras we can construct. Thus we can avoid the
1118 expensive ``_charpoly_coefficients()`` call unless we truly
1119 need to compute the whole characteristic polynomial.
1123 sage: from mjo.eja.eja_algebra import (HadamardEJA,
1124 ....: JordanSpinEJA,
1125 ....: RealSymmetricEJA,
1126 ....: ComplexHermitianEJA,
1127 ....: QuaternionHermitianEJA,
1132 The rank of the Jordan spin algebra is always two::
1134 sage: JordanSpinEJA(2).rank()
1136 sage: JordanSpinEJA(3).rank()
1138 sage: JordanSpinEJA(4).rank()
1141 The rank of the `n`-by-`n` Hermitian real, complex, or
1142 quaternion matrices is `n`::
1144 sage: RealSymmetricEJA(4).rank()
1146 sage: ComplexHermitianEJA(3).rank()
1148 sage: QuaternionHermitianEJA(2).rank()
1153 Ensure that every EJA that we know how to construct has a
1154 positive integer rank, unless the algebra is trivial in
1155 which case its rank will be zero::
1157 sage: set_random_seed()
1158 sage: J = random_eja()
1162 sage: r > 0 or (r == 0 and J.is_trivial())
1165 Ensure that computing the rank actually works, since the ranks
1166 of all simple algebras are known and will be cached by default::
1168 sage: set_random_seed() # long time
1169 sage: J = random_eja() # long time
1170 sage: caches = J.rank() # long time
1171 sage: J.rank.clear_cache() # long time
1172 sage: J.rank() == cached # long time
1176 return len(self
._charpoly
_coefficients
())
1179 def vector_space(self
):
1181 Return the vector space that underlies this algebra.
1185 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1189 sage: J = RealSymmetricEJA(2)
1190 sage: J.vector_space()
1191 Vector space of dimension 3 over...
1194 return self
.zero().to_vector().parent().ambient_vector_space()
1197 Element
= FiniteDimensionalEJAElement
1199 class RationalBasisEJA(FiniteDimensionalEJA
):
1201 New class for algebras whose supplied basis elements have all rational entries.
1205 sage: from mjo.eja.eja_algebra import BilinearFormEJA
1209 The supplied basis is orthonormalized by default::
1211 sage: B = matrix(QQ, [[1, 0, 0], [0, 25, -32], [0, -32, 41]])
1212 sage: J = BilinearFormEJA(B)
1213 sage: J.matrix_basis()
1226 orthonormalize
=True,
1232 # Abuse the check_field parameter to check that the entries of
1233 # out basis (in ambient coordinates) are in the field QQ.
1234 if not all( all(b_i
in QQ
for b_i
in b
.list()) for b
in basis
):
1235 raise TypeError("basis not rational")
1238 # There's no point in constructing the extra algebra if this
1239 # one is already rational.
1241 # Note: the same Jordan and inner-products work here,
1242 # because they are necessarily defined with respect to
1243 # ambient coordinates and not any particular basis.
1244 self
._rational
_algebra
= FiniteDimensionalEJA(
1249 orthonormalize
=False,
1254 super().__init
__(basis
,
1258 check_field
=check_field
,
1259 check_axioms
=check_axioms
,
1263 def _charpoly_coefficients(self
):
1267 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
1268 ....: JordanSpinEJA)
1272 The base ring of the resulting polynomial coefficients is what
1273 it should be, and not the rationals (unless the algebra was
1274 already over the rationals)::
1276 sage: J = JordanSpinEJA(3)
1277 sage: J._charpoly_coefficients()
1278 (X1^2 - X2^2 - X3^2, -2*X1)
1279 sage: a0 = J._charpoly_coefficients()[0]
1281 Algebraic Real Field
1282 sage: a0.base_ring()
1283 Algebraic Real Field
1286 if self
._rational
_algebra
is None:
1287 # There's no need to construct *another* algebra over the
1288 # rationals if this one is already over the
1289 # rationals. Likewise, if we never orthonormalized our
1290 # basis, we might as well just use the given one.
1291 return super()._charpoly
_coefficients
()
1293 # Do the computation over the rationals. The answer will be
1294 # the same, because all we've done is a change of basis.
1295 # Then, change back from QQ to our real base ring
1296 a
= ( a_i
.change_ring(self
.base_ring())
1297 for a_i
in self
._rational
_algebra
._charpoly
_coefficients
() )
1299 # Now convert the coordinate variables back to the
1300 # deorthonormalized ones.
1301 R
= self
.coordinate_polynomial_ring()
1302 from sage
.modules
.free_module_element
import vector
1303 X
= vector(R
, R
.gens())
1304 BX
= self
._deortho
_matrix
*X
1306 subs_dict
= { X[i]: BX[i] for i in range(len(X)) }
1307 return tuple( a_i
.subs(subs_dict
) for a_i
in a
)
1309 class ConcreteEJA(RationalBasisEJA
):
1311 A class for the Euclidean Jordan algebras that we know by name.
1313 These are the Jordan algebras whose basis, multiplication table,
1314 rank, and so on are known a priori. More to the point, they are
1315 the Euclidean Jordan algebras for which we are able to conjure up
1316 a "random instance."
1320 sage: from mjo.eja.eja_algebra import ConcreteEJA
1324 Our basis is normalized with respect to the algebra's inner
1325 product, unless we specify otherwise::
1327 sage: set_random_seed()
1328 sage: J = ConcreteEJA.random_instance()
1329 sage: all( b.norm() == 1 for b in J.gens() )
1332 Since our basis is orthonormal with respect to the algebra's inner
1333 product, and since we know that this algebra is an EJA, any
1334 left-multiplication operator's matrix will be symmetric because
1335 natural->EJA basis representation is an isometry and within the
1336 EJA the operator is self-adjoint by the Jordan axiom::
1338 sage: set_random_seed()
1339 sage: J = ConcreteEJA.random_instance()
1340 sage: x = J.random_element()
1341 sage: x.operator().is_self_adjoint()
1346 def _max_random_instance_size():
1348 Return an integer "size" that is an upper bound on the size of
1349 this algebra when it is used in a random test
1350 case. Unfortunately, the term "size" is ambiguous -- when
1351 dealing with `R^n` under either the Hadamard or Jordan spin
1352 product, the "size" refers to the dimension `n`. When dealing
1353 with a matrix algebra (real symmetric or complex/quaternion
1354 Hermitian), it refers to the size of the matrix, which is far
1355 less than the dimension of the underlying vector space.
1357 This method must be implemented in each subclass.
1359 raise NotImplementedError
1362 def random_instance(cls
, *args
, **kwargs
):
1364 Return a random instance of this type of algebra.
1366 This method should be implemented in each subclass.
1368 from sage
.misc
.prandom
import choice
1369 eja_class
= choice(cls
.__subclasses
__())
1371 # These all bubble up to the RationalBasisEJA superclass
1372 # constructor, so any (kw)args valid there are also valid
1374 return eja_class
.random_instance(*args
, **kwargs
)
1379 def dimension_over_reals():
1381 The dimension of this matrix's base ring over the reals.
1383 The reals are dimension one over themselves, obviously; that's
1384 just `\mathbb{R}^{1}`. Likewise, the complex numbers `a + bi`
1385 have dimension two. Finally, the quaternions have dimension
1386 four over the reals.
1388 This is used to determine the size of the matrix returned from
1389 :meth:`real_embed`, among other things.
1391 raise NotImplementedError
1394 def real_embed(cls
,M
):
1396 Embed the matrix ``M`` into a space of real matrices.
1398 The matrix ``M`` can have entries in any field at the moment:
1399 the real numbers, complex numbers, or quaternions. And although
1400 they are not a field, we can probably support octonions at some
1401 point, too. This function returns a real matrix that "acts like"
1402 the original with respect to matrix multiplication; i.e.
1404 real_embed(M*N) = real_embed(M)*real_embed(N)
1407 if M
.ncols() != M
.nrows():
1408 raise ValueError("the matrix 'M' must be square")
1413 def real_unembed(cls
,M
):
1415 The inverse of :meth:`real_embed`.
1417 if M
.ncols() != M
.nrows():
1418 raise ValueError("the matrix 'M' must be square")
1419 if not ZZ(M
.nrows()).mod(cls
.dimension_over_reals()).is_zero():
1420 raise ValueError("the matrix 'M' must be a real embedding")
1424 def jordan_product(X
,Y
):
1425 return (X
*Y
+ Y
*X
)/2
1428 def trace_inner_product(cls
,X
,Y
):
1430 Compute the trace inner-product of two real-embeddings.
1434 sage: from mjo.eja.eja_algebra import (RealSymmetricEJA,
1435 ....: ComplexHermitianEJA,
1436 ....: QuaternionHermitianEJA)
1440 This gives the same answer as it would if we computed the trace
1441 from the unembedded (original) matrices::
1443 sage: set_random_seed()
1444 sage: J = RealSymmetricEJA.random_instance()
1445 sage: x,y = J.random_elements(2)
1446 sage: Xe = x.to_matrix()
1447 sage: Ye = y.to_matrix()
1448 sage: X = J.real_unembed(Xe)
1449 sage: Y = J.real_unembed(Ye)
1450 sage: expected = (X*Y).trace()
1451 sage: actual = J.trace_inner_product(Xe,Ye)
1452 sage: actual == expected
1457 sage: set_random_seed()
1458 sage: J = ComplexHermitianEJA.random_instance()
1459 sage: x,y = J.random_elements(2)
1460 sage: Xe = x.to_matrix()
1461 sage: Ye = y.to_matrix()
1462 sage: X = J.real_unembed(Xe)
1463 sage: Y = J.real_unembed(Ye)
1464 sage: expected = (X*Y).trace().real()
1465 sage: actual = J.trace_inner_product(Xe,Ye)
1466 sage: actual == expected
1471 sage: set_random_seed()
1472 sage: J = QuaternionHermitianEJA.random_instance()
1473 sage: x,y = J.random_elements(2)
1474 sage: Xe = x.to_matrix()
1475 sage: Ye = y.to_matrix()
1476 sage: X = J.real_unembed(Xe)
1477 sage: Y = J.real_unembed(Ye)
1478 sage: expected = (X*Y).trace().coefficient_tuple()[0]
1479 sage: actual = J.trace_inner_product(Xe,Ye)
1480 sage: actual == expected
1484 Xu
= cls
.real_unembed(X
)
1485 Yu
= cls
.real_unembed(Y
)
1486 tr
= (Xu
*Yu
).trace()
1489 # Works in QQ, AA, RDF, et cetera.
1491 except AttributeError:
1492 # A quaternion doesn't have a real() method, but does
1493 # have coefficient_tuple() method that returns the
1494 # coefficients of 1, i, j, and k -- in that order.
1495 return tr
.coefficient_tuple()[0]
1498 class RealMatrixEJA(MatrixEJA
):
1500 def dimension_over_reals():
1504 class RealSymmetricEJA(ConcreteEJA
, RealMatrixEJA
):
1506 The rank-n simple EJA consisting of real symmetric n-by-n
1507 matrices, the usual symmetric Jordan product, and the trace inner
1508 product. It has dimension `(n^2 + n)/2` over the reals.
1512 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1516 sage: J = RealSymmetricEJA(2)
1517 sage: e0, e1, e2 = J.gens()
1525 In theory, our "field" can be any subfield of the reals::
1527 sage: RealSymmetricEJA(2, field=RDF)
1528 Euclidean Jordan algebra of dimension 3 over Real Double Field
1529 sage: RealSymmetricEJA(2, field=RR)
1530 Euclidean Jordan algebra of dimension 3 over Real Field with
1531 53 bits of precision
1535 The dimension of this algebra is `(n^2 + n) / 2`::
1537 sage: set_random_seed()
1538 sage: n_max = RealSymmetricEJA._max_random_instance_size()
1539 sage: n = ZZ.random_element(1, n_max)
1540 sage: J = RealSymmetricEJA(n)
1541 sage: J.dimension() == (n^2 + n)/2
1544 The Jordan multiplication is what we think it is::
1546 sage: set_random_seed()
1547 sage: J = RealSymmetricEJA.random_instance()
1548 sage: x,y = J.random_elements(2)
1549 sage: actual = (x*y).to_matrix()
1550 sage: X = x.to_matrix()
1551 sage: Y = y.to_matrix()
1552 sage: expected = (X*Y + Y*X)/2
1553 sage: actual == expected
1555 sage: J(expected) == x*y
1558 We can change the generator prefix::
1560 sage: RealSymmetricEJA(3, prefix='q').gens()
1561 (q0, q1, q2, q3, q4, q5)
1563 We can construct the (trivial) algebra of rank zero::
1565 sage: RealSymmetricEJA(0)
1566 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1570 def _denormalized_basis(cls
, n
):
1572 Return a basis for the space of real symmetric n-by-n matrices.
1576 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1580 sage: set_random_seed()
1581 sage: n = ZZ.random_element(1,5)
1582 sage: B = RealSymmetricEJA._denormalized_basis(n)
1583 sage: all( M.is_symmetric() for M in B)
1587 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1591 for j
in range(i
+1):
1592 Eij
= matrix(ZZ
, n
, lambda k
,l
: k
==i
and l
==j
)
1596 Sij
= Eij
+ Eij
.transpose()
1602 def _max_random_instance_size():
1603 return 4 # Dimension 10
1606 def random_instance(cls
, **kwargs
):
1608 Return a random instance of this type of algebra.
1610 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1611 return cls(n
, **kwargs
)
1613 def __init__(self
, n
, **kwargs
):
1614 # We know this is a valid EJA, but will double-check
1615 # if the user passes check_axioms=True.
1616 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
1618 super(RealSymmetricEJA
, self
).__init
__(self
._denormalized
_basis
(n
),
1619 self
.jordan_product
,
1620 self
.trace_inner_product
,
1623 # TODO: this could be factored out somehow, but is left here
1624 # because the MatrixEJA is not presently a subclass of the
1625 # FDEJA class that defines rank() and one().
1626 self
.rank
.set_cache(n
)
1627 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
1628 self
.one
.set_cache(self(idV
))
1632 class ComplexMatrixEJA(MatrixEJA
):
1634 def dimension_over_reals():
1638 def real_embed(cls
,M
):
1640 Embed the n-by-n complex matrix ``M`` into the space of real
1641 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1642 bi` to the block matrix ``[[a,b],[-b,a]]``.
1646 sage: from mjo.eja.eja_algebra import ComplexMatrixEJA
1650 sage: F = QuadraticField(-1, 'I')
1651 sage: x1 = F(4 - 2*i)
1652 sage: x2 = F(1 + 2*i)
1655 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1656 sage: ComplexMatrixEJA.real_embed(M)
1665 Embedding is a homomorphism (isomorphism, in fact)::
1667 sage: set_random_seed()
1668 sage: n = ZZ.random_element(3)
1669 sage: F = QuadraticField(-1, 'I')
1670 sage: X = random_matrix(F, n)
1671 sage: Y = random_matrix(F, n)
1672 sage: Xe = ComplexMatrixEJA.real_embed(X)
1673 sage: Ye = ComplexMatrixEJA.real_embed(Y)
1674 sage: XYe = ComplexMatrixEJA.real_embed(X*Y)
1679 super(ComplexMatrixEJA
,cls
).real_embed(M
)
1682 # We don't need any adjoined elements...
1683 field
= M
.base_ring().base_ring()
1687 a
= z
.list()[0] # real part, I guess
1688 b
= z
.list()[1] # imag part, I guess
1689 blocks
.append(matrix(field
, 2, [[a
,b
],[-b
,a
]]))
1691 return matrix
.block(field
, n
, blocks
)
1695 def real_unembed(cls
,M
):
1697 The inverse of _embed_complex_matrix().
1701 sage: from mjo.eja.eja_algebra import ComplexMatrixEJA
1705 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1706 ....: [-2, 1, -4, 3],
1707 ....: [ 9, 10, 11, 12],
1708 ....: [-10, 9, -12, 11] ])
1709 sage: ComplexMatrixEJA.real_unembed(A)
1711 [ 10*I + 9 12*I + 11]
1715 Unembedding is the inverse of embedding::
1717 sage: set_random_seed()
1718 sage: F = QuadraticField(-1, 'I')
1719 sage: M = random_matrix(F, 3)
1720 sage: Me = ComplexMatrixEJA.real_embed(M)
1721 sage: ComplexMatrixEJA.real_unembed(Me) == M
1725 super(ComplexMatrixEJA
,cls
).real_unembed(M
)
1727 d
= cls
.dimension_over_reals()
1729 # If "M" was normalized, its base ring might have roots
1730 # adjoined and they can stick around after unembedding.
1731 field
= M
.base_ring()
1732 R
= PolynomialRing(field
, 'z')
1735 # Sage doesn't know how to adjoin the complex "i" (the root of
1736 # x^2 + 1) to a field in a general way. Here, we just enumerate
1737 # all of the cases that I have cared to support so far.
1739 # Sage doesn't know how to embed AA into QQbar, i.e. how
1740 # to adjoin sqrt(-1) to AA.
1742 elif not field
.is_exact():
1744 F
= field
.complex_field()
1746 # Works for QQ and... maybe some other fields.
1747 F
= field
.extension(z
**2 + 1, 'I', embedding
=CLF(-1).sqrt())
1750 # Go top-left to bottom-right (reading order), converting every
1751 # 2-by-2 block we see to a single complex element.
1753 for k
in range(n
/d
):
1754 for j
in range(n
/d
):
1755 submat
= M
[d
*k
:d
*k
+d
,d
*j
:d
*j
+d
]
1756 if submat
[0,0] != submat
[1,1]:
1757 raise ValueError('bad on-diagonal submatrix')
1758 if submat
[0,1] != -submat
[1,0]:
1759 raise ValueError('bad off-diagonal submatrix')
1760 z
= submat
[0,0] + submat
[0,1]*i
1763 return matrix(F
, n
/d
, elements
)
1766 class ComplexHermitianEJA(ConcreteEJA
, ComplexMatrixEJA
):
1768 The rank-n simple EJA consisting of complex Hermitian n-by-n
1769 matrices over the real numbers, the usual symmetric Jordan product,
1770 and the real-part-of-trace inner product. It has dimension `n^2` over
1775 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1779 In theory, our "field" can be any subfield of the reals::
1781 sage: ComplexHermitianEJA(2, field=RDF)
1782 Euclidean Jordan algebra of dimension 4 over Real Double Field
1783 sage: ComplexHermitianEJA(2, field=RR)
1784 Euclidean Jordan algebra of dimension 4 over Real Field with
1785 53 bits of precision
1789 The dimension of this algebra is `n^2`::
1791 sage: set_random_seed()
1792 sage: n_max = ComplexHermitianEJA._max_random_instance_size()
1793 sage: n = ZZ.random_element(1, n_max)
1794 sage: J = ComplexHermitianEJA(n)
1795 sage: J.dimension() == n^2
1798 The Jordan multiplication is what we think it is::
1800 sage: set_random_seed()
1801 sage: J = ComplexHermitianEJA.random_instance()
1802 sage: x,y = J.random_elements(2)
1803 sage: actual = (x*y).to_matrix()
1804 sage: X = x.to_matrix()
1805 sage: Y = y.to_matrix()
1806 sage: expected = (X*Y + Y*X)/2
1807 sage: actual == expected
1809 sage: J(expected) == x*y
1812 We can change the generator prefix::
1814 sage: ComplexHermitianEJA(2, prefix='z').gens()
1817 We can construct the (trivial) algebra of rank zero::
1819 sage: ComplexHermitianEJA(0)
1820 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1825 def _denormalized_basis(cls
, n
):
1827 Returns a basis for the space of complex Hermitian n-by-n matrices.
1829 Why do we embed these? Basically, because all of numerical linear
1830 algebra assumes that you're working with vectors consisting of `n`
1831 entries from a field and scalars from the same field. There's no way
1832 to tell SageMath that (for example) the vectors contain complex
1833 numbers, while the scalar field is real.
1837 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1841 sage: set_random_seed()
1842 sage: n = ZZ.random_element(1,5)
1843 sage: field = QuadraticField(2, 'sqrt2')
1844 sage: B = ComplexHermitianEJA._denormalized_basis(n)
1845 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
, **kwargs
):
1880 # We know this is a valid EJA, but will double-check
1881 # if the user passes check_axioms=True.
1882 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
1884 super(ComplexHermitianEJA
, self
).__init
__(self
._denormalized
_basis
(n
),
1885 self
.jordan_product
,
1886 self
.trace_inner_product
,
1888 # TODO: this could be factored out somehow, but is left here
1889 # because the MatrixEJA is not presently a subclass of the
1890 # FDEJA class that defines rank() and one().
1891 self
.rank
.set_cache(n
)
1892 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
1893 self
.one
.set_cache(self(idV
))
1896 def _max_random_instance_size():
1897 return 3 # Dimension 9
1900 def random_instance(cls
, **kwargs
):
1902 Return a random instance of this type of algebra.
1904 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1905 return cls(n
, **kwargs
)
1907 class QuaternionMatrixEJA(MatrixEJA
):
1909 def dimension_over_reals():
1913 def real_embed(cls
,M
):
1915 Embed the n-by-n quaternion matrix ``M`` into the space of real
1916 matrices of size 4n-by-4n by first sending each quaternion entry `z
1917 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
1918 c+di],[-c + di, a-bi]]`, and then embedding those into a real
1923 sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA
1927 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1928 sage: i,j,k = Q.gens()
1929 sage: x = 1 + 2*i + 3*j + 4*k
1930 sage: M = matrix(Q, 1, [[x]])
1931 sage: QuaternionMatrixEJA.real_embed(M)
1937 Embedding is a homomorphism (isomorphism, in fact)::
1939 sage: set_random_seed()
1940 sage: n = ZZ.random_element(2)
1941 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1942 sage: X = random_matrix(Q, n)
1943 sage: Y = random_matrix(Q, n)
1944 sage: Xe = QuaternionMatrixEJA.real_embed(X)
1945 sage: Ye = QuaternionMatrixEJA.real_embed(Y)
1946 sage: XYe = QuaternionMatrixEJA.real_embed(X*Y)
1951 super(QuaternionMatrixEJA
,cls
).real_embed(M
)
1952 quaternions
= M
.base_ring()
1955 F
= QuadraticField(-1, 'I')
1960 t
= z
.coefficient_tuple()
1965 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
1966 [-c
+ d
*i
, a
- b
*i
]])
1967 realM
= ComplexMatrixEJA
.real_embed(cplxM
)
1968 blocks
.append(realM
)
1970 # We should have real entries by now, so use the realest field
1971 # we've got for the return value.
1972 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
1977 def real_unembed(cls
,M
):
1979 The inverse of _embed_quaternion_matrix().
1983 sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA
1987 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
1988 ....: [-2, 1, -4, 3],
1989 ....: [-3, 4, 1, -2],
1990 ....: [-4, -3, 2, 1]])
1991 sage: QuaternionMatrixEJA.real_unembed(M)
1992 [1 + 2*i + 3*j + 4*k]
1996 Unembedding is the inverse of embedding::
1998 sage: set_random_seed()
1999 sage: Q = QuaternionAlgebra(QQ, -1, -1)
2000 sage: M = random_matrix(Q, 3)
2001 sage: Me = QuaternionMatrixEJA.real_embed(M)
2002 sage: QuaternionMatrixEJA.real_unembed(Me) == M
2006 super(QuaternionMatrixEJA
,cls
).real_unembed(M
)
2008 d
= cls
.dimension_over_reals()
2010 # Use the base ring of the matrix to ensure that its entries can be
2011 # multiplied by elements of the quaternion algebra.
2012 field
= M
.base_ring()
2013 Q
= QuaternionAlgebra(field
,-1,-1)
2016 # Go top-left to bottom-right (reading order), converting every
2017 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
2020 for l
in range(n
/d
):
2021 for m
in range(n
/d
):
2022 submat
= ComplexMatrixEJA
.real_unembed(
2023 M
[d
*l
:d
*l
+d
,d
*m
:d
*m
+d
] )
2024 if submat
[0,0] != submat
[1,1].conjugate():
2025 raise ValueError('bad on-diagonal submatrix')
2026 if submat
[0,1] != -submat
[1,0].conjugate():
2027 raise ValueError('bad off-diagonal submatrix')
2028 z
= submat
[0,0].real()
2029 z
+= submat
[0,0].imag()*i
2030 z
+= submat
[0,1].real()*j
2031 z
+= submat
[0,1].imag()*k
2034 return matrix(Q
, n
/d
, elements
)
2037 class QuaternionHermitianEJA(ConcreteEJA
, QuaternionMatrixEJA
):
2039 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
2040 matrices, the usual symmetric Jordan product, and the
2041 real-part-of-trace inner product. It has dimension `2n^2 - n` over
2046 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2050 In theory, our "field" can be any subfield of the reals::
2052 sage: QuaternionHermitianEJA(2, field=RDF)
2053 Euclidean Jordan algebra of dimension 6 over Real Double Field
2054 sage: QuaternionHermitianEJA(2, field=RR)
2055 Euclidean Jordan algebra of dimension 6 over Real Field with
2056 53 bits of precision
2060 The dimension of this algebra is `2*n^2 - n`::
2062 sage: set_random_seed()
2063 sage: n_max = QuaternionHermitianEJA._max_random_instance_size()
2064 sage: n = ZZ.random_element(1, n_max)
2065 sage: J = QuaternionHermitianEJA(n)
2066 sage: J.dimension() == 2*(n^2) - n
2069 The Jordan multiplication is what we think it is::
2071 sage: set_random_seed()
2072 sage: J = QuaternionHermitianEJA.random_instance()
2073 sage: x,y = J.random_elements(2)
2074 sage: actual = (x*y).to_matrix()
2075 sage: X = x.to_matrix()
2076 sage: Y = y.to_matrix()
2077 sage: expected = (X*Y + Y*X)/2
2078 sage: actual == expected
2080 sage: J(expected) == x*y
2083 We can change the generator prefix::
2085 sage: QuaternionHermitianEJA(2, prefix='a').gens()
2086 (a0, a1, a2, a3, a4, a5)
2088 We can construct the (trivial) algebra of rank zero::
2090 sage: QuaternionHermitianEJA(0)
2091 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2095 def _denormalized_basis(cls
, n
):
2097 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
2099 Why do we embed these? Basically, because all of numerical
2100 linear algebra assumes that you're working with vectors consisting
2101 of `n` entries from a field and scalars from the same field. There's
2102 no way to tell SageMath that (for example) the vectors contain
2103 complex numbers, while the scalar field is real.
2107 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2111 sage: set_random_seed()
2112 sage: n = ZZ.random_element(1,5)
2113 sage: B = QuaternionHermitianEJA._denormalized_basis(n)
2114 sage: all( M.is_symmetric() for M in B )
2119 Q
= QuaternionAlgebra(QQ
,-1,-1)
2122 # This is like the symmetric case, but we need to be careful:
2124 # * We want conjugate-symmetry, not just symmetry.
2125 # * The diagonal will (as a result) be real.
2129 for j
in range(i
+1):
2130 Eij
= matrix(Q
, n
, lambda k
,l
: k
==i
and l
==j
)
2132 Sij
= cls
.real_embed(Eij
)
2135 # The second, third, and fourth ones have a minus
2136 # because they're conjugated.
2137 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
2139 Sij_I
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
2141 Sij_J
= cls
.real_embed(J
*Eij
- J
*Eij
.transpose())
2143 Sij_K
= cls
.real_embed(K
*Eij
- K
*Eij
.transpose())
2146 # Since we embedded these, we can drop back to the "field" that we
2147 # started with instead of the quaternion algebra "Q".
2148 return tuple( s
.change_ring(field
) for s
in S
)
2151 def __init__(self
, n
, **kwargs
):
2152 # We know this is a valid EJA, but will double-check
2153 # if the user passes check_axioms=True.
2154 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2156 super(QuaternionHermitianEJA
, self
).__init
__(self
._denormalized
_basis
(n
),
2157 self
.jordan_product
,
2158 self
.trace_inner_product
,
2160 # TODO: this could be factored out somehow, but is left here
2161 # because the MatrixEJA is not presently a subclass of the
2162 # FDEJA class that defines rank() and one().
2163 self
.rank
.set_cache(n
)
2164 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
2165 self
.one
.set_cache(self(idV
))
2169 def _max_random_instance_size():
2171 The maximum rank of a random QuaternionHermitianEJA.
2173 return 2 # Dimension 6
2176 def random_instance(cls
, **kwargs
):
2178 Return a random instance of this type of algebra.
2180 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2181 return cls(n
, **kwargs
)
2184 class HadamardEJA(ConcreteEJA
):
2186 Return the Euclidean Jordan Algebra corresponding to the set
2187 `R^n` under the Hadamard product.
2189 Note: this is nothing more than the Cartesian product of ``n``
2190 copies of the spin algebra. Once Cartesian product algebras
2191 are implemented, this can go.
2195 sage: from mjo.eja.eja_algebra import HadamardEJA
2199 This multiplication table can be verified by hand::
2201 sage: J = HadamardEJA(3)
2202 sage: e0,e1,e2 = J.gens()
2218 We can change the generator prefix::
2220 sage: HadamardEJA(3, prefix='r').gens()
2224 def __init__(self
, n
, **kwargs
):
2225 def jordan_product(x
,y
):
2227 return P(tuple( xi
*yi
for (xi
,yi
) in zip(x
,y
) ))
2228 def inner_product(x
,y
):
2229 return x
.inner_product(y
)
2231 # New defaults for keyword arguments. Don't orthonormalize
2232 # because our basis is already orthonormal with respect to our
2233 # inner-product. Don't check the axioms, because we know this
2234 # is a valid EJA... but do double-check if the user passes
2235 # check_axioms=True. Note: we DON'T override the "check_field"
2236 # default here, because the user can pass in a field!
2237 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2238 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2241 standard_basis
= FreeModule(ZZ
, n
).basis()
2242 super(HadamardEJA
, self
).__init
__(standard_basis
,
2246 self
.rank
.set_cache(n
)
2249 self
.one
.set_cache( self
.zero() )
2251 self
.one
.set_cache( sum(self
.gens()) )
2254 def _max_random_instance_size():
2256 The maximum dimension of a random HadamardEJA.
2261 def random_instance(cls
, **kwargs
):
2263 Return a random instance of this type of algebra.
2265 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2266 return cls(n
, **kwargs
)
2269 class BilinearFormEJA(ConcreteEJA
):
2271 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2272 with the half-trace inner product and jordan product ``x*y =
2273 (<Bx,y>,y_bar>, x0*y_bar + y0*x_bar)`` where `B = 1 \times B22` is
2274 a symmetric positive-definite "bilinear form" matrix. Its
2275 dimension is the size of `B`, and it has rank two in dimensions
2276 larger than two. It reduces to the ``JordanSpinEJA`` when `B` is
2277 the identity matrix of order ``n``.
2279 We insist that the one-by-one upper-left identity block of `B` be
2280 passed in as well so that we can be passed a matrix of size zero
2281 to construct a trivial algebra.
2285 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
2286 ....: JordanSpinEJA)
2290 When no bilinear form is specified, the identity matrix is used,
2291 and the resulting algebra is the Jordan spin algebra::
2293 sage: B = matrix.identity(AA,3)
2294 sage: J0 = BilinearFormEJA(B)
2295 sage: J1 = JordanSpinEJA(3)
2296 sage: J0.multiplication_table() == J0.multiplication_table()
2299 An error is raised if the matrix `B` does not correspond to a
2300 positive-definite bilinear form::
2302 sage: B = matrix.random(QQ,2,3)
2303 sage: J = BilinearFormEJA(B)
2304 Traceback (most recent call last):
2306 ValueError: bilinear form is not positive-definite
2307 sage: B = matrix.zero(QQ,3)
2308 sage: J = BilinearFormEJA(B)
2309 Traceback (most recent call last):
2311 ValueError: bilinear form is not positive-definite
2315 We can create a zero-dimensional algebra::
2317 sage: B = matrix.identity(AA,0)
2318 sage: J = BilinearFormEJA(B)
2322 We can check the multiplication condition given in the Jordan, von
2323 Neumann, and Wigner paper (and also discussed on my "On the
2324 symmetry..." paper). Note that this relies heavily on the standard
2325 choice of basis, as does anything utilizing the bilinear form
2326 matrix. We opt not to orthonormalize the basis, because if we
2327 did, we would have to normalize the `s_{i}` in a similar manner::
2329 sage: set_random_seed()
2330 sage: n = ZZ.random_element(5)
2331 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2332 sage: B11 = matrix.identity(QQ,1)
2333 sage: B22 = M.transpose()*M
2334 sage: B = block_matrix(2,2,[ [B11,0 ],
2336 sage: J = BilinearFormEJA(B, orthonormalize=False)
2337 sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
2338 sage: V = J.vector_space()
2339 sage: sis = [ J( V([0] + (M.inverse()*ei).list()).column() )
2340 ....: for ei in eis ]
2341 sage: actual = [ sis[i]*sis[j]
2342 ....: for i in range(n-1)
2343 ....: for j in range(n-1) ]
2344 sage: expected = [ J.one() if i == j else J.zero()
2345 ....: for i in range(n-1)
2346 ....: for j in range(n-1) ]
2347 sage: actual == expected
2350 def __init__(self
, B
, **kwargs
):
2351 if not B
.is_positive_definite():
2352 raise ValueError("bilinear form is not positive-definite")
2354 def inner_product(x
,y
):
2355 return (B
*x
).inner_product(y
)
2357 def jordan_product(x
,y
):
2363 z0
= inner_product(x
,y
)
2364 zbar
= y0
*xbar
+ x0
*ybar
2365 return P((z0
,) + tuple(zbar
))
2367 # We know this is a valid EJA, but will double-check
2368 # if the user passes check_axioms=True.
2369 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2372 standard_basis
= FreeModule(ZZ
, n
).basis()
2373 super(BilinearFormEJA
, self
).__init
__(standard_basis
,
2378 # The rank of this algebra is two, unless we're in a
2379 # one-dimensional ambient space (because the rank is bounded
2380 # by the ambient dimension).
2381 self
.rank
.set_cache(min(n
,2))
2384 self
.one
.set_cache( self
.zero() )
2386 self
.one
.set_cache( self
.monomial(0) )
2389 def _max_random_instance_size():
2391 The maximum dimension of a random BilinearFormEJA.
2396 def random_instance(cls
, **kwargs
):
2398 Return a random instance of this algebra.
2400 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2402 B
= matrix
.identity(ZZ
, n
)
2403 return cls(B
, **kwargs
)
2405 B11
= matrix
.identity(ZZ
, 1)
2406 M
= matrix
.random(ZZ
, n
-1)
2407 I
= matrix
.identity(ZZ
, n
-1)
2409 while alpha
.is_zero():
2410 alpha
= ZZ
.random_element().abs()
2411 B22
= M
.transpose()*M
+ alpha
*I
2413 from sage
.matrix
.special
import block_matrix
2414 B
= block_matrix(2,2, [ [B11
, ZZ(0) ],
2417 return cls(B
, **kwargs
)
2420 class JordanSpinEJA(BilinearFormEJA
):
2422 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2423 with the usual inner product and jordan product ``x*y =
2424 (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2429 sage: from mjo.eja.eja_algebra import JordanSpinEJA
2433 This multiplication table can be verified by hand::
2435 sage: J = JordanSpinEJA(4)
2436 sage: e0,e1,e2,e3 = J.gens()
2452 We can change the generator prefix::
2454 sage: JordanSpinEJA(2, prefix='B').gens()
2459 Ensure that we have the usual inner product on `R^n`::
2461 sage: set_random_seed()
2462 sage: J = JordanSpinEJA.random_instance()
2463 sage: x,y = J.random_elements(2)
2464 sage: actual = x.inner_product(y)
2465 sage: expected = x.to_vector().inner_product(y.to_vector())
2466 sage: actual == expected
2470 def __init__(self
, n
, **kwargs
):
2471 # This is a special case of the BilinearFormEJA with the
2472 # identity matrix as its bilinear form.
2473 B
= matrix
.identity(ZZ
, n
)
2475 # Don't orthonormalize because our basis is already
2476 # orthonormal with respect to our inner-product.
2477 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2479 # But also don't pass check_field=False here, because the user
2480 # can pass in a field!
2481 super(JordanSpinEJA
, self
).__init
__(B
, **kwargs
)
2484 def _max_random_instance_size():
2486 The maximum dimension of a random JordanSpinEJA.
2491 def random_instance(cls
, **kwargs
):
2493 Return a random instance of this type of algebra.
2495 Needed here to override the implementation for ``BilinearFormEJA``.
2497 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2498 return cls(n
, **kwargs
)
2501 class TrivialEJA(ConcreteEJA
):
2503 The trivial Euclidean Jordan algebra consisting of only a zero element.
2507 sage: from mjo.eja.eja_algebra import TrivialEJA
2511 sage: J = TrivialEJA()
2518 sage: 7*J.one()*12*J.one()
2520 sage: J.one().inner_product(J.one())
2522 sage: J.one().norm()
2524 sage: J.one().subalgebra_generated_by()
2525 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2530 def __init__(self
, **kwargs
):
2531 jordan_product
= lambda x
,y
: x
2532 inner_product
= lambda x
,y
: 0
2535 # New defaults for keyword arguments
2536 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2537 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2539 super(TrivialEJA
, self
).__init
__(basis
,
2543 # The rank is zero using my definition, namely the dimension of the
2544 # largest subalgebra generated by any element.
2545 self
.rank
.set_cache(0)
2546 self
.one
.set_cache( self
.zero() )
2549 def random_instance(cls
, **kwargs
):
2550 # We don't take a "size" argument so the superclass method is
2551 # inappropriate for us.
2552 return cls(**kwargs
)
2554 class DirectSumEJA(ConcreteEJA
):
2556 The external (orthogonal) direct sum of two other Euclidean Jordan
2557 algebras. Essentially the Cartesian product of its two factors.
2558 Every Euclidean Jordan algebra decomposes into an orthogonal
2559 direct sum of simple Euclidean Jordan algebras, so no generality
2560 is lost by providing only this construction.
2564 sage: from mjo.eja.eja_algebra import (random_eja,
2566 ....: RealSymmetricEJA,
2571 sage: J1 = HadamardEJA(2)
2572 sage: J2 = RealSymmetricEJA(3)
2573 sage: J = DirectSumEJA(J1,J2)
2581 The external direct sum construction is only valid when the two factors
2582 have the same base ring; an error is raised otherwise::
2584 sage: set_random_seed()
2585 sage: J1 = random_eja(field=AA)
2586 sage: J2 = random_eja(field=QQ,orthonormalize=False)
2587 sage: J = DirectSumEJA(J1,J2)
2588 Traceback (most recent call last):
2590 ValueError: algebras must share the same base field
2593 def __init__(self
, J1
, J2
, **kwargs
):
2594 if J1
.base_ring() != J2
.base_ring():
2595 raise ValueError("algebras must share the same base field")
2596 field
= J1
.base_ring()
2598 self
._factors
= (J1
, J2
)
2602 V
= VectorSpace(field
, n
)
2603 mult_table
= [ [ V
.zero() for j
in range(i
+1) ]
2606 for j
in range(i
+1):
2607 p
= (J1
.monomial(i
)*J1
.monomial(j
)).to_vector()
2608 mult_table
[i
][j
] = V(p
.list() + [field
.zero()]*n2
)
2611 for j
in range(i
+1):
2612 p
= (J2
.monomial(i
)*J2
.monomial(j
)).to_vector()
2613 mult_table
[n1
+i
][n1
+j
] = V([field
.zero()]*n1
+ p
.list())
2615 # TODO: build the IP table here from the two constituent IP
2616 # matrices (it'll be block diagonal, I think).
2617 ip_table
= [ [ field
.zero() for j
in range(i
+1) ]
2619 super(DirectSumEJA
, self
).__init
__(field
,
2624 self
.rank
.set_cache(J1
.rank() + J2
.rank())
2629 Return the pair of this algebra's factors.
2633 sage: from mjo.eja.eja_algebra import (HadamardEJA,
2634 ....: JordanSpinEJA,
2639 sage: J1 = HadamardEJA(2, field=QQ)
2640 sage: J2 = JordanSpinEJA(3, field=QQ)
2641 sage: J = DirectSumEJA(J1,J2)
2643 (Euclidean Jordan algebra of dimension 2 over Rational Field,
2644 Euclidean Jordan algebra of dimension 3 over Rational Field)
2647 return self
._factors
2649 def projections(self
):
2651 Return a pair of projections onto this algebra's factors.
2655 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
2656 ....: ComplexHermitianEJA,
2661 sage: J1 = JordanSpinEJA(2)
2662 sage: J2 = ComplexHermitianEJA(2)
2663 sage: J = DirectSumEJA(J1,J2)
2664 sage: (pi_left, pi_right) = J.projections()
2665 sage: J.one().to_vector()
2667 sage: pi_left(J.one()).to_vector()
2669 sage: pi_right(J.one()).to_vector()
2673 (J1
,J2
) = self
.factors()
2676 V_basis
= self
.vector_space().basis()
2677 # Need to specify the dimensions explicitly so that we don't
2678 # wind up with a zero-by-zero matrix when we want e.g. a
2679 # zero-by-two matrix (important for composing things).
2680 P1
= matrix(self
.base_ring(), m
, m
+n
, V_basis
[:m
])
2681 P2
= matrix(self
.base_ring(), n
, m
+n
, V_basis
[m
:])
2682 pi_left
= FiniteDimensionalEJAOperator(self
,J1
,P1
)
2683 pi_right
= FiniteDimensionalEJAOperator(self
,J2
,P2
)
2684 return (pi_left
, pi_right
)
2686 def inclusions(self
):
2688 Return the pair of inclusion maps from our factors into us.
2692 sage: from mjo.eja.eja_algebra import (random_eja,
2693 ....: JordanSpinEJA,
2694 ....: RealSymmetricEJA,
2699 sage: J1 = JordanSpinEJA(3)
2700 sage: J2 = RealSymmetricEJA(2)
2701 sage: J = DirectSumEJA(J1,J2)
2702 sage: (iota_left, iota_right) = J.inclusions()
2703 sage: iota_left(J1.zero()) == J.zero()
2705 sage: iota_right(J2.zero()) == J.zero()
2707 sage: J1.one().to_vector()
2709 sage: iota_left(J1.one()).to_vector()
2711 sage: J2.one().to_vector()
2713 sage: iota_right(J2.one()).to_vector()
2715 sage: J.one().to_vector()
2720 Composing a projection with the corresponding inclusion should
2721 produce the identity map, and mismatching them should produce
2724 sage: set_random_seed()
2725 sage: J1 = random_eja()
2726 sage: J2 = random_eja()
2727 sage: J = DirectSumEJA(J1,J2)
2728 sage: (iota_left, iota_right) = J.inclusions()
2729 sage: (pi_left, pi_right) = J.projections()
2730 sage: pi_left*iota_left == J1.one().operator()
2732 sage: pi_right*iota_right == J2.one().operator()
2734 sage: (pi_left*iota_right).is_zero()
2736 sage: (pi_right*iota_left).is_zero()
2740 (J1
,J2
) = self
.factors()
2743 V_basis
= self
.vector_space().basis()
2744 # Need to specify the dimensions explicitly so that we don't
2745 # wind up with a zero-by-zero matrix when we want e.g. a
2746 # two-by-zero matrix (important for composing things).
2747 I1
= matrix
.column(self
.base_ring(), m
, m
+n
, V_basis
[:m
])
2748 I2
= matrix
.column(self
.base_ring(), n
, m
+n
, V_basis
[m
:])
2749 iota_left
= FiniteDimensionalEJAOperator(J1
,self
,I1
)
2750 iota_right
= FiniteDimensionalEJAOperator(J2
,self
,I2
)
2751 return (iota_left
, iota_right
)
2753 def inner_product(self
, x
, y
):
2755 The standard Cartesian inner-product.
2757 We project ``x`` and ``y`` onto our factors, and add up the
2758 inner-products from the subalgebras.
2763 sage: from mjo.eja.eja_algebra import (HadamardEJA,
2764 ....: QuaternionHermitianEJA,
2769 sage: J1 = HadamardEJA(3,field=QQ)
2770 sage: J2 = QuaternionHermitianEJA(2,field=QQ,orthonormalize=False)
2771 sage: J = DirectSumEJA(J1,J2)
2776 sage: x1.inner_product(x2)
2778 sage: y1.inner_product(y2)
2780 sage: J.one().inner_product(J.one())
2784 (pi_left
, pi_right
) = self
.projections()
2790 return (x1
.inner_product(y1
) + x2
.inner_product(y2
))
2794 random_eja
= ConcreteEJA
.random_instance