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.
118 # Works on both column and square matrices...
119 degree
= len(basis
[0].list())
121 # Build an ambient space that fits our matrix basis when
122 # written out as "long vectors."
123 V
= VectorSpace(field
, degree
)
125 # The matrix that will hole the orthonormal -> unorthonormal
126 # coordinate transformation.
127 self
._deortho
_matrix
= None
130 # Save a copy of the un-orthonormalized basis for later.
131 # Convert it to ambient V (vector) coordinates while we're
132 # at it, because we'd have to do it later anyway.
133 deortho_vector_basis
= tuple( V(b
.list()) for b
in basis
)
135 from mjo
.eja
.eja_utils
import gram_schmidt
136 basis
= tuple(gram_schmidt(basis
, inner_product
))
138 # Save the (possibly orthonormalized) matrix basis for
140 self
._matrix
_basis
= basis
142 # Now create the vector space for the algebra, which will have
143 # its own set of non-ambient coordinates (in terms of the
145 vector_basis
= tuple( V(b
.list()) for b
in basis
)
146 W
= V
.span_of_basis( vector_basis
, check
=check_axioms
)
149 # Now "W" is the vector space of our algebra coordinates. The
150 # variables "X1", "X2",... refer to the entries of vectors in
151 # W. Thus to convert back and forth between the orthonormal
152 # coordinates and the given ones, we need to stick the original
154 U
= V
.span_of_basis( deortho_vector_basis
, check
=check_axioms
)
155 self
._deortho
_matrix
= matrix( U
.coordinate_vector(q
)
156 for q
in vector_basis
)
159 # Now we actually compute the multiplication and inner-product
160 # tables/matrices using the possibly-orthonormalized basis.
161 self
._inner
_product
_matrix
= matrix
.identity(field
, n
)
162 self
._multiplication
_table
= [ [0 for j
in range(i
+1)]
165 # Note: the Jordan and inner-products are defined in terms
166 # of the ambient basis. It's important that their arguments
167 # are in ambient coordinates as well.
170 # ortho basis w.r.t. ambient coords
174 # The jordan product returns a matrixy answer, so we
175 # have to convert it to the algebra coordinates.
176 elt
= jordan_product(q_i
, q_j
)
177 elt
= W
.coordinate_vector(V(elt
.list()))
178 self
._multiplication
_table
[i
][j
] = self
.from_vector(elt
)
180 if not orthonormalize
:
181 # If we're orthonormalizing the basis with respect
182 # to an inner-product, then the inner-product
183 # matrix with respect to the resulting basis is
184 # just going to be the identity.
185 ip
= inner_product(q_i
, q_j
)
186 self
._inner
_product
_matrix
[i
,j
] = ip
187 self
._inner
_product
_matrix
[j
,i
] = ip
189 self
._inner
_product
_matrix
._cache
= {'hermitian': True}
190 self
._inner
_product
_matrix
.set_immutable()
193 if not self
._is
_jordanian
():
194 raise ValueError("Jordan identity does not hold")
195 if not self
._inner
_product
_is
_associative
():
196 raise ValueError("inner product is not associative")
199 def _coerce_map_from_base_ring(self
):
201 Disable the map from the base ring into the algebra.
203 Performing a nonsense conversion like this automatically
204 is counterpedagogical. The fallback is to try the usual
205 element constructor, which should also fail.
209 sage: from mjo.eja.eja_algebra import random_eja
213 sage: set_random_seed()
214 sage: J = random_eja()
216 Traceback (most recent call last):
218 ValueError: not an element of this algebra
224 def product_on_basis(self
, i
, j
):
225 # We only stored the lower-triangular portion of the
226 # multiplication table.
228 return self
._multiplication
_table
[i
][j
]
230 return self
._multiplication
_table
[j
][i
]
232 def inner_product(self
, x
, y
):
234 The inner product associated with this Euclidean Jordan algebra.
236 Defaults to the trace inner product, but can be overridden by
237 subclasses if they are sure that the necessary properties are
242 sage: from mjo.eja.eja_algebra import (random_eja,
244 ....: BilinearFormEJA)
248 Our inner product is "associative," which means the following for
249 a symmetric bilinear form::
251 sage: set_random_seed()
252 sage: J = random_eja()
253 sage: x,y,z = J.random_elements(3)
254 sage: (x*y).inner_product(z) == y.inner_product(x*z)
259 Ensure that this is the usual inner product for the algebras
262 sage: set_random_seed()
263 sage: J = HadamardEJA.random_instance()
264 sage: x,y = J.random_elements(2)
265 sage: actual = x.inner_product(y)
266 sage: expected = x.to_vector().inner_product(y.to_vector())
267 sage: actual == expected
270 Ensure that this is one-half of the trace inner-product in a
271 BilinearFormEJA that isn't just the reals (when ``n`` isn't
272 one). This is in Faraut and Koranyi, and also my "On the
275 sage: set_random_seed()
276 sage: J = BilinearFormEJA.random_instance()
277 sage: n = J.dimension()
278 sage: x = J.random_element()
279 sage: y = J.random_element()
280 sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
283 B
= self
._inner
_product
_matrix
284 return (B
*x
.to_vector()).inner_product(y
.to_vector())
287 def _is_commutative(self
):
289 Whether or not this algebra's multiplication table is commutative.
291 This method should of course always return ``True``, unless
292 this algebra was constructed with ``check_axioms=False`` and
293 passed an invalid multiplication table.
295 return all( self
.product_on_basis(i
,j
) == self
.product_on_basis(i
,j
)
296 for i
in range(self
.dimension())
297 for j
in range(self
.dimension()) )
299 def _is_jordanian(self
):
301 Whether or not this algebra's multiplication table respects the
302 Jordan identity `(x^{2})(xy) = x(x^{2}y)`.
304 We only check one arrangement of `x` and `y`, so for a
305 ``True`` result to be truly true, you should also check
306 :meth:`_is_commutative`. This method should of course always
307 return ``True``, unless this algebra was constructed with
308 ``check_axioms=False`` and passed an invalid multiplication table.
310 return all( (self
.monomial(i
)**2)*(self
.monomial(i
)*self
.monomial(j
))
312 (self
.monomial(i
))*((self
.monomial(i
)**2)*self
.monomial(j
))
313 for i
in range(self
.dimension())
314 for j
in range(self
.dimension()) )
316 def _inner_product_is_associative(self
):
318 Return whether or not this algebra's inner product `B` is
319 associative; that is, whether or not `B(xy,z) = B(x,yz)`.
321 This method should of course always return ``True``, unless
322 this algebra was constructed with ``check_axioms=False`` and
323 passed an invalid Jordan or inner-product.
326 # Used to check whether or not something is zero in an inexact
327 # ring. This number is sufficient to allow the construction of
328 # QuaternionHermitianEJA(2, field=RDF) with check_axioms=True.
331 for i
in range(self
.dimension()):
332 for j
in range(self
.dimension()):
333 for k
in range(self
.dimension()):
337 diff
= (x
*y
).inner_product(z
) - x
.inner_product(y
*z
)
339 if self
.base_ring().is_exact():
343 if diff
.abs() > epsilon
:
348 def _element_constructor_(self
, elt
):
350 Construct an element of this algebra from its vector or matrix
353 This gets called only after the parent element _call_ method
354 fails to find a coercion for the argument.
358 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
360 ....: RealSymmetricEJA)
364 The identity in `S^n` is converted to the identity in the EJA::
366 sage: J = RealSymmetricEJA(3)
367 sage: I = matrix.identity(QQ,3)
368 sage: J(I) == J.one()
371 This skew-symmetric matrix can't be represented in the EJA::
373 sage: J = RealSymmetricEJA(3)
374 sage: A = matrix(QQ,3, lambda i,j: i-j)
376 Traceback (most recent call last):
378 ValueError: not an element of this algebra
382 Ensure that we can convert any element of the two non-matrix
383 simple algebras (whose matrix representations are columns)
384 back and forth faithfully::
386 sage: set_random_seed()
387 sage: J = HadamardEJA.random_instance()
388 sage: x = J.random_element()
389 sage: J(x.to_vector().column()) == x
391 sage: J = JordanSpinEJA.random_instance()
392 sage: x = J.random_element()
393 sage: J(x.to_vector().column()) == x
397 msg
= "not an element of this algebra"
399 # The superclass implementation of random_element()
400 # needs to be able to coerce "0" into the algebra.
402 elif elt
in self
.base_ring():
403 # Ensure that no base ring -> algebra coercion is performed
404 # by this method. There's some stupidity in sage that would
405 # otherwise propagate to this method; for example, sage thinks
406 # that the integer 3 belongs to the space of 2-by-2 matrices.
407 raise ValueError(msg
)
411 except (AttributeError, TypeError):
412 # Try to convert a vector into a column-matrix
415 if elt
not in self
.matrix_space():
416 raise ValueError(msg
)
418 # Thanks for nothing! Matrix spaces aren't vector spaces in
419 # Sage, so we have to figure out its matrix-basis coordinates
420 # ourselves. We use the basis space's ring instead of the
421 # element's ring because the basis space might be an algebraic
422 # closure whereas the base ring of the 3-by-3 identity matrix
423 # could be QQ instead of QQbar.
425 # We pass check=False because the matrix basis is "guaranteed"
426 # to be linearly independent... right? Ha ha.
427 V
= VectorSpace(self
.base_ring(), elt
.nrows()*elt
.ncols())
428 W
= V
.span_of_basis( (_mat2vec(s
) for s
in self
.matrix_basis()),
432 coords
= W
.coordinate_vector(_mat2vec(elt
))
433 except ArithmeticError: # vector is not in free module
434 raise ValueError(msg
)
436 return self
.from_vector(coords
)
440 Return a string representation of ``self``.
444 sage: from mjo.eja.eja_algebra import JordanSpinEJA
448 Ensure that it says what we think it says::
450 sage: JordanSpinEJA(2, field=AA)
451 Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
452 sage: JordanSpinEJA(3, field=RDF)
453 Euclidean Jordan algebra of dimension 3 over Real Double Field
456 fmt
= "Euclidean Jordan algebra of dimension {} over {}"
457 return fmt
.format(self
.dimension(), self
.base_ring())
461 def characteristic_polynomial_of(self
):
463 Return the algebra's "characteristic polynomial of" function,
464 which is itself a multivariate polynomial that, when evaluated
465 at the coordinates of some algebra element, returns that
466 element's characteristic polynomial.
468 The resulting polynomial has `n+1` variables, where `n` is the
469 dimension of this algebra. The first `n` variables correspond to
470 the coordinates of an algebra element: when evaluated at the
471 coordinates of an algebra element with respect to a certain
472 basis, the result is a univariate polynomial (in the one
473 remaining variable ``t``), namely the characteristic polynomial
478 sage: from mjo.eja.eja_algebra import JordanSpinEJA, TrivialEJA
482 The characteristic polynomial in the spin algebra is given in
483 Alizadeh, Example 11.11::
485 sage: J = JordanSpinEJA(3)
486 sage: p = J.characteristic_polynomial_of(); p
487 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
488 sage: xvec = J.one().to_vector()
492 By definition, the characteristic polynomial is a monic
493 degree-zero polynomial in a rank-zero algebra. Note that
494 Cayley-Hamilton is indeed satisfied since the polynomial
495 ``1`` evaluates to the identity element of the algebra on
498 sage: J = TrivialEJA()
499 sage: J.characteristic_polynomial_of()
506 # The list of coefficient polynomials a_0, a_1, a_2, ..., a_(r-1).
507 a
= self
._charpoly
_coefficients
()
509 # We go to a bit of trouble here to reorder the
510 # indeterminates, so that it's easier to evaluate the
511 # characteristic polynomial at x's coordinates and get back
512 # something in terms of t, which is what we want.
513 S
= PolynomialRing(self
.base_ring(),'t')
517 S
= PolynomialRing(S
, R
.variable_names())
520 return (t
**r
+ sum( a
[k
]*(t
**k
) for k
in range(r
) ))
522 def coordinate_polynomial_ring(self
):
524 The multivariate polynomial ring in which this algebra's
525 :meth:`characteristic_polynomial_of` lives.
529 sage: from mjo.eja.eja_algebra import (HadamardEJA,
530 ....: RealSymmetricEJA)
534 sage: J = HadamardEJA(2)
535 sage: J.coordinate_polynomial_ring()
536 Multivariate Polynomial Ring in X1, X2...
537 sage: J = RealSymmetricEJA(3,field=QQ,orthonormalize=False)
538 sage: J.coordinate_polynomial_ring()
539 Multivariate Polynomial Ring in X1, X2, X3, X4, X5, X6...
542 var_names
= tuple( "X%d" % z
for z
in range(1, self
.dimension()+1) )
543 return PolynomialRing(self
.base_ring(), var_names
)
545 def inner_product(self
, x
, y
):
547 The inner product associated with this Euclidean Jordan algebra.
549 Defaults to the trace inner product, but can be overridden by
550 subclasses if they are sure that the necessary properties are
555 sage: from mjo.eja.eja_algebra import (random_eja,
557 ....: BilinearFormEJA)
561 Our inner product is "associative," which means the following for
562 a symmetric bilinear form::
564 sage: set_random_seed()
565 sage: J = random_eja()
566 sage: x,y,z = J.random_elements(3)
567 sage: (x*y).inner_product(z) == y.inner_product(x*z)
572 Ensure that this is the usual inner product for the algebras
575 sage: set_random_seed()
576 sage: J = HadamardEJA.random_instance()
577 sage: x,y = J.random_elements(2)
578 sage: actual = x.inner_product(y)
579 sage: expected = x.to_vector().inner_product(y.to_vector())
580 sage: actual == expected
583 Ensure that this is one-half of the trace inner-product in a
584 BilinearFormEJA that isn't just the reals (when ``n`` isn't
585 one). This is in Faraut and Koranyi, and also my "On the
588 sage: set_random_seed()
589 sage: J = BilinearFormEJA.random_instance()
590 sage: n = J.dimension()
591 sage: x = J.random_element()
592 sage: y = J.random_element()
593 sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
596 B
= self
._inner
_product
_matrix
597 return (B
*x
.to_vector()).inner_product(y
.to_vector())
600 def is_trivial(self
):
602 Return whether or not this algebra is trivial.
604 A trivial algebra contains only the zero element.
608 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
613 sage: J = ComplexHermitianEJA(3)
619 sage: J = TrivialEJA()
624 return self
.dimension() == 0
627 def multiplication_table(self
):
629 Return a visual representation of this algebra's multiplication
630 table (on basis elements).
634 sage: from mjo.eja.eja_algebra import JordanSpinEJA
638 sage: J = JordanSpinEJA(4)
639 sage: J.multiplication_table()
640 +----++----+----+----+----+
641 | * || e0 | e1 | e2 | e3 |
642 +====++====+====+====+====+
643 | e0 || e0 | e1 | e2 | e3 |
644 +----++----+----+----+----+
645 | e1 || e1 | e0 | 0 | 0 |
646 +----++----+----+----+----+
647 | e2 || e2 | 0 | e0 | 0 |
648 +----++----+----+----+----+
649 | e3 || e3 | 0 | 0 | e0 |
650 +----++----+----+----+----+
654 # Prepend the header row.
655 M
= [["*"] + list(self
.gens())]
657 # And to each subsequent row, prepend an entry that belongs to
658 # the left-side "header column."
659 M
+= [ [self
.monomial(i
)] + [ self
.product_on_basis(i
,j
)
663 return table(M
, header_row
=True, header_column
=True, frame
=True)
666 def matrix_basis(self
):
668 Return an (often more natural) representation of this algebras
669 basis as an ordered tuple of matrices.
671 Every finite-dimensional Euclidean Jordan Algebra is a, up to
672 Jordan isomorphism, a direct sum of five simple
673 algebras---four of which comprise Hermitian matrices. And the
674 last type of algebra can of course be thought of as `n`-by-`1`
675 column matrices (ambiguusly called column vectors) to avoid
676 special cases. As a result, matrices (and column vectors) are
677 a natural representation format for Euclidean Jordan algebra
680 But, when we construct an algebra from a basis of matrices,
681 those matrix representations are lost in favor of coordinate
682 vectors *with respect to* that basis. We could eventually
683 convert back if we tried hard enough, but having the original
684 representations handy is valuable enough that we simply store
685 them and return them from this method.
687 Why implement this for non-matrix algebras? Avoiding special
688 cases for the :class:`BilinearFormEJA` pays with simplicity in
689 its own right. But mainly, we would like to be able to assume
690 that elements of a :class:`DirectSumEJA` can be displayed
691 nicely, without having to have special classes for direct sums
692 one of whose components was a matrix algebra.
696 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
697 ....: RealSymmetricEJA)
701 sage: J = RealSymmetricEJA(2)
703 Finite family {0: e0, 1: e1, 2: e2}
704 sage: J.matrix_basis()
706 [1 0] [ 0 0.7071067811865475?] [0 0]
707 [0 0], [0.7071067811865475? 0], [0 1]
712 sage: J = JordanSpinEJA(2)
714 Finite family {0: e0, 1: e1}
715 sage: J.matrix_basis()
721 return self
._matrix
_basis
724 def matrix_space(self
):
726 Return the matrix space in which this algebra's elements live, if
727 we think of them as matrices (including column vectors of the
730 Generally this will be an `n`-by-`1` column-vector space,
731 except when the algebra is trivial. There it's `n`-by-`n`
732 (where `n` is zero), to ensure that two elements of the matrix
733 space (empty matrices) can be multiplied.
735 Matrix algebras override this with something more useful.
737 if self
.is_trivial():
738 return MatrixSpace(self
.base_ring(), 0)
740 return self
._matrix
_basis
[0].matrix_space()
746 Return the unit element of this algebra.
750 sage: from mjo.eja.eja_algebra import (HadamardEJA,
755 We can compute unit element in the Hadamard EJA::
757 sage: J = HadamardEJA(5)
759 e0 + e1 + e2 + e3 + e4
761 The unit element in the Hadamard EJA is inherited in the
762 subalgebras generated by its elements::
764 sage: J = HadamardEJA(5)
766 e0 + e1 + e2 + e3 + e4
767 sage: x = sum(J.gens())
768 sage: A = x.subalgebra_generated_by(orthonormalize=False)
771 sage: A.one().superalgebra_element()
772 e0 + e1 + e2 + e3 + e4
776 The identity element acts like the identity, regardless of
777 whether or not we orthonormalize::
779 sage: set_random_seed()
780 sage: J = random_eja()
781 sage: x = J.random_element()
782 sage: J.one()*x == x and x*J.one() == x
784 sage: A = x.subalgebra_generated_by()
785 sage: y = A.random_element()
786 sage: A.one()*y == y and y*A.one() == y
791 sage: set_random_seed()
792 sage: J = random_eja(field=QQ, orthonormalize=False)
793 sage: x = J.random_element()
794 sage: J.one()*x == x and x*J.one() == x
796 sage: A = x.subalgebra_generated_by(orthonormalize=False)
797 sage: y = A.random_element()
798 sage: A.one()*y == y and y*A.one() == y
801 The matrix of the unit element's operator is the identity,
802 regardless of the base field and whether or not we
805 sage: set_random_seed()
806 sage: J = random_eja()
807 sage: actual = J.one().operator().matrix()
808 sage: expected = matrix.identity(J.base_ring(), J.dimension())
809 sage: actual == expected
811 sage: x = J.random_element()
812 sage: A = x.subalgebra_generated_by()
813 sage: actual = A.one().operator().matrix()
814 sage: expected = matrix.identity(A.base_ring(), A.dimension())
815 sage: actual == expected
820 sage: set_random_seed()
821 sage: J = random_eja(field=QQ, orthonormalize=False)
822 sage: actual = J.one().operator().matrix()
823 sage: expected = matrix.identity(J.base_ring(), J.dimension())
824 sage: actual == expected
826 sage: x = J.random_element()
827 sage: A = x.subalgebra_generated_by(orthonormalize=False)
828 sage: actual = A.one().operator().matrix()
829 sage: expected = matrix.identity(A.base_ring(), A.dimension())
830 sage: actual == expected
833 Ensure that the cached unit element (often precomputed by
834 hand) agrees with the computed one::
836 sage: set_random_seed()
837 sage: J = random_eja()
838 sage: cached = J.one()
839 sage: J.one.clear_cache()
840 sage: J.one() == cached
845 sage: set_random_seed()
846 sage: J = random_eja(field=QQ, orthonormalize=False)
847 sage: cached = J.one()
848 sage: J.one.clear_cache()
849 sage: J.one() == cached
853 # We can brute-force compute the matrices of the operators
854 # that correspond to the basis elements of this algebra.
855 # If some linear combination of those basis elements is the
856 # algebra identity, then the same linear combination of
857 # their matrices has to be the identity matrix.
859 # Of course, matrices aren't vectors in sage, so we have to
860 # appeal to the "long vectors" isometry.
861 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
863 # Now we use basic linear algebra to find the coefficients,
864 # of the matrices-as-vectors-linear-combination, which should
865 # work for the original algebra basis too.
866 A
= matrix(self
.base_ring(), oper_vecs
)
868 # We used the isometry on the left-hand side already, but we
869 # still need to do it for the right-hand side. Recall that we
870 # wanted something that summed to the identity matrix.
871 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
873 # Now if there's an identity element in the algebra, this
874 # should work. We solve on the left to avoid having to
875 # transpose the matrix "A".
876 return self
.from_vector(A
.solve_left(b
))
879 def peirce_decomposition(self
, c
):
881 The Peirce decomposition of this algebra relative to the
884 In the future, this can be extended to a complete system of
885 orthogonal idempotents.
889 - ``c`` -- an idempotent of this algebra.
893 A triple (J0, J5, J1) containing two subalgebras and one subspace
896 - ``J0`` -- the algebra on the eigenspace of ``c.operator()``
897 corresponding to the eigenvalue zero.
899 - ``J5`` -- the eigenspace (NOT a subalgebra) of ``c.operator()``
900 corresponding to the eigenvalue one-half.
902 - ``J1`` -- the algebra on the eigenspace of ``c.operator()``
903 corresponding to the eigenvalue one.
905 These are the only possible eigenspaces for that operator, and this
906 algebra is a direct sum of them. The spaces ``J0`` and ``J1`` are
907 orthogonal, and are subalgebras of this algebra with the appropriate
912 sage: from mjo.eja.eja_algebra import random_eja, RealSymmetricEJA
916 The canonical example comes from the symmetric matrices, which
917 decompose into diagonal and off-diagonal parts::
919 sage: J = RealSymmetricEJA(3)
920 sage: C = matrix(QQ, [ [1,0,0],
924 sage: J0,J5,J1 = J.peirce_decomposition(c)
926 Euclidean Jordan algebra of dimension 1...
928 Vector space of degree 6 and dimension 2...
930 Euclidean Jordan algebra of dimension 3...
931 sage: J0.one().to_matrix()
935 sage: orig_df = AA.options.display_format
936 sage: AA.options.display_format = 'radical'
937 sage: J.from_vector(J5.basis()[0]).to_matrix()
941 sage: J.from_vector(J5.basis()[1]).to_matrix()
945 sage: AA.options.display_format = orig_df
946 sage: J1.one().to_matrix()
953 Every algebra decomposes trivially with respect to its identity
956 sage: set_random_seed()
957 sage: J = random_eja()
958 sage: J0,J5,J1 = J.peirce_decomposition(J.one())
959 sage: J0.dimension() == 0 and J5.dimension() == 0
961 sage: J1.superalgebra() == J and J1.dimension() == J.dimension()
964 The decomposition is into eigenspaces, and its components are
965 therefore necessarily orthogonal. Moreover, the identity
966 elements in the two subalgebras are the projections onto their
967 respective subspaces of the superalgebra's identity element::
969 sage: set_random_seed()
970 sage: J = random_eja()
971 sage: x = J.random_element()
972 sage: if not J.is_trivial():
973 ....: while x.is_nilpotent():
974 ....: x = J.random_element()
975 sage: c = x.subalgebra_idempotent()
976 sage: J0,J5,J1 = J.peirce_decomposition(c)
978 sage: for (w,y,z) in zip(J0.basis(), J5.basis(), J1.basis()):
979 ....: w = w.superalgebra_element()
980 ....: y = J.from_vector(y)
981 ....: z = z.superalgebra_element()
982 ....: ipsum += w.inner_product(y).abs()
983 ....: ipsum += w.inner_product(z).abs()
984 ....: ipsum += y.inner_product(z).abs()
987 sage: J1(c) == J1.one()
989 sage: J0(J.one() - c) == J0.one()
993 if not c
.is_idempotent():
994 raise ValueError("element is not idempotent: %s" % c
)
996 from mjo
.eja
.eja_subalgebra
import FiniteDimensionalEJASubalgebra
998 # Default these to what they should be if they turn out to be
999 # trivial, because eigenspaces_left() won't return eigenvalues
1000 # corresponding to trivial spaces (e.g. it returns only the
1001 # eigenspace corresponding to lambda=1 if you take the
1002 # decomposition relative to the identity element).
1003 trivial
= FiniteDimensionalEJASubalgebra(self
, ())
1004 J0
= trivial
# eigenvalue zero
1005 J5
= VectorSpace(self
.base_ring(), 0) # eigenvalue one-half
1006 J1
= trivial
# eigenvalue one
1008 for (eigval
, eigspace
) in c
.operator().matrix().right_eigenspaces():
1009 if eigval
== ~
(self
.base_ring()(2)):
1012 gens
= tuple( self
.from_vector(b
) for b
in eigspace
.basis() )
1013 subalg
= FiniteDimensionalEJASubalgebra(self
,
1021 raise ValueError("unexpected eigenvalue: %s" % eigval
)
1026 def random_element(self
, thorough
=False):
1028 Return a random element of this algebra.
1030 Our algebra superclass method only returns a linear
1031 combination of at most two basis elements. We instead
1032 want the vector space "random element" method that
1033 returns a more diverse selection.
1037 - ``thorough`` -- (boolean; default False) whether or not we
1038 should generate irrational coefficients for the random
1039 element when our base ring is irrational; this slows the
1040 algebra operations to a crawl, but any truly random method
1044 # For a general base ring... maybe we can trust this to do the
1045 # right thing? Unlikely, but.
1046 V
= self
.vector_space()
1047 v
= V
.random_element()
1049 if self
.base_ring() is AA
:
1050 # The "random element" method of the algebraic reals is
1051 # stupid at the moment, and only returns integers between
1052 # -2 and 2, inclusive:
1054 # https://trac.sagemath.org/ticket/30875
1056 # Instead, we implement our own "random vector" method,
1057 # and then coerce that into the algebra. We use the vector
1058 # space degree here instead of the dimension because a
1059 # subalgebra could (for example) be spanned by only two
1060 # vectors, each with five coordinates. We need to
1061 # generate all five coordinates.
1063 v
*= QQbar
.random_element().real()
1065 v
*= QQ
.random_element()
1067 return self
.from_vector(V
.coordinate_vector(v
))
1069 def random_elements(self
, count
, thorough
=False):
1071 Return ``count`` random elements as a tuple.
1075 - ``thorough`` -- (boolean; default False) whether or not we
1076 should generate irrational coefficients for the random
1077 elements when our base ring is irrational; this slows the
1078 algebra operations to a crawl, but any truly random method
1083 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1087 sage: J = JordanSpinEJA(3)
1088 sage: x,y,z = J.random_elements(3)
1089 sage: all( [ x in J, y in J, z in J ])
1091 sage: len( J.random_elements(10) ) == 10
1095 return tuple( self
.random_element(thorough
)
1096 for idx
in range(count
) )
1100 def _charpoly_coefficients(self
):
1102 The `r` polynomial coefficients of the "characteristic polynomial
1105 n
= self
.dimension()
1106 R
= self
.coordinate_polynomial_ring()
1108 F
= R
.fraction_field()
1111 # From a result in my book, these are the entries of the
1112 # basis representation of L_x.
1113 return sum( vars[k
]*self
.monomial(k
).operator().matrix()[i
,j
]
1116 L_x
= matrix(F
, n
, n
, L_x_i_j
)
1119 if self
.rank
.is_in_cache():
1121 # There's no need to pad the system with redundant
1122 # columns if we *know* they'll be redundant.
1125 # Compute an extra power in case the rank is equal to
1126 # the dimension (otherwise, we would stop at x^(r-1)).
1127 x_powers
= [ (L_x
**k
)*self
.one().to_vector()
1128 for k
in range(n
+1) ]
1129 A
= matrix
.column(F
, x_powers
[:n
])
1130 AE
= A
.extended_echelon_form()
1137 # The theory says that only the first "r" coefficients are
1138 # nonzero, and they actually live in the original polynomial
1139 # ring and not the fraction field. We negate them because
1140 # in the actual characteristic polynomial, they get moved
1141 # to the other side where x^r lives.
1142 return -A_rref
.solve_right(E
*b
).change_ring(R
)[:r
]
1147 Return the rank of this EJA.
1149 This is a cached method because we know the rank a priori for
1150 all of the algebras we can construct. Thus we can avoid the
1151 expensive ``_charpoly_coefficients()`` call unless we truly
1152 need to compute the whole characteristic polynomial.
1156 sage: from mjo.eja.eja_algebra import (HadamardEJA,
1157 ....: JordanSpinEJA,
1158 ....: RealSymmetricEJA,
1159 ....: ComplexHermitianEJA,
1160 ....: QuaternionHermitianEJA,
1165 The rank of the Jordan spin algebra is always two::
1167 sage: JordanSpinEJA(2).rank()
1169 sage: JordanSpinEJA(3).rank()
1171 sage: JordanSpinEJA(4).rank()
1174 The rank of the `n`-by-`n` Hermitian real, complex, or
1175 quaternion matrices is `n`::
1177 sage: RealSymmetricEJA(4).rank()
1179 sage: ComplexHermitianEJA(3).rank()
1181 sage: QuaternionHermitianEJA(2).rank()
1186 Ensure that every EJA that we know how to construct has a
1187 positive integer rank, unless the algebra is trivial in
1188 which case its rank will be zero::
1190 sage: set_random_seed()
1191 sage: J = random_eja()
1195 sage: r > 0 or (r == 0 and J.is_trivial())
1198 Ensure that computing the rank actually works, since the ranks
1199 of all simple algebras are known and will be cached by default::
1201 sage: set_random_seed() # long time
1202 sage: J = random_eja() # long time
1203 sage: caches = J.rank() # long time
1204 sage: J.rank.clear_cache() # long time
1205 sage: J.rank() == cached # long time
1209 return len(self
._charpoly
_coefficients
())
1212 def vector_space(self
):
1214 Return the vector space that underlies this algebra.
1218 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1222 sage: J = RealSymmetricEJA(2)
1223 sage: J.vector_space()
1224 Vector space of dimension 3 over...
1227 return self
.zero().to_vector().parent().ambient_vector_space()
1230 Element
= FiniteDimensionalEJAElement
1232 class RationalBasisEJA(FiniteDimensionalEJA
):
1234 New class for algebras whose supplied basis elements have all rational entries.
1238 sage: from mjo.eja.eja_algebra import BilinearFormEJA
1242 The supplied basis is orthonormalized by default::
1244 sage: B = matrix(QQ, [[1, 0, 0], [0, 25, -32], [0, -32, 41]])
1245 sage: J = BilinearFormEJA(B)
1246 sage: J.matrix_basis()
1263 # Abuse the check_field parameter to check that the entries of
1264 # out basis (in ambient coordinates) are in the field QQ.
1265 if not all( all(b_i
in QQ
for b_i
in b
.list()) for b
in basis
):
1266 raise TypeError("basis not rational")
1268 self
._rational
_algebra
= None
1270 # There's no point in constructing the extra algebra if this
1271 # one is already rational.
1273 # Note: the same Jordan and inner-products work here,
1274 # because they are necessarily defined with respect to
1275 # ambient coordinates and not any particular basis.
1276 self
._rational
_algebra
= FiniteDimensionalEJA(
1281 orthonormalize
=False,
1285 super().__init
__(basis
,
1289 check_field
=check_field
,
1293 def _charpoly_coefficients(self
):
1297 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
1298 ....: JordanSpinEJA)
1302 The base ring of the resulting polynomial coefficients is what
1303 it should be, and not the rationals (unless the algebra was
1304 already over the rationals)::
1306 sage: J = JordanSpinEJA(3)
1307 sage: J._charpoly_coefficients()
1308 (X1^2 - X2^2 - X3^2, -2*X1)
1309 sage: a0 = J._charpoly_coefficients()[0]
1311 Algebraic Real Field
1312 sage: a0.base_ring()
1313 Algebraic Real Field
1316 if self
._rational
_algebra
is None:
1317 # There's no need to construct *another* algebra over the
1318 # rationals if this one is already over the
1319 # rationals. Likewise, if we never orthonormalized our
1320 # basis, we might as well just use the given one.
1321 return super()._charpoly
_coefficients
()
1323 # Do the computation over the rationals. The answer will be
1324 # the same, because all we've done is a change of basis.
1325 # Then, change back from QQ to our real base ring
1326 a
= ( a_i
.change_ring(self
.base_ring())
1327 for a_i
in self
._rational
_algebra
._charpoly
_coefficients
() )
1329 if self
._deortho
_matrix
is None:
1330 # This can happen if our base ring was, say, AA and we
1331 # chose not to (or didn't need to) orthonormalize. It's
1332 # still faster to do the computations over QQ even if
1333 # the numbers in the boxes stay the same.
1336 # Otherwise, convert the coordinate variables back to the
1337 # deorthonormalized ones.
1338 R
= self
.coordinate_polynomial_ring()
1339 from sage
.modules
.free_module_element
import vector
1340 X
= vector(R
, R
.gens())
1341 BX
= self
._deortho
_matrix
*X
1343 subs_dict
= { X[i]: BX[i] for i in range(len(X)) }
1344 return tuple( a_i
.subs(subs_dict
) for a_i
in a
)
1346 class ConcreteEJA(RationalBasisEJA
):
1348 A class for the Euclidean Jordan algebras that we know by name.
1350 These are the Jordan algebras whose basis, multiplication table,
1351 rank, and so on are known a priori. More to the point, they are
1352 the Euclidean Jordan algebras for which we are able to conjure up
1353 a "random instance."
1357 sage: from mjo.eja.eja_algebra import ConcreteEJA
1361 Our basis is normalized with respect to the algebra's inner
1362 product, unless we specify otherwise::
1364 sage: set_random_seed()
1365 sage: J = ConcreteEJA.random_instance()
1366 sage: all( b.norm() == 1 for b in J.gens() )
1369 Since our basis is orthonormal with respect to the algebra's inner
1370 product, and since we know that this algebra is an EJA, any
1371 left-multiplication operator's matrix will be symmetric because
1372 natural->EJA basis representation is an isometry and within the
1373 EJA the operator is self-adjoint by the Jordan axiom::
1375 sage: set_random_seed()
1376 sage: J = ConcreteEJA.random_instance()
1377 sage: x = J.random_element()
1378 sage: x.operator().is_self_adjoint()
1383 def _max_random_instance_size():
1385 Return an integer "size" that is an upper bound on the size of
1386 this algebra when it is used in a random test
1387 case. Unfortunately, the term "size" is ambiguous -- when
1388 dealing with `R^n` under either the Hadamard or Jordan spin
1389 product, the "size" refers to the dimension `n`. When dealing
1390 with a matrix algebra (real symmetric or complex/quaternion
1391 Hermitian), it refers to the size of the matrix, which is far
1392 less than the dimension of the underlying vector space.
1394 This method must be implemented in each subclass.
1396 raise NotImplementedError
1399 def random_instance(cls
, *args
, **kwargs
):
1401 Return a random instance of this type of algebra.
1403 This method should be implemented in each subclass.
1405 from sage
.misc
.prandom
import choice
1406 eja_class
= choice(cls
.__subclasses
__())
1408 # These all bubble up to the RationalBasisEJA superclass
1409 # constructor, so any (kw)args valid there are also valid
1411 return eja_class
.random_instance(*args
, **kwargs
)
1416 def dimension_over_reals():
1418 The dimension of this matrix's base ring over the reals.
1420 The reals are dimension one over themselves, obviously; that's
1421 just `\mathbb{R}^{1}`. Likewise, the complex numbers `a + bi`
1422 have dimension two. Finally, the quaternions have dimension
1423 four over the reals.
1425 This is used to determine the size of the matrix returned from
1426 :meth:`real_embed`, among other things.
1428 raise NotImplementedError
1431 def real_embed(cls
,M
):
1433 Embed the matrix ``M`` into a space of real matrices.
1435 The matrix ``M`` can have entries in any field at the moment:
1436 the real numbers, complex numbers, or quaternions. And although
1437 they are not a field, we can probably support octonions at some
1438 point, too. This function returns a real matrix that "acts like"
1439 the original with respect to matrix multiplication; i.e.
1441 real_embed(M*N) = real_embed(M)*real_embed(N)
1444 if M
.ncols() != M
.nrows():
1445 raise ValueError("the matrix 'M' must be square")
1450 def real_unembed(cls
,M
):
1452 The inverse of :meth:`real_embed`.
1454 if M
.ncols() != M
.nrows():
1455 raise ValueError("the matrix 'M' must be square")
1456 if not ZZ(M
.nrows()).mod(cls
.dimension_over_reals()).is_zero():
1457 raise ValueError("the matrix 'M' must be a real embedding")
1461 def jordan_product(X
,Y
):
1462 return (X
*Y
+ Y
*X
)/2
1465 def trace_inner_product(cls
,X
,Y
):
1467 Compute the trace inner-product of two real-embeddings.
1471 sage: from mjo.eja.eja_algebra import (RealSymmetricEJA,
1472 ....: ComplexHermitianEJA,
1473 ....: QuaternionHermitianEJA)
1477 This gives the same answer as it would if we computed the trace
1478 from the unembedded (original) matrices::
1480 sage: set_random_seed()
1481 sage: J = RealSymmetricEJA.random_instance()
1482 sage: x,y = J.random_elements(2)
1483 sage: Xe = x.to_matrix()
1484 sage: Ye = y.to_matrix()
1485 sage: X = J.real_unembed(Xe)
1486 sage: Y = J.real_unembed(Ye)
1487 sage: expected = (X*Y).trace()
1488 sage: actual = J.trace_inner_product(Xe,Ye)
1489 sage: actual == expected
1494 sage: set_random_seed()
1495 sage: J = ComplexHermitianEJA.random_instance()
1496 sage: x,y = J.random_elements(2)
1497 sage: Xe = x.to_matrix()
1498 sage: Ye = y.to_matrix()
1499 sage: X = J.real_unembed(Xe)
1500 sage: Y = J.real_unembed(Ye)
1501 sage: expected = (X*Y).trace().real()
1502 sage: actual = J.trace_inner_product(Xe,Ye)
1503 sage: actual == expected
1508 sage: set_random_seed()
1509 sage: J = QuaternionHermitianEJA.random_instance()
1510 sage: x,y = J.random_elements(2)
1511 sage: Xe = x.to_matrix()
1512 sage: Ye = y.to_matrix()
1513 sage: X = J.real_unembed(Xe)
1514 sage: Y = J.real_unembed(Ye)
1515 sage: expected = (X*Y).trace().coefficient_tuple()[0]
1516 sage: actual = J.trace_inner_product(Xe,Ye)
1517 sage: actual == expected
1521 Xu
= cls
.real_unembed(X
)
1522 Yu
= cls
.real_unembed(Y
)
1523 tr
= (Xu
*Yu
).trace()
1526 # Works in QQ, AA, RDF, et cetera.
1528 except AttributeError:
1529 # A quaternion doesn't have a real() method, but does
1530 # have coefficient_tuple() method that returns the
1531 # coefficients of 1, i, j, and k -- in that order.
1532 return tr
.coefficient_tuple()[0]
1535 class RealMatrixEJA(MatrixEJA
):
1537 def dimension_over_reals():
1541 class RealSymmetricEJA(ConcreteEJA
, RealMatrixEJA
):
1543 The rank-n simple EJA consisting of real symmetric n-by-n
1544 matrices, the usual symmetric Jordan product, and the trace inner
1545 product. It has dimension `(n^2 + n)/2` over the reals.
1549 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1553 sage: J = RealSymmetricEJA(2)
1554 sage: e0, e1, e2 = J.gens()
1562 In theory, our "field" can be any subfield of the reals::
1564 sage: RealSymmetricEJA(2, field=RDF)
1565 Euclidean Jordan algebra of dimension 3 over Real Double Field
1566 sage: RealSymmetricEJA(2, field=RR)
1567 Euclidean Jordan algebra of dimension 3 over Real Field with
1568 53 bits of precision
1572 The dimension of this algebra is `(n^2 + n) / 2`::
1574 sage: set_random_seed()
1575 sage: n_max = RealSymmetricEJA._max_random_instance_size()
1576 sage: n = ZZ.random_element(1, n_max)
1577 sage: J = RealSymmetricEJA(n)
1578 sage: J.dimension() == (n^2 + n)/2
1581 The Jordan multiplication is what we think it is::
1583 sage: set_random_seed()
1584 sage: J = RealSymmetricEJA.random_instance()
1585 sage: x,y = J.random_elements(2)
1586 sage: actual = (x*y).to_matrix()
1587 sage: X = x.to_matrix()
1588 sage: Y = y.to_matrix()
1589 sage: expected = (X*Y + Y*X)/2
1590 sage: actual == expected
1592 sage: J(expected) == x*y
1595 We can change the generator prefix::
1597 sage: RealSymmetricEJA(3, prefix='q').gens()
1598 (q0, q1, q2, q3, q4, q5)
1600 We can construct the (trivial) algebra of rank zero::
1602 sage: RealSymmetricEJA(0)
1603 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1607 def _denormalized_basis(cls
, n
):
1609 Return a basis for the space of real symmetric n-by-n matrices.
1613 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1617 sage: set_random_seed()
1618 sage: n = ZZ.random_element(1,5)
1619 sage: B = RealSymmetricEJA._denormalized_basis(n)
1620 sage: all( M.is_symmetric() for M in B)
1624 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1628 for j
in range(i
+1):
1629 Eij
= matrix(ZZ
, n
, lambda k
,l
: k
==i
and l
==j
)
1633 Sij
= Eij
+ Eij
.transpose()
1639 def _max_random_instance_size():
1640 return 4 # Dimension 10
1643 def random_instance(cls
, **kwargs
):
1645 Return a random instance of this type of algebra.
1647 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1648 return cls(n
, **kwargs
)
1650 def __init__(self
, n
, **kwargs
):
1651 # We know this is a valid EJA, but will double-check
1652 # if the user passes check_axioms=True.
1653 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
1655 super(RealSymmetricEJA
, self
).__init
__(self
._denormalized
_basis
(n
),
1656 self
.jordan_product
,
1657 self
.trace_inner_product
,
1660 # TODO: this could be factored out somehow, but is left here
1661 # because the MatrixEJA is not presently a subclass of the
1662 # FDEJA class that defines rank() and one().
1663 self
.rank
.set_cache(n
)
1664 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
1665 self
.one
.set_cache(self(idV
))
1669 class ComplexMatrixEJA(MatrixEJA
):
1670 # A manual dictionary-cache for the complex_extension() method,
1671 # since apparently @classmethods can't also be @cached_methods.
1672 _complex_extension
= {}
1675 def complex_extension(cls
,field
):
1677 The complex field that we embed/unembed, as an extension
1678 of the given ``field``.
1680 if field
in cls
._complex
_extension
:
1681 return cls
._complex
_extension
[field
]
1683 # Sage doesn't know how to adjoin the complex "i" (the root of
1684 # x^2 + 1) to a field in a general way. Here, we just enumerate
1685 # all of the cases that I have cared to support so far.
1687 # Sage doesn't know how to embed AA into QQbar, i.e. how
1688 # to adjoin sqrt(-1) to AA.
1690 elif not field
.is_exact():
1692 F
= field
.complex_field()
1694 # Works for QQ and... maybe some other fields.
1695 R
= PolynomialRing(field
, 'z')
1697 F
= field
.extension(z
**2 + 1, 'I', embedding
=CLF(-1).sqrt())
1699 cls
._complex
_extension
[field
] = F
1703 def dimension_over_reals():
1707 def real_embed(cls
,M
):
1709 Embed the n-by-n complex matrix ``M`` into the space of real
1710 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1711 bi` to the block matrix ``[[a,b],[-b,a]]``.
1715 sage: from mjo.eja.eja_algebra import ComplexMatrixEJA
1719 sage: F = QuadraticField(-1, 'I')
1720 sage: x1 = F(4 - 2*i)
1721 sage: x2 = F(1 + 2*i)
1724 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1725 sage: ComplexMatrixEJA.real_embed(M)
1734 Embedding is a homomorphism (isomorphism, in fact)::
1736 sage: set_random_seed()
1737 sage: n = ZZ.random_element(3)
1738 sage: F = QuadraticField(-1, 'I')
1739 sage: X = random_matrix(F, n)
1740 sage: Y = random_matrix(F, n)
1741 sage: Xe = ComplexMatrixEJA.real_embed(X)
1742 sage: Ye = ComplexMatrixEJA.real_embed(Y)
1743 sage: XYe = ComplexMatrixEJA.real_embed(X*Y)
1748 super(ComplexMatrixEJA
,cls
).real_embed(M
)
1751 # We don't need any adjoined elements...
1752 field
= M
.base_ring().base_ring()
1758 blocks
.append(matrix(field
, 2, [ [ a
, b
],
1761 return matrix
.block(field
, n
, blocks
)
1765 def real_unembed(cls
,M
):
1767 The inverse of _embed_complex_matrix().
1771 sage: from mjo.eja.eja_algebra import ComplexMatrixEJA
1775 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1776 ....: [-2, 1, -4, 3],
1777 ....: [ 9, 10, 11, 12],
1778 ....: [-10, 9, -12, 11] ])
1779 sage: ComplexMatrixEJA.real_unembed(A)
1781 [ 10*I + 9 12*I + 11]
1785 Unembedding is the inverse of embedding::
1787 sage: set_random_seed()
1788 sage: F = QuadraticField(-1, 'I')
1789 sage: M = random_matrix(F, 3)
1790 sage: Me = ComplexMatrixEJA.real_embed(M)
1791 sage: ComplexMatrixEJA.real_unembed(Me) == M
1795 super(ComplexMatrixEJA
,cls
).real_unembed(M
)
1797 d
= cls
.dimension_over_reals()
1798 F
= cls
.complex_extension(M
.base_ring())
1801 # Go top-left to bottom-right (reading order), converting every
1802 # 2-by-2 block we see to a single complex element.
1804 for k
in range(n
/d
):
1805 for j
in range(n
/d
):
1806 submat
= M
[d
*k
:d
*k
+d
,d
*j
:d
*j
+d
]
1807 if submat
[0,0] != submat
[1,1]:
1808 raise ValueError('bad on-diagonal submatrix')
1809 if submat
[0,1] != -submat
[1,0]:
1810 raise ValueError('bad off-diagonal submatrix')
1811 z
= submat
[0,0] + submat
[0,1]*i
1814 return matrix(F
, n
/d
, elements
)
1817 class ComplexHermitianEJA(ConcreteEJA
, ComplexMatrixEJA
):
1819 The rank-n simple EJA consisting of complex Hermitian n-by-n
1820 matrices over the real numbers, the usual symmetric Jordan product,
1821 and the real-part-of-trace inner product. It has dimension `n^2` over
1826 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1830 In theory, our "field" can be any subfield of the reals::
1832 sage: ComplexHermitianEJA(2, field=RDF)
1833 Euclidean Jordan algebra of dimension 4 over Real Double Field
1834 sage: ComplexHermitianEJA(2, field=RR)
1835 Euclidean Jordan algebra of dimension 4 over Real Field with
1836 53 bits of precision
1840 The dimension of this algebra is `n^2`::
1842 sage: set_random_seed()
1843 sage: n_max = ComplexHermitianEJA._max_random_instance_size()
1844 sage: n = ZZ.random_element(1, n_max)
1845 sage: J = ComplexHermitianEJA(n)
1846 sage: J.dimension() == n^2
1849 The Jordan multiplication is what we think it is::
1851 sage: set_random_seed()
1852 sage: J = ComplexHermitianEJA.random_instance()
1853 sage: x,y = J.random_elements(2)
1854 sage: actual = (x*y).to_matrix()
1855 sage: X = x.to_matrix()
1856 sage: Y = y.to_matrix()
1857 sage: expected = (X*Y + Y*X)/2
1858 sage: actual == expected
1860 sage: J(expected) == x*y
1863 We can change the generator prefix::
1865 sage: ComplexHermitianEJA(2, prefix='z').gens()
1868 We can construct the (trivial) algebra of rank zero::
1870 sage: ComplexHermitianEJA(0)
1871 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1876 def _denormalized_basis(cls
, n
):
1878 Returns a basis for the space of complex Hermitian n-by-n matrices.
1880 Why do we embed these? Basically, because all of numerical linear
1881 algebra assumes that you're working with vectors consisting of `n`
1882 entries from a field and scalars from the same field. There's no way
1883 to tell SageMath that (for example) the vectors contain complex
1884 numbers, while the scalar field is real.
1888 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1892 sage: set_random_seed()
1893 sage: n = ZZ.random_element(1,5)
1894 sage: B = ComplexHermitianEJA._denormalized_basis(n)
1895 sage: all( M.is_symmetric() for M in B)
1900 R
= PolynomialRing(field
, 'z')
1902 F
= field
.extension(z
**2 + 1, 'I')
1905 # This is like the symmetric case, but we need to be careful:
1907 # * We want conjugate-symmetry, not just symmetry.
1908 # * The diagonal will (as a result) be real.
1911 Eij
= matrix
.zero(F
,n
)
1913 for j
in range(i
+1):
1917 Sij
= cls
.real_embed(Eij
)
1920 # The second one has a minus because it's conjugated.
1921 Eij
[j
,i
] = 1 # Eij = Eij + Eij.transpose()
1922 Sij_real
= cls
.real_embed(Eij
)
1924 # Eij = I*Eij - I*Eij.transpose()
1927 Sij_imag
= cls
.real_embed(Eij
)
1933 # Since we embedded these, we can drop back to the "field" that we
1934 # started with instead of the complex extension "F".
1935 return tuple( s
.change_ring(field
) for s
in S
)
1938 def __init__(self
, n
, **kwargs
):
1939 # We know this is a valid EJA, but will double-check
1940 # if the user passes check_axioms=True.
1941 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
1943 super(ComplexHermitianEJA
, self
).__init
__(self
._denormalized
_basis
(n
),
1944 self
.jordan_product
,
1945 self
.trace_inner_product
,
1947 # TODO: this could be factored out somehow, but is left here
1948 # because the MatrixEJA is not presently a subclass of the
1949 # FDEJA class that defines rank() and one().
1950 self
.rank
.set_cache(n
)
1951 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
1952 self
.one
.set_cache(self(idV
))
1955 def _max_random_instance_size():
1956 return 3 # Dimension 9
1959 def random_instance(cls
, **kwargs
):
1961 Return a random instance of this type of algebra.
1963 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1964 return cls(n
, **kwargs
)
1966 class QuaternionMatrixEJA(MatrixEJA
):
1968 # A manual dictionary-cache for the quaternion_extension() method,
1969 # since apparently @classmethods can't also be @cached_methods.
1970 _quaternion_extension
= {}
1973 def quaternion_extension(cls
,field
):
1975 The quaternion field that we embed/unembed, as an extension
1976 of the given ``field``.
1978 if field
in cls
._quaternion
_extension
:
1979 return cls
._quaternion
_extension
[field
]
1981 Q
= QuaternionAlgebra(field
,-1,-1)
1983 cls
._quaternion
_extension
[field
] = Q
1987 def dimension_over_reals():
1991 def real_embed(cls
,M
):
1993 Embed the n-by-n quaternion matrix ``M`` into the space of real
1994 matrices of size 4n-by-4n by first sending each quaternion entry `z
1995 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
1996 c+di],[-c + di, a-bi]]`, and then embedding those into a real
2001 sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA
2005 sage: Q = QuaternionAlgebra(QQ,-1,-1)
2006 sage: i,j,k = Q.gens()
2007 sage: x = 1 + 2*i + 3*j + 4*k
2008 sage: M = matrix(Q, 1, [[x]])
2009 sage: QuaternionMatrixEJA.real_embed(M)
2015 Embedding is a homomorphism (isomorphism, in fact)::
2017 sage: set_random_seed()
2018 sage: n = ZZ.random_element(2)
2019 sage: Q = QuaternionAlgebra(QQ,-1,-1)
2020 sage: X = random_matrix(Q, n)
2021 sage: Y = random_matrix(Q, n)
2022 sage: Xe = QuaternionMatrixEJA.real_embed(X)
2023 sage: Ye = QuaternionMatrixEJA.real_embed(Y)
2024 sage: XYe = QuaternionMatrixEJA.real_embed(X*Y)
2029 super(QuaternionMatrixEJA
,cls
).real_embed(M
)
2030 quaternions
= M
.base_ring()
2033 F
= QuadraticField(-1, 'I')
2038 t
= z
.coefficient_tuple()
2043 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
2044 [-c
+ d
*i
, a
- b
*i
]])
2045 realM
= ComplexMatrixEJA
.real_embed(cplxM
)
2046 blocks
.append(realM
)
2048 # We should have real entries by now, so use the realest field
2049 # we've got for the return value.
2050 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
2055 def real_unembed(cls
,M
):
2057 The inverse of _embed_quaternion_matrix().
2061 sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA
2065 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
2066 ....: [-2, 1, -4, 3],
2067 ....: [-3, 4, 1, -2],
2068 ....: [-4, -3, 2, 1]])
2069 sage: QuaternionMatrixEJA.real_unembed(M)
2070 [1 + 2*i + 3*j + 4*k]
2074 Unembedding is the inverse of embedding::
2076 sage: set_random_seed()
2077 sage: Q = QuaternionAlgebra(QQ, -1, -1)
2078 sage: M = random_matrix(Q, 3)
2079 sage: Me = QuaternionMatrixEJA.real_embed(M)
2080 sage: QuaternionMatrixEJA.real_unembed(Me) == M
2084 super(QuaternionMatrixEJA
,cls
).real_unembed(M
)
2086 d
= cls
.dimension_over_reals()
2088 # Use the base ring of the matrix to ensure that its entries can be
2089 # multiplied by elements of the quaternion algebra.
2090 Q
= cls
.quaternion_extension(M
.base_ring())
2093 # Go top-left to bottom-right (reading order), converting every
2094 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
2097 for l
in range(n
/d
):
2098 for m
in range(n
/d
):
2099 submat
= ComplexMatrixEJA
.real_unembed(
2100 M
[d
*l
:d
*l
+d
,d
*m
:d
*m
+d
] )
2101 if submat
[0,0] != submat
[1,1].conjugate():
2102 raise ValueError('bad on-diagonal submatrix')
2103 if submat
[0,1] != -submat
[1,0].conjugate():
2104 raise ValueError('bad off-diagonal submatrix')
2105 z
= submat
[0,0].real()
2106 z
+= submat
[0,0].imag()*i
2107 z
+= submat
[0,1].real()*j
2108 z
+= submat
[0,1].imag()*k
2111 return matrix(Q
, n
/d
, elements
)
2114 class QuaternionHermitianEJA(ConcreteEJA
, QuaternionMatrixEJA
):
2116 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
2117 matrices, the usual symmetric Jordan product, and the
2118 real-part-of-trace inner product. It has dimension `2n^2 - n` over
2123 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2127 In theory, our "field" can be any subfield of the reals::
2129 sage: QuaternionHermitianEJA(2, field=RDF)
2130 Euclidean Jordan algebra of dimension 6 over Real Double Field
2131 sage: QuaternionHermitianEJA(2, field=RR)
2132 Euclidean Jordan algebra of dimension 6 over Real Field with
2133 53 bits of precision
2137 The dimension of this algebra is `2*n^2 - n`::
2139 sage: set_random_seed()
2140 sage: n_max = QuaternionHermitianEJA._max_random_instance_size()
2141 sage: n = ZZ.random_element(1, n_max)
2142 sage: J = QuaternionHermitianEJA(n)
2143 sage: J.dimension() == 2*(n^2) - n
2146 The Jordan multiplication is what we think it is::
2148 sage: set_random_seed()
2149 sage: J = QuaternionHermitianEJA.random_instance()
2150 sage: x,y = J.random_elements(2)
2151 sage: actual = (x*y).to_matrix()
2152 sage: X = x.to_matrix()
2153 sage: Y = y.to_matrix()
2154 sage: expected = (X*Y + Y*X)/2
2155 sage: actual == expected
2157 sage: J(expected) == x*y
2160 We can change the generator prefix::
2162 sage: QuaternionHermitianEJA(2, prefix='a').gens()
2163 (a0, a1, a2, a3, a4, a5)
2165 We can construct the (trivial) algebra of rank zero::
2167 sage: QuaternionHermitianEJA(0)
2168 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2172 def _denormalized_basis(cls
, n
):
2174 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
2176 Why do we embed these? Basically, because all of numerical
2177 linear algebra assumes that you're working with vectors consisting
2178 of `n` entries from a field and scalars from the same field. There's
2179 no way to tell SageMath that (for example) the vectors contain
2180 complex numbers, while the scalar field is real.
2184 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2188 sage: set_random_seed()
2189 sage: n = ZZ.random_element(1,5)
2190 sage: B = QuaternionHermitianEJA._denormalized_basis(n)
2191 sage: all( M.is_symmetric() for M in B )
2196 Q
= QuaternionAlgebra(QQ
,-1,-1)
2199 # This is like the symmetric case, but we need to be careful:
2201 # * We want conjugate-symmetry, not just symmetry.
2202 # * The diagonal will (as a result) be real.
2205 Eij
= matrix
.zero(Q
,n
)
2207 for j
in range(i
+1):
2211 Sij
= cls
.real_embed(Eij
)
2214 # The second, third, and fourth ones have a minus
2215 # because they're conjugated.
2216 # Eij = Eij + Eij.transpose()
2218 Sij_real
= cls
.real_embed(Eij
)
2220 # Eij = I*(Eij - Eij.transpose())
2223 Sij_I
= cls
.real_embed(Eij
)
2225 # Eij = J*(Eij - Eij.transpose())
2228 Sij_J
= cls
.real_embed(Eij
)
2230 # Eij = K*(Eij - Eij.transpose())
2233 Sij_K
= cls
.real_embed(Eij
)
2239 # Since we embedded these, we can drop back to the "field" that we
2240 # started with instead of the quaternion algebra "Q".
2241 return tuple( s
.change_ring(field
) for s
in S
)
2244 def __init__(self
, n
, **kwargs
):
2245 # We know this is a valid EJA, but will double-check
2246 # if the user passes check_axioms=True.
2247 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2249 super(QuaternionHermitianEJA
, self
).__init
__(self
._denormalized
_basis
(n
),
2250 self
.jordan_product
,
2251 self
.trace_inner_product
,
2253 # TODO: this could be factored out somehow, but is left here
2254 # because the MatrixEJA is not presently a subclass of the
2255 # FDEJA class that defines rank() and one().
2256 self
.rank
.set_cache(n
)
2257 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
2258 self
.one
.set_cache(self(idV
))
2262 def _max_random_instance_size():
2264 The maximum rank of a random QuaternionHermitianEJA.
2266 return 2 # Dimension 6
2269 def random_instance(cls
, **kwargs
):
2271 Return a random instance of this type of algebra.
2273 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2274 return cls(n
, **kwargs
)
2277 class HadamardEJA(ConcreteEJA
):
2279 Return the Euclidean Jordan Algebra corresponding to the set
2280 `R^n` under the Hadamard product.
2282 Note: this is nothing more than the Cartesian product of ``n``
2283 copies of the spin algebra. Once Cartesian product algebras
2284 are implemented, this can go.
2288 sage: from mjo.eja.eja_algebra import HadamardEJA
2292 This multiplication table can be verified by hand::
2294 sage: J = HadamardEJA(3)
2295 sage: e0,e1,e2 = J.gens()
2311 We can change the generator prefix::
2313 sage: HadamardEJA(3, prefix='r').gens()
2317 def __init__(self
, n
, **kwargs
):
2319 jordan_product
= lambda x
,y
: x
2320 inner_product
= lambda x
,y
: x
2322 def jordan_product(x
,y
):
2324 return P( xi
*yi
for (xi
,yi
) in zip(x
,y
) )
2326 def inner_product(x
,y
):
2329 # New defaults for keyword arguments. Don't orthonormalize
2330 # because our basis is already orthonormal with respect to our
2331 # inner-product. Don't check the axioms, because we know this
2332 # is a valid EJA... but do double-check if the user passes
2333 # check_axioms=True. Note: we DON'T override the "check_field"
2334 # default here, because the user can pass in a field!
2335 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2336 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2338 column_basis
= tuple( b
.column() for b
in FreeModule(ZZ
, n
).basis() )
2339 super().__init
__(column_basis
, jordan_product
, inner_product
, **kwargs
)
2340 self
.rank
.set_cache(n
)
2343 self
.one
.set_cache( self
.zero() )
2345 self
.one
.set_cache( sum(self
.gens()) )
2348 def _max_random_instance_size():
2350 The maximum dimension of a random HadamardEJA.
2355 def random_instance(cls
, **kwargs
):
2357 Return a random instance of this type of algebra.
2359 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2360 return cls(n
, **kwargs
)
2363 class BilinearFormEJA(ConcreteEJA
):
2365 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2366 with the half-trace inner product and jordan product ``x*y =
2367 (<Bx,y>,y_bar>, x0*y_bar + y0*x_bar)`` where `B = 1 \times B22` is
2368 a symmetric positive-definite "bilinear form" matrix. Its
2369 dimension is the size of `B`, and it has rank two in dimensions
2370 larger than two. It reduces to the ``JordanSpinEJA`` when `B` is
2371 the identity matrix of order ``n``.
2373 We insist that the one-by-one upper-left identity block of `B` be
2374 passed in as well so that we can be passed a matrix of size zero
2375 to construct a trivial algebra.
2379 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
2380 ....: JordanSpinEJA)
2384 When no bilinear form is specified, the identity matrix is used,
2385 and the resulting algebra is the Jordan spin algebra::
2387 sage: B = matrix.identity(AA,3)
2388 sage: J0 = BilinearFormEJA(B)
2389 sage: J1 = JordanSpinEJA(3)
2390 sage: J0.multiplication_table() == J0.multiplication_table()
2393 An error is raised if the matrix `B` does not correspond to a
2394 positive-definite bilinear form::
2396 sage: B = matrix.random(QQ,2,3)
2397 sage: J = BilinearFormEJA(B)
2398 Traceback (most recent call last):
2400 ValueError: bilinear form is not positive-definite
2401 sage: B = matrix.zero(QQ,3)
2402 sage: J = BilinearFormEJA(B)
2403 Traceback (most recent call last):
2405 ValueError: bilinear form is not positive-definite
2409 We can create a zero-dimensional algebra::
2411 sage: B = matrix.identity(AA,0)
2412 sage: J = BilinearFormEJA(B)
2416 We can check the multiplication condition given in the Jordan, von
2417 Neumann, and Wigner paper (and also discussed on my "On the
2418 symmetry..." paper). Note that this relies heavily on the standard
2419 choice of basis, as does anything utilizing the bilinear form
2420 matrix. We opt not to orthonormalize the basis, because if we
2421 did, we would have to normalize the `s_{i}` in a similar manner::
2423 sage: set_random_seed()
2424 sage: n = ZZ.random_element(5)
2425 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2426 sage: B11 = matrix.identity(QQ,1)
2427 sage: B22 = M.transpose()*M
2428 sage: B = block_matrix(2,2,[ [B11,0 ],
2430 sage: J = BilinearFormEJA(B, orthonormalize=False)
2431 sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
2432 sage: V = J.vector_space()
2433 sage: sis = [ J( V([0] + (M.inverse()*ei).list()).column() )
2434 ....: for ei in eis ]
2435 sage: actual = [ sis[i]*sis[j]
2436 ....: for i in range(n-1)
2437 ....: for j in range(n-1) ]
2438 sage: expected = [ J.one() if i == j else J.zero()
2439 ....: for i in range(n-1)
2440 ....: for j in range(n-1) ]
2441 sage: actual == expected
2445 def __init__(self
, B
, **kwargs
):
2446 # The matrix "B" is supplied by the user in most cases,
2447 # so it makes sense to check whether or not its positive-
2448 # definite unless we are specifically asked not to...
2449 if ("check_axioms" not in kwargs
) or kwargs
["check_axioms"]:
2450 if not B
.is_positive_definite():
2451 raise ValueError("bilinear form is not positive-definite")
2453 # However, all of the other data for this EJA is computed
2454 # by us in manner that guarantees the axioms are
2455 # satisfied. So, again, unless we are specifically asked to
2456 # verify things, we'll skip the rest of the checks.
2457 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2459 def inner_product(x
,y
):
2460 return (y
.T
*B
*x
)[0,0]
2462 def jordan_product(x
,y
):
2468 z0
= inner_product(y
,x
)
2469 zbar
= y0
*xbar
+ x0
*ybar
2470 return P([z0
] + zbar
.list())
2473 column_basis
= tuple( b
.column() for b
in FreeModule(ZZ
, n
).basis() )
2474 super(BilinearFormEJA
, self
).__init
__(column_basis
,
2479 # The rank of this algebra is two, unless we're in a
2480 # one-dimensional ambient space (because the rank is bounded
2481 # by the ambient dimension).
2482 self
.rank
.set_cache(min(n
,2))
2485 self
.one
.set_cache( self
.zero() )
2487 self
.one
.set_cache( self
.monomial(0) )
2490 def _max_random_instance_size():
2492 The maximum dimension of a random BilinearFormEJA.
2497 def random_instance(cls
, **kwargs
):
2499 Return a random instance of this algebra.
2501 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2503 B
= matrix
.identity(ZZ
, n
)
2504 return cls(B
, **kwargs
)
2506 B11
= matrix
.identity(ZZ
, 1)
2507 M
= matrix
.random(ZZ
, n
-1)
2508 I
= matrix
.identity(ZZ
, n
-1)
2510 while alpha
.is_zero():
2511 alpha
= ZZ
.random_element().abs()
2512 B22
= M
.transpose()*M
+ alpha
*I
2514 from sage
.matrix
.special
import block_matrix
2515 B
= block_matrix(2,2, [ [B11
, ZZ(0) ],
2518 return cls(B
, **kwargs
)
2521 class JordanSpinEJA(BilinearFormEJA
):
2523 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2524 with the usual inner product and jordan product ``x*y =
2525 (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2530 sage: from mjo.eja.eja_algebra import JordanSpinEJA
2534 This multiplication table can be verified by hand::
2536 sage: J = JordanSpinEJA(4)
2537 sage: e0,e1,e2,e3 = J.gens()
2553 We can change the generator prefix::
2555 sage: JordanSpinEJA(2, prefix='B').gens()
2560 Ensure that we have the usual inner product on `R^n`::
2562 sage: set_random_seed()
2563 sage: J = JordanSpinEJA.random_instance()
2564 sage: x,y = J.random_elements(2)
2565 sage: actual = x.inner_product(y)
2566 sage: expected = x.to_vector().inner_product(y.to_vector())
2567 sage: actual == expected
2571 def __init__(self
, n
, **kwargs
):
2572 # This is a special case of the BilinearFormEJA with the
2573 # identity matrix as its bilinear form.
2574 B
= matrix
.identity(ZZ
, n
)
2576 # Don't orthonormalize because our basis is already
2577 # orthonormal with respect to our inner-product.
2578 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2580 # But also don't pass check_field=False here, because the user
2581 # can pass in a field!
2582 super(JordanSpinEJA
, self
).__init
__(B
, **kwargs
)
2585 def _max_random_instance_size():
2587 The maximum dimension of a random JordanSpinEJA.
2592 def random_instance(cls
, **kwargs
):
2594 Return a random instance of this type of algebra.
2596 Needed here to override the implementation for ``BilinearFormEJA``.
2598 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2599 return cls(n
, **kwargs
)
2602 class TrivialEJA(ConcreteEJA
):
2604 The trivial Euclidean Jordan algebra consisting of only a zero element.
2608 sage: from mjo.eja.eja_algebra import TrivialEJA
2612 sage: J = TrivialEJA()
2619 sage: 7*J.one()*12*J.one()
2621 sage: J.one().inner_product(J.one())
2623 sage: J.one().norm()
2625 sage: J.one().subalgebra_generated_by()
2626 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2631 def __init__(self
, **kwargs
):
2632 jordan_product
= lambda x
,y
: x
2633 inner_product
= lambda x
,y
: 0
2636 # New defaults for keyword arguments
2637 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2638 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2640 super(TrivialEJA
, self
).__init
__(basis
,
2644 # The rank is zero using my definition, namely the dimension of the
2645 # largest subalgebra generated by any element.
2646 self
.rank
.set_cache(0)
2647 self
.one
.set_cache( self
.zero() )
2650 def random_instance(cls
, **kwargs
):
2651 # We don't take a "size" argument so the superclass method is
2652 # inappropriate for us.
2653 return cls(**kwargs
)
2655 # class DirectSumEJA(ConcreteEJA):
2657 # The external (orthogonal) direct sum of two other Euclidean Jordan
2658 # algebras. Essentially the Cartesian product of its two factors.
2659 # Every Euclidean Jordan algebra decomposes into an orthogonal
2660 # direct sum of simple Euclidean Jordan algebras, so no generality
2661 # is lost by providing only this construction.
2665 # sage: from mjo.eja.eja_algebra import (random_eja,
2666 # ....: HadamardEJA,
2667 # ....: RealSymmetricEJA,
2668 # ....: DirectSumEJA)
2672 # sage: J1 = HadamardEJA(2)
2673 # sage: J2 = RealSymmetricEJA(3)
2674 # sage: J = DirectSumEJA(J1,J2)
2675 # sage: J.dimension()
2682 # The external direct sum construction is only valid when the two factors
2683 # have the same base ring; an error is raised otherwise::
2685 # sage: set_random_seed()
2686 # sage: J1 = random_eja(field=AA)
2687 # sage: J2 = random_eja(field=QQ,orthonormalize=False)
2688 # sage: J = DirectSumEJA(J1,J2)
2689 # Traceback (most recent call last):
2691 # ValueError: algebras must share the same base field
2694 # def __init__(self, J1, J2, **kwargs):
2695 # if J1.base_ring() != J2.base_ring():
2696 # raise ValueError("algebras must share the same base field")
2697 # field = J1.base_ring()
2699 # self._factors = (J1, J2)
2700 # n1 = J1.dimension()
2701 # n2 = J2.dimension()
2703 # V = VectorSpace(field, n)
2704 # mult_table = [ [ V.zero() for j in range(i+1) ]
2705 # for i in range(n) ]
2706 # for i in range(n1):
2707 # for j in range(i+1):
2708 # p = (J1.monomial(i)*J1.monomial(j)).to_vector()
2709 # mult_table[i][j] = V(p.list() + [field.zero()]*n2)
2711 # for i in range(n2):
2712 # for j in range(i+1):
2713 # p = (J2.monomial(i)*J2.monomial(j)).to_vector()
2714 # mult_table[n1+i][n1+j] = V([field.zero()]*n1 + p.list())
2716 # # TODO: build the IP table here from the two constituent IP
2717 # # matrices (it'll be block diagonal, I think).
2718 # ip_table = [ [ field.zero() for j in range(i+1) ]
2719 # for i in range(n) ]
2720 # super(DirectSumEJA, self).__init__(field,
2723 # check_axioms=False,
2725 # self.rank.set_cache(J1.rank() + J2.rank())
2728 # def factors(self):
2730 # Return the pair of this algebra's factors.
2734 # sage: from mjo.eja.eja_algebra import (HadamardEJA,
2735 # ....: JordanSpinEJA,
2736 # ....: DirectSumEJA)
2740 # sage: J1 = HadamardEJA(2, field=QQ)
2741 # sage: J2 = JordanSpinEJA(3, field=QQ)
2742 # sage: J = DirectSumEJA(J1,J2)
2744 # (Euclidean Jordan algebra of dimension 2 over Rational Field,
2745 # Euclidean Jordan algebra of dimension 3 over Rational Field)
2748 # return self._factors
2750 # def projections(self):
2752 # Return a pair of projections onto this algebra's factors.
2756 # sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
2757 # ....: ComplexHermitianEJA,
2758 # ....: DirectSumEJA)
2762 # sage: J1 = JordanSpinEJA(2)
2763 # sage: J2 = ComplexHermitianEJA(2)
2764 # sage: J = DirectSumEJA(J1,J2)
2765 # sage: (pi_left, pi_right) = J.projections()
2766 # sage: J.one().to_vector()
2767 # (1, 0, 1, 0, 0, 1)
2768 # sage: pi_left(J.one()).to_vector()
2770 # sage: pi_right(J.one()).to_vector()
2774 # (J1,J2) = self.factors()
2775 # m = J1.dimension()
2776 # n = J2.dimension()
2777 # V_basis = self.vector_space().basis()
2778 # # Need to specify the dimensions explicitly so that we don't
2779 # # wind up with a zero-by-zero matrix when we want e.g. a
2780 # # zero-by-two matrix (important for composing things).
2781 # P1 = matrix(self.base_ring(), m, m+n, V_basis[:m])
2782 # P2 = matrix(self.base_ring(), n, m+n, V_basis[m:])
2783 # pi_left = FiniteDimensionalEJAOperator(self,J1,P1)
2784 # pi_right = FiniteDimensionalEJAOperator(self,J2,P2)
2785 # return (pi_left, pi_right)
2787 # def inclusions(self):
2789 # Return the pair of inclusion maps from our factors into us.
2793 # sage: from mjo.eja.eja_algebra import (random_eja,
2794 # ....: JordanSpinEJA,
2795 # ....: RealSymmetricEJA,
2796 # ....: DirectSumEJA)
2800 # sage: J1 = JordanSpinEJA(3)
2801 # sage: J2 = RealSymmetricEJA(2)
2802 # sage: J = DirectSumEJA(J1,J2)
2803 # sage: (iota_left, iota_right) = J.inclusions()
2804 # sage: iota_left(J1.zero()) == J.zero()
2806 # sage: iota_right(J2.zero()) == J.zero()
2808 # sage: J1.one().to_vector()
2810 # sage: iota_left(J1.one()).to_vector()
2811 # (1, 0, 0, 0, 0, 0)
2812 # sage: J2.one().to_vector()
2814 # sage: iota_right(J2.one()).to_vector()
2815 # (0, 0, 0, 1, 0, 1)
2816 # sage: J.one().to_vector()
2817 # (1, 0, 0, 1, 0, 1)
2821 # Composing a projection with the corresponding inclusion should
2822 # produce the identity map, and mismatching them should produce
2825 # sage: set_random_seed()
2826 # sage: J1 = random_eja()
2827 # sage: J2 = random_eja()
2828 # sage: J = DirectSumEJA(J1,J2)
2829 # sage: (iota_left, iota_right) = J.inclusions()
2830 # sage: (pi_left, pi_right) = J.projections()
2831 # sage: pi_left*iota_left == J1.one().operator()
2833 # sage: pi_right*iota_right == J2.one().operator()
2835 # sage: (pi_left*iota_right).is_zero()
2837 # sage: (pi_right*iota_left).is_zero()
2841 # (J1,J2) = self.factors()
2842 # m = J1.dimension()
2843 # n = J2.dimension()
2844 # V_basis = self.vector_space().basis()
2845 # # Need to specify the dimensions explicitly so that we don't
2846 # # wind up with a zero-by-zero matrix when we want e.g. a
2847 # # two-by-zero matrix (important for composing things).
2848 # I1 = matrix.column(self.base_ring(), m, m+n, V_basis[:m])
2849 # I2 = matrix.column(self.base_ring(), n, m+n, V_basis[m:])
2850 # iota_left = FiniteDimensionalEJAOperator(J1,self,I1)
2851 # iota_right = FiniteDimensionalEJAOperator(J2,self,I2)
2852 # return (iota_left, iota_right)
2854 # def inner_product(self, x, y):
2856 # The standard Cartesian inner-product.
2858 # We project ``x`` and ``y`` onto our factors, and add up the
2859 # inner-products from the subalgebras.
2864 # sage: from mjo.eja.eja_algebra import (HadamardEJA,
2865 # ....: QuaternionHermitianEJA,
2866 # ....: DirectSumEJA)
2870 # sage: J1 = HadamardEJA(3,field=QQ)
2871 # sage: J2 = QuaternionHermitianEJA(2,field=QQ,orthonormalize=False)
2872 # sage: J = DirectSumEJA(J1,J2)
2873 # sage: x1 = J1.one()
2875 # sage: y1 = J2.one()
2877 # sage: x1.inner_product(x2)
2879 # sage: y1.inner_product(y2)
2881 # sage: J.one().inner_product(J.one())
2885 # (pi_left, pi_right) = self.projections()
2891 # return (x1.inner_product(y1) + x2.inner_product(y2))
2895 random_eja
= ConcreteEJA
.random_instance