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
.categories
.sets_cat
import cartesian_product
24 from sage
.combinat
.free_module
import (CombinatorialFreeModule
,
25 CombinatorialFreeModule_CartesianProduct
)
26 from sage
.matrix
.constructor
import matrix
27 from sage
.matrix
.matrix_space
import MatrixSpace
28 from sage
.misc
.cachefunc
import cached_method
29 from sage
.misc
.table
import table
30 from sage
.modules
.free_module
import FreeModule
, VectorSpace
31 from sage
.rings
.all
import (ZZ
, QQ
, AA
, QQbar
, RR
, RLF
, CLF
,
34 from mjo
.eja
.eja_element
import (CartesianProductEJAElement
,
35 FiniteDimensionalEJAElement
)
36 from mjo
.eja
.eja_operator
import FiniteDimensionalEJAOperator
37 from mjo
.eja
.eja_utils
import _mat2vec
39 class FiniteDimensionalEJA(CombinatorialFreeModule
):
41 A finite-dimensional Euclidean Jordan algebra.
45 - basis -- a tuple of basis elements in "matrix form," which
46 must be the same form as the arguments to ``jordan_product``
47 and ``inner_product``. In reality, "matrix form" can be either
48 vectors, matrices, or a Cartesian product (ordered tuple)
49 of vectors or matrices. All of these would ideally be vector
50 spaces in sage with no special-casing needed; but in reality
51 we turn vectors into column-matrices and Cartesian products
52 `(a,b)` into column matrices `(a,b)^{T}` after converting
53 `a` and `b` themselves.
55 - jordan_product -- function of two elements (in matrix form)
56 that returns their jordan product in this algebra; this will
57 be applied to ``basis`` to compute a multiplication table for
60 - inner_product -- function of two elements (in matrix form) that
61 returns their inner product. This will be applied to ``basis`` to
62 compute an inner-product table (basically a matrix) for this algebra.
65 Element
= FiniteDimensionalEJAElement
74 cartesian_product
=False,
80 if not field
.is_subring(RR
):
81 # Note: this does return true for the real algebraic
82 # field, the rationals, and any quadratic field where
83 # we've specified a real embedding.
84 raise ValueError("scalar field is not real")
86 # If the basis given to us wasn't over the field that it's
87 # supposed to be over, fix that. Or, you know, crash.
88 basis
= tuple( b
.change_ring(field
) for b
in basis
)
91 # Check commutativity of the Jordan and inner-products.
92 # This has to be done before we build the multiplication
93 # and inner-product tables/matrices, because we take
94 # advantage of symmetry in the process.
95 if not all( jordan_product(bi
,bj
) == jordan_product(bj
,bi
)
98 raise ValueError("Jordan product is not commutative")
100 if not all( inner_product(bi
,bj
) == inner_product(bj
,bi
)
103 raise ValueError("inner-product is not commutative")
106 category
= MagmaticAlgebras(field
).FiniteDimensional()
107 category
= category
.WithBasis().Unital()
109 # Element subalgebras can take advantage of this.
110 category
= category
.Associative()
111 if cartesian_product
:
112 category
= category
.CartesianProducts()
114 # Call the superclass constructor so that we can use its from_vector()
115 # method to build our multiplication table.
117 super().__init
__(field
,
123 # Now comes all of the hard work. We'll be constructing an
124 # ambient vector space V that our (vectorized) basis lives in,
125 # as well as a subspace W of V spanned by those (vectorized)
126 # basis elements. The W-coordinates are the coefficients that
127 # we see in things like x = 1*e1 + 2*e2.
132 # Works on both column and square matrices...
133 degree
= len(basis
[0].list())
135 # Build an ambient space that fits our matrix basis when
136 # written out as "long vectors."
137 V
= VectorSpace(field
, degree
)
139 # The matrix that will hole the orthonormal -> unorthonormal
140 # coordinate transformation.
141 self
._deortho
_matrix
= None
144 # Save a copy of the un-orthonormalized basis for later.
145 # Convert it to ambient V (vector) coordinates while we're
146 # at it, because we'd have to do it later anyway.
147 deortho_vector_basis
= tuple( V(b
.list()) for b
in basis
)
149 from mjo
.eja
.eja_utils
import gram_schmidt
150 basis
= tuple(gram_schmidt(basis
, inner_product
))
152 # Save the (possibly orthonormalized) matrix basis for
154 self
._matrix
_basis
= basis
156 # Now create the vector space for the algebra, which will have
157 # its own set of non-ambient coordinates (in terms of the
159 vector_basis
= tuple( V(b
.list()) for b
in basis
)
160 W
= V
.span_of_basis( vector_basis
, check
=check_axioms
)
163 # Now "W" is the vector space of our algebra coordinates. The
164 # variables "X1", "X2",... refer to the entries of vectors in
165 # W. Thus to convert back and forth between the orthonormal
166 # coordinates and the given ones, we need to stick the original
168 U
= V
.span_of_basis( deortho_vector_basis
, check
=check_axioms
)
169 self
._deortho
_matrix
= matrix( U
.coordinate_vector(q
)
170 for q
in vector_basis
)
173 # Now we actually compute the multiplication and inner-product
174 # tables/matrices using the possibly-orthonormalized basis.
175 self
._inner
_product
_matrix
= matrix
.identity(field
, n
)
176 self
._multiplication
_table
= [ [0 for j
in range(i
+1)]
179 # Note: the Jordan and inner-products are defined in terms
180 # of the ambient basis. It's important that their arguments
181 # are in ambient coordinates as well.
184 # ortho basis w.r.t. ambient coords
188 # The jordan product returns a matrixy answer, so we
189 # have to convert it to the algebra coordinates.
190 elt
= jordan_product(q_i
, q_j
)
191 elt
= W
.coordinate_vector(V(elt
.list()))
192 self
._multiplication
_table
[i
][j
] = self
.from_vector(elt
)
194 if not orthonormalize
:
195 # If we're orthonormalizing the basis with respect
196 # to an inner-product, then the inner-product
197 # matrix with respect to the resulting basis is
198 # just going to be the identity.
199 ip
= inner_product(q_i
, q_j
)
200 self
._inner
_product
_matrix
[i
,j
] = ip
201 self
._inner
_product
_matrix
[j
,i
] = ip
203 self
._inner
_product
_matrix
._cache
= {'hermitian': True}
204 self
._inner
_product
_matrix
.set_immutable()
207 if not self
._is
_jordanian
():
208 raise ValueError("Jordan identity does not hold")
209 if not self
._inner
_product
_is
_associative
():
210 raise ValueError("inner product is not associative")
213 def _coerce_map_from_base_ring(self
):
215 Disable the map from the base ring into the algebra.
217 Performing a nonsense conversion like this automatically
218 is counterpedagogical. The fallback is to try the usual
219 element constructor, which should also fail.
223 sage: from mjo.eja.eja_algebra import random_eja
227 sage: set_random_seed()
228 sage: J = random_eja()
230 Traceback (most recent call last):
232 ValueError: not an element of this algebra
238 def product_on_basis(self
, i
, j
):
239 # We only stored the lower-triangular portion of the
240 # multiplication table.
242 return self
._multiplication
_table
[i
][j
]
244 return self
._multiplication
_table
[j
][i
]
246 def inner_product(self
, x
, y
):
248 The inner product associated with this Euclidean Jordan algebra.
250 Defaults to the trace inner product, but can be overridden by
251 subclasses if they are sure that the necessary properties are
256 sage: from mjo.eja.eja_algebra import (random_eja,
258 ....: BilinearFormEJA)
262 Our inner product is "associative," which means the following for
263 a symmetric bilinear form::
265 sage: set_random_seed()
266 sage: J = random_eja()
267 sage: x,y,z = J.random_elements(3)
268 sage: (x*y).inner_product(z) == y.inner_product(x*z)
273 Ensure that this is the usual inner product for the algebras
276 sage: set_random_seed()
277 sage: J = HadamardEJA.random_instance()
278 sage: x,y = J.random_elements(2)
279 sage: actual = x.inner_product(y)
280 sage: expected = x.to_vector().inner_product(y.to_vector())
281 sage: actual == expected
284 Ensure that this is one-half of the trace inner-product in a
285 BilinearFormEJA that isn't just the reals (when ``n`` isn't
286 one). This is in Faraut and Koranyi, and also my "On the
289 sage: set_random_seed()
290 sage: J = BilinearFormEJA.random_instance()
291 sage: n = J.dimension()
292 sage: x = J.random_element()
293 sage: y = J.random_element()
294 sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
297 B
= self
._inner
_product
_matrix
298 return (B
*x
.to_vector()).inner_product(y
.to_vector())
301 def _is_commutative(self
):
303 Whether or not this algebra's multiplication table is commutative.
305 This method should of course always return ``True``, unless
306 this algebra was constructed with ``check_axioms=False`` and
307 passed an invalid multiplication table.
309 return all( self
.product_on_basis(i
,j
) == self
.product_on_basis(i
,j
)
310 for i
in range(self
.dimension())
311 for j
in range(self
.dimension()) )
313 def _is_jordanian(self
):
315 Whether or not this algebra's multiplication table respects the
316 Jordan identity `(x^{2})(xy) = x(x^{2}y)`.
318 We only check one arrangement of `x` and `y`, so for a
319 ``True`` result to be truly true, you should also check
320 :meth:`_is_commutative`. This method should of course always
321 return ``True``, unless this algebra was constructed with
322 ``check_axioms=False`` and passed an invalid multiplication table.
324 return all( (self
.gens()[i
]**2)*(self
.gens()[i
]*self
.gens()[j
])
326 (self
.gens()[i
])*((self
.gens()[i
]**2)*self
.gens()[j
])
327 for i
in range(self
.dimension())
328 for j
in range(self
.dimension()) )
330 def _inner_product_is_associative(self
):
332 Return whether or not this algebra's inner product `B` is
333 associative; that is, whether or not `B(xy,z) = B(x,yz)`.
335 This method should of course always return ``True``, unless
336 this algebra was constructed with ``check_axioms=False`` and
337 passed an invalid Jordan or inner-product.
340 # Used to check whether or not something is zero in an inexact
341 # ring. This number is sufficient to allow the construction of
342 # QuaternionHermitianEJA(2, field=RDF) with check_axioms=True.
345 for i
in range(self
.dimension()):
346 for j
in range(self
.dimension()):
347 for k
in range(self
.dimension()):
351 diff
= (x
*y
).inner_product(z
) - x
.inner_product(y
*z
)
353 if self
.base_ring().is_exact():
357 if diff
.abs() > epsilon
:
362 def _element_constructor_(self
, elt
):
364 Construct an element of this algebra from its vector or matrix
367 This gets called only after the parent element _call_ method
368 fails to find a coercion for the argument.
372 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
374 ....: RealSymmetricEJA)
378 The identity in `S^n` is converted to the identity in the EJA::
380 sage: J = RealSymmetricEJA(3)
381 sage: I = matrix.identity(QQ,3)
382 sage: J(I) == J.one()
385 This skew-symmetric matrix can't be represented in the EJA::
387 sage: J = RealSymmetricEJA(3)
388 sage: A = matrix(QQ,3, lambda i,j: i-j)
390 Traceback (most recent call last):
392 ValueError: not an element of this algebra
396 Ensure that we can convert any element of the two non-matrix
397 simple algebras (whose matrix representations are columns)
398 back and forth faithfully::
400 sage: set_random_seed()
401 sage: J = HadamardEJA.random_instance()
402 sage: x = J.random_element()
403 sage: J(x.to_vector().column()) == x
405 sage: J = JordanSpinEJA.random_instance()
406 sage: x = J.random_element()
407 sage: J(x.to_vector().column()) == x
411 msg
= "not an element of this algebra"
413 # The superclass implementation of random_element()
414 # needs to be able to coerce "0" into the algebra.
416 elif elt
in self
.base_ring():
417 # Ensure that no base ring -> algebra coercion is performed
418 # by this method. There's some stupidity in sage that would
419 # otherwise propagate to this method; for example, sage thinks
420 # that the integer 3 belongs to the space of 2-by-2 matrices.
421 raise ValueError(msg
)
425 except (AttributeError, TypeError):
426 # Try to convert a vector into a column-matrix
429 if elt
not in self
.matrix_space():
430 raise ValueError(msg
)
432 # Thanks for nothing! Matrix spaces aren't vector spaces in
433 # Sage, so we have to figure out its matrix-basis coordinates
434 # ourselves. We use the basis space's ring instead of the
435 # element's ring because the basis space might be an algebraic
436 # closure whereas the base ring of the 3-by-3 identity matrix
437 # could be QQ instead of QQbar.
439 # We pass check=False because the matrix basis is "guaranteed"
440 # to be linearly independent... right? Ha ha.
441 V
= VectorSpace(self
.base_ring(), elt
.nrows()*elt
.ncols())
442 W
= V
.span_of_basis( (_mat2vec(s
) for s
in self
.matrix_basis()),
446 coords
= W
.coordinate_vector(_mat2vec(elt
))
447 except ArithmeticError: # vector is not in free module
448 raise ValueError(msg
)
450 return self
.from_vector(coords
)
454 Return a string representation of ``self``.
458 sage: from mjo.eja.eja_algebra import JordanSpinEJA
462 Ensure that it says what we think it says::
464 sage: JordanSpinEJA(2, field=AA)
465 Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
466 sage: JordanSpinEJA(3, field=RDF)
467 Euclidean Jordan algebra of dimension 3 over Real Double Field
470 fmt
= "Euclidean Jordan algebra of dimension {} over {}"
471 return fmt
.format(self
.dimension(), self
.base_ring())
475 def characteristic_polynomial_of(self
):
477 Return the algebra's "characteristic polynomial of" function,
478 which is itself a multivariate polynomial that, when evaluated
479 at the coordinates of some algebra element, returns that
480 element's characteristic polynomial.
482 The resulting polynomial has `n+1` variables, where `n` is the
483 dimension of this algebra. The first `n` variables correspond to
484 the coordinates of an algebra element: when evaluated at the
485 coordinates of an algebra element with respect to a certain
486 basis, the result is a univariate polynomial (in the one
487 remaining variable ``t``), namely the characteristic polynomial
492 sage: from mjo.eja.eja_algebra import JordanSpinEJA, TrivialEJA
496 The characteristic polynomial in the spin algebra is given in
497 Alizadeh, Example 11.11::
499 sage: J = JordanSpinEJA(3)
500 sage: p = J.characteristic_polynomial_of(); p
501 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
502 sage: xvec = J.one().to_vector()
506 By definition, the characteristic polynomial is a monic
507 degree-zero polynomial in a rank-zero algebra. Note that
508 Cayley-Hamilton is indeed satisfied since the polynomial
509 ``1`` evaluates to the identity element of the algebra on
512 sage: J = TrivialEJA()
513 sage: J.characteristic_polynomial_of()
520 # The list of coefficient polynomials a_0, a_1, a_2, ..., a_(r-1).
521 a
= self
._charpoly
_coefficients
()
523 # We go to a bit of trouble here to reorder the
524 # indeterminates, so that it's easier to evaluate the
525 # characteristic polynomial at x's coordinates and get back
526 # something in terms of t, which is what we want.
527 S
= PolynomialRing(self
.base_ring(),'t')
531 S
= PolynomialRing(S
, R
.variable_names())
534 return (t
**r
+ sum( a
[k
]*(t
**k
) for k
in range(r
) ))
536 def coordinate_polynomial_ring(self
):
538 The multivariate polynomial ring in which this algebra's
539 :meth:`characteristic_polynomial_of` lives.
543 sage: from mjo.eja.eja_algebra import (HadamardEJA,
544 ....: RealSymmetricEJA)
548 sage: J = HadamardEJA(2)
549 sage: J.coordinate_polynomial_ring()
550 Multivariate Polynomial Ring in X1, X2...
551 sage: J = RealSymmetricEJA(3,field=QQ,orthonormalize=False)
552 sage: J.coordinate_polynomial_ring()
553 Multivariate Polynomial Ring in X1, X2, X3, X4, X5, X6...
556 var_names
= tuple( "X%d" % z
for z
in range(1, self
.dimension()+1) )
557 return PolynomialRing(self
.base_ring(), var_names
)
559 def inner_product(self
, x
, y
):
561 The inner product associated with this Euclidean Jordan algebra.
563 Defaults to the trace inner product, but can be overridden by
564 subclasses if they are sure that the necessary properties are
569 sage: from mjo.eja.eja_algebra import (random_eja,
571 ....: BilinearFormEJA)
575 Our inner product is "associative," which means the following for
576 a symmetric bilinear form::
578 sage: set_random_seed()
579 sage: J = random_eja()
580 sage: x,y,z = J.random_elements(3)
581 sage: (x*y).inner_product(z) == y.inner_product(x*z)
586 Ensure that this is the usual inner product for the algebras
589 sage: set_random_seed()
590 sage: J = HadamardEJA.random_instance()
591 sage: x,y = J.random_elements(2)
592 sage: actual = x.inner_product(y)
593 sage: expected = x.to_vector().inner_product(y.to_vector())
594 sage: actual == expected
597 Ensure that this is one-half of the trace inner-product in a
598 BilinearFormEJA that isn't just the reals (when ``n`` isn't
599 one). This is in Faraut and Koranyi, and also my "On the
602 sage: set_random_seed()
603 sage: J = BilinearFormEJA.random_instance()
604 sage: n = J.dimension()
605 sage: x = J.random_element()
606 sage: y = J.random_element()
607 sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
610 B
= self
._inner
_product
_matrix
611 return (B
*x
.to_vector()).inner_product(y
.to_vector())
614 def is_trivial(self
):
616 Return whether or not this algebra is trivial.
618 A trivial algebra contains only the zero element.
622 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
627 sage: J = ComplexHermitianEJA(3)
633 sage: J = TrivialEJA()
638 return self
.dimension() == 0
641 def multiplication_table(self
):
643 Return a visual representation of this algebra's multiplication
644 table (on basis elements).
648 sage: from mjo.eja.eja_algebra import JordanSpinEJA
652 sage: J = JordanSpinEJA(4)
653 sage: J.multiplication_table()
654 +----++----+----+----+----+
655 | * || e0 | e1 | e2 | e3 |
656 +====++====+====+====+====+
657 | e0 || e0 | e1 | e2 | e3 |
658 +----++----+----+----+----+
659 | e1 || e1 | e0 | 0 | 0 |
660 +----++----+----+----+----+
661 | e2 || e2 | 0 | e0 | 0 |
662 +----++----+----+----+----+
663 | e3 || e3 | 0 | 0 | e0 |
664 +----++----+----+----+----+
668 # Prepend the header row.
669 M
= [["*"] + list(self
.gens())]
671 # And to each subsequent row, prepend an entry that belongs to
672 # the left-side "header column."
673 M
+= [ [self
.gens()[i
]] + [ self
.product_on_basis(i
,j
)
677 return table(M
, header_row
=True, header_column
=True, frame
=True)
680 def matrix_basis(self
):
682 Return an (often more natural) representation of this algebras
683 basis as an ordered tuple of matrices.
685 Every finite-dimensional Euclidean Jordan Algebra is a, up to
686 Jordan isomorphism, a direct sum of five simple
687 algebras---four of which comprise Hermitian matrices. And the
688 last type of algebra can of course be thought of as `n`-by-`1`
689 column matrices (ambiguusly called column vectors) to avoid
690 special cases. As a result, matrices (and column vectors) are
691 a natural representation format for Euclidean Jordan algebra
694 But, when we construct an algebra from a basis of matrices,
695 those matrix representations are lost in favor of coordinate
696 vectors *with respect to* that basis. We could eventually
697 convert back if we tried hard enough, but having the original
698 representations handy is valuable enough that we simply store
699 them and return them from this method.
701 Why implement this for non-matrix algebras? Avoiding special
702 cases for the :class:`BilinearFormEJA` pays with simplicity in
703 its own right. But mainly, we would like to be able to assume
704 that elements of a :class:`CartesianProductEJA` can be displayed
705 nicely, without having to have special classes for direct sums
706 one of whose components was a matrix algebra.
710 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
711 ....: RealSymmetricEJA)
715 sage: J = RealSymmetricEJA(2)
717 Finite family {0: e0, 1: e1, 2: e2}
718 sage: J.matrix_basis()
720 [1 0] [ 0 0.7071067811865475?] [0 0]
721 [0 0], [0.7071067811865475? 0], [0 1]
726 sage: J = JordanSpinEJA(2)
728 Finite family {0: e0, 1: e1}
729 sage: J.matrix_basis()
735 return self
._matrix
_basis
738 def matrix_space(self
):
740 Return the matrix space in which this algebra's elements live, if
741 we think of them as matrices (including column vectors of the
744 Generally this will be an `n`-by-`1` column-vector space,
745 except when the algebra is trivial. There it's `n`-by-`n`
746 (where `n` is zero), to ensure that two elements of the matrix
747 space (empty matrices) can be multiplied.
749 Matrix algebras override this with something more useful.
751 if self
.is_trivial():
752 return MatrixSpace(self
.base_ring(), 0)
754 return self
.matrix_basis()[0].parent()
760 Return the unit element of this algebra.
764 sage: from mjo.eja.eja_algebra import (HadamardEJA,
769 We can compute unit element in the Hadamard EJA::
771 sage: J = HadamardEJA(5)
773 e0 + e1 + e2 + e3 + e4
775 The unit element in the Hadamard EJA is inherited in the
776 subalgebras generated by its elements::
778 sage: J = HadamardEJA(5)
780 e0 + e1 + e2 + e3 + e4
781 sage: x = sum(J.gens())
782 sage: A = x.subalgebra_generated_by(orthonormalize=False)
785 sage: A.one().superalgebra_element()
786 e0 + e1 + e2 + e3 + e4
790 The identity element acts like the identity, regardless of
791 whether or not we orthonormalize::
793 sage: set_random_seed()
794 sage: J = random_eja()
795 sage: x = J.random_element()
796 sage: J.one()*x == x and x*J.one() == x
798 sage: A = x.subalgebra_generated_by()
799 sage: y = A.random_element()
800 sage: A.one()*y == y and y*A.one() == y
805 sage: set_random_seed()
806 sage: J = random_eja(field=QQ, orthonormalize=False)
807 sage: x = J.random_element()
808 sage: J.one()*x == x and x*J.one() == x
810 sage: A = x.subalgebra_generated_by(orthonormalize=False)
811 sage: y = A.random_element()
812 sage: A.one()*y == y and y*A.one() == y
815 The matrix of the unit element's operator is the identity,
816 regardless of the base field and whether or not we
819 sage: set_random_seed()
820 sage: J = random_eja()
821 sage: actual = J.one().operator().matrix()
822 sage: expected = matrix.identity(J.base_ring(), J.dimension())
823 sage: actual == expected
825 sage: x = J.random_element()
826 sage: A = x.subalgebra_generated_by()
827 sage: actual = A.one().operator().matrix()
828 sage: expected = matrix.identity(A.base_ring(), A.dimension())
829 sage: actual == expected
834 sage: set_random_seed()
835 sage: J = random_eja(field=QQ, orthonormalize=False)
836 sage: actual = J.one().operator().matrix()
837 sage: expected = matrix.identity(J.base_ring(), J.dimension())
838 sage: actual == expected
840 sage: x = J.random_element()
841 sage: A = x.subalgebra_generated_by(orthonormalize=False)
842 sage: actual = A.one().operator().matrix()
843 sage: expected = matrix.identity(A.base_ring(), A.dimension())
844 sage: actual == expected
847 Ensure that the cached unit element (often precomputed by
848 hand) agrees with the computed one::
850 sage: set_random_seed()
851 sage: J = random_eja()
852 sage: cached = J.one()
853 sage: J.one.clear_cache()
854 sage: J.one() == cached
859 sage: set_random_seed()
860 sage: J = random_eja(field=QQ, orthonormalize=False)
861 sage: cached = J.one()
862 sage: J.one.clear_cache()
863 sage: J.one() == cached
867 # We can brute-force compute the matrices of the operators
868 # that correspond to the basis elements of this algebra.
869 # If some linear combination of those basis elements is the
870 # algebra identity, then the same linear combination of
871 # their matrices has to be the identity matrix.
873 # Of course, matrices aren't vectors in sage, so we have to
874 # appeal to the "long vectors" isometry.
875 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
877 # Now we use basic linear algebra to find the coefficients,
878 # of the matrices-as-vectors-linear-combination, which should
879 # work for the original algebra basis too.
880 A
= matrix(self
.base_ring(), oper_vecs
)
882 # We used the isometry on the left-hand side already, but we
883 # still need to do it for the right-hand side. Recall that we
884 # wanted something that summed to the identity matrix.
885 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
887 # Now if there's an identity element in the algebra, this
888 # should work. We solve on the left to avoid having to
889 # transpose the matrix "A".
890 return self
.from_vector(A
.solve_left(b
))
893 def peirce_decomposition(self
, c
):
895 The Peirce decomposition of this algebra relative to the
898 In the future, this can be extended to a complete system of
899 orthogonal idempotents.
903 - ``c`` -- an idempotent of this algebra.
907 A triple (J0, J5, J1) containing two subalgebras and one subspace
910 - ``J0`` -- the algebra on the eigenspace of ``c.operator()``
911 corresponding to the eigenvalue zero.
913 - ``J5`` -- the eigenspace (NOT a subalgebra) of ``c.operator()``
914 corresponding to the eigenvalue one-half.
916 - ``J1`` -- the algebra on the eigenspace of ``c.operator()``
917 corresponding to the eigenvalue one.
919 These are the only possible eigenspaces for that operator, and this
920 algebra is a direct sum of them. The spaces ``J0`` and ``J1`` are
921 orthogonal, and are subalgebras of this algebra with the appropriate
926 sage: from mjo.eja.eja_algebra import random_eja, RealSymmetricEJA
930 The canonical example comes from the symmetric matrices, which
931 decompose into diagonal and off-diagonal parts::
933 sage: J = RealSymmetricEJA(3)
934 sage: C = matrix(QQ, [ [1,0,0],
938 sage: J0,J5,J1 = J.peirce_decomposition(c)
940 Euclidean Jordan algebra of dimension 1...
942 Vector space of degree 6 and dimension 2...
944 Euclidean Jordan algebra of dimension 3...
945 sage: J0.one().to_matrix()
949 sage: orig_df = AA.options.display_format
950 sage: AA.options.display_format = 'radical'
951 sage: J.from_vector(J5.basis()[0]).to_matrix()
955 sage: J.from_vector(J5.basis()[1]).to_matrix()
959 sage: AA.options.display_format = orig_df
960 sage: J1.one().to_matrix()
967 Every algebra decomposes trivially with respect to its identity
970 sage: set_random_seed()
971 sage: J = random_eja()
972 sage: J0,J5,J1 = J.peirce_decomposition(J.one())
973 sage: J0.dimension() == 0 and J5.dimension() == 0
975 sage: J1.superalgebra() == J and J1.dimension() == J.dimension()
978 The decomposition is into eigenspaces, and its components are
979 therefore necessarily orthogonal. Moreover, the identity
980 elements in the two subalgebras are the projections onto their
981 respective subspaces of the superalgebra's identity element::
983 sage: set_random_seed()
984 sage: J = random_eja()
985 sage: x = J.random_element()
986 sage: if not J.is_trivial():
987 ....: while x.is_nilpotent():
988 ....: x = J.random_element()
989 sage: c = x.subalgebra_idempotent()
990 sage: J0,J5,J1 = J.peirce_decomposition(c)
992 sage: for (w,y,z) in zip(J0.basis(), J5.basis(), J1.basis()):
993 ....: w = w.superalgebra_element()
994 ....: y = J.from_vector(y)
995 ....: z = z.superalgebra_element()
996 ....: ipsum += w.inner_product(y).abs()
997 ....: ipsum += w.inner_product(z).abs()
998 ....: ipsum += y.inner_product(z).abs()
1001 sage: J1(c) == J1.one()
1003 sage: J0(J.one() - c) == J0.one()
1007 if not c
.is_idempotent():
1008 raise ValueError("element is not idempotent: %s" % c
)
1010 from mjo
.eja
.eja_subalgebra
import FiniteDimensionalEJASubalgebra
1012 # Default these to what they should be if they turn out to be
1013 # trivial, because eigenspaces_left() won't return eigenvalues
1014 # corresponding to trivial spaces (e.g. it returns only the
1015 # eigenspace corresponding to lambda=1 if you take the
1016 # decomposition relative to the identity element).
1017 trivial
= FiniteDimensionalEJASubalgebra(self
, ())
1018 J0
= trivial
# eigenvalue zero
1019 J5
= VectorSpace(self
.base_ring(), 0) # eigenvalue one-half
1020 J1
= trivial
# eigenvalue one
1022 for (eigval
, eigspace
) in c
.operator().matrix().right_eigenspaces():
1023 if eigval
== ~
(self
.base_ring()(2)):
1026 gens
= tuple( self
.from_vector(b
) for b
in eigspace
.basis() )
1027 subalg
= FiniteDimensionalEJASubalgebra(self
,
1035 raise ValueError("unexpected eigenvalue: %s" % eigval
)
1040 def random_element(self
, thorough
=False):
1042 Return a random element of this algebra.
1044 Our algebra superclass method only returns a linear
1045 combination of at most two basis elements. We instead
1046 want the vector space "random element" method that
1047 returns a more diverse selection.
1051 - ``thorough`` -- (boolean; default False) whether or not we
1052 should generate irrational coefficients for the random
1053 element when our base ring is irrational; this slows the
1054 algebra operations to a crawl, but any truly random method
1058 # For a general base ring... maybe we can trust this to do the
1059 # right thing? Unlikely, but.
1060 V
= self
.vector_space()
1061 v
= V
.random_element()
1063 if self
.base_ring() is AA
:
1064 # The "random element" method of the algebraic reals is
1065 # stupid at the moment, and only returns integers between
1066 # -2 and 2, inclusive:
1068 # https://trac.sagemath.org/ticket/30875
1070 # Instead, we implement our own "random vector" method,
1071 # and then coerce that into the algebra. We use the vector
1072 # space degree here instead of the dimension because a
1073 # subalgebra could (for example) be spanned by only two
1074 # vectors, each with five coordinates. We need to
1075 # generate all five coordinates.
1077 v
*= QQbar
.random_element().real()
1079 v
*= QQ
.random_element()
1081 return self
.from_vector(V
.coordinate_vector(v
))
1083 def random_elements(self
, count
, thorough
=False):
1085 Return ``count`` random elements as a tuple.
1089 - ``thorough`` -- (boolean; default False) whether or not we
1090 should generate irrational coefficients for the random
1091 elements when our base ring is irrational; this slows the
1092 algebra operations to a crawl, but any truly random method
1097 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1101 sage: J = JordanSpinEJA(3)
1102 sage: x,y,z = J.random_elements(3)
1103 sage: all( [ x in J, y in J, z in J ])
1105 sage: len( J.random_elements(10) ) == 10
1109 return tuple( self
.random_element(thorough
)
1110 for idx
in range(count
) )
1114 def _charpoly_coefficients(self
):
1116 The `r` polynomial coefficients of the "characteristic polynomial
1121 sage: from mjo.eja.eja_algebra import random_eja
1125 The theory shows that these are all homogeneous polynomials of
1128 sage: set_random_seed()
1129 sage: J = random_eja()
1130 sage: all(p.is_homogeneous() for p in J._charpoly_coefficients())
1134 n
= self
.dimension()
1135 R
= self
.coordinate_polynomial_ring()
1137 F
= R
.fraction_field()
1140 # From a result in my book, these are the entries of the
1141 # basis representation of L_x.
1142 return sum( vars[k
]*self
.gens()[k
].operator().matrix()[i
,j
]
1145 L_x
= matrix(F
, n
, n
, L_x_i_j
)
1148 if self
.rank
.is_in_cache():
1150 # There's no need to pad the system with redundant
1151 # columns if we *know* they'll be redundant.
1154 # Compute an extra power in case the rank is equal to
1155 # the dimension (otherwise, we would stop at x^(r-1)).
1156 x_powers
= [ (L_x
**k
)*self
.one().to_vector()
1157 for k
in range(n
+1) ]
1158 A
= matrix
.column(F
, x_powers
[:n
])
1159 AE
= A
.extended_echelon_form()
1166 # The theory says that only the first "r" coefficients are
1167 # nonzero, and they actually live in the original polynomial
1168 # ring and not the fraction field. We negate them because in
1169 # the actual characteristic polynomial, they get moved to the
1170 # other side where x^r lives. We don't bother to trim A_rref
1171 # down to a square matrix and solve the resulting system,
1172 # because the upper-left r-by-r portion of A_rref is
1173 # guaranteed to be the identity matrix, so e.g.
1175 # A_rref.solve_right(Y)
1177 # would just be returning Y.
1178 return (-E
*b
)[:r
].change_ring(R
)
1183 Return the rank of this EJA.
1185 This is a cached method because we know the rank a priori for
1186 all of the algebras we can construct. Thus we can avoid the
1187 expensive ``_charpoly_coefficients()`` call unless we truly
1188 need to compute the whole characteristic polynomial.
1192 sage: from mjo.eja.eja_algebra import (HadamardEJA,
1193 ....: JordanSpinEJA,
1194 ....: RealSymmetricEJA,
1195 ....: ComplexHermitianEJA,
1196 ....: QuaternionHermitianEJA,
1201 The rank of the Jordan spin algebra is always two::
1203 sage: JordanSpinEJA(2).rank()
1205 sage: JordanSpinEJA(3).rank()
1207 sage: JordanSpinEJA(4).rank()
1210 The rank of the `n`-by-`n` Hermitian real, complex, or
1211 quaternion matrices is `n`::
1213 sage: RealSymmetricEJA(4).rank()
1215 sage: ComplexHermitianEJA(3).rank()
1217 sage: QuaternionHermitianEJA(2).rank()
1222 Ensure that every EJA that we know how to construct has a
1223 positive integer rank, unless the algebra is trivial in
1224 which case its rank will be zero::
1226 sage: set_random_seed()
1227 sage: J = random_eja()
1231 sage: r > 0 or (r == 0 and J.is_trivial())
1234 Ensure that computing the rank actually works, since the ranks
1235 of all simple algebras are known and will be cached by default::
1237 sage: set_random_seed() # long time
1238 sage: J = random_eja() # long time
1239 sage: cached = J.rank() # long time
1240 sage: J.rank.clear_cache() # long time
1241 sage: J.rank() == cached # long time
1245 return len(self
._charpoly
_coefficients
())
1248 def vector_space(self
):
1250 Return the vector space that underlies this algebra.
1254 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1258 sage: J = RealSymmetricEJA(2)
1259 sage: J.vector_space()
1260 Vector space of dimension 3 over...
1263 return self
.zero().to_vector().parent().ambient_vector_space()
1266 Element
= FiniteDimensionalEJAElement
1268 class RationalBasisEJA(FiniteDimensionalEJA
):
1270 New class for algebras whose supplied basis elements have all rational entries.
1274 sage: from mjo.eja.eja_algebra import BilinearFormEJA
1278 The supplied basis is orthonormalized by default::
1280 sage: B = matrix(QQ, [[1, 0, 0], [0, 25, -32], [0, -32, 41]])
1281 sage: J = BilinearFormEJA(B)
1282 sage: J.matrix_basis()
1299 # Abuse the check_field parameter to check that the entries of
1300 # out basis (in ambient coordinates) are in the field QQ.
1301 if not all( all(b_i
in QQ
for b_i
in b
.list()) for b
in basis
):
1302 raise TypeError("basis not rational")
1304 self
._rational
_algebra
= None
1306 # There's no point in constructing the extra algebra if this
1307 # one is already rational.
1309 # Note: the same Jordan and inner-products work here,
1310 # because they are necessarily defined with respect to
1311 # ambient coordinates and not any particular basis.
1312 self
._rational
_algebra
= FiniteDimensionalEJA(
1317 orthonormalize
=False,
1321 super().__init
__(basis
,
1325 check_field
=check_field
,
1329 def _charpoly_coefficients(self
):
1333 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
1334 ....: JordanSpinEJA)
1338 The base ring of the resulting polynomial coefficients is what
1339 it should be, and not the rationals (unless the algebra was
1340 already over the rationals)::
1342 sage: J = JordanSpinEJA(3)
1343 sage: J._charpoly_coefficients()
1344 (X1^2 - X2^2 - X3^2, -2*X1)
1345 sage: a0 = J._charpoly_coefficients()[0]
1347 Algebraic Real Field
1348 sage: a0.base_ring()
1349 Algebraic Real Field
1352 if self
._rational
_algebra
is None:
1353 # There's no need to construct *another* algebra over the
1354 # rationals if this one is already over the
1355 # rationals. Likewise, if we never orthonormalized our
1356 # basis, we might as well just use the given one.
1357 return super()._charpoly
_coefficients
()
1359 # Do the computation over the rationals. The answer will be
1360 # the same, because all we've done is a change of basis.
1361 # Then, change back from QQ to our real base ring
1362 a
= ( a_i
.change_ring(self
.base_ring())
1363 for a_i
in self
._rational
_algebra
._charpoly
_coefficients
() )
1365 if self
._deortho
_matrix
is None:
1366 # This can happen if our base ring was, say, AA and we
1367 # chose not to (or didn't need to) orthonormalize. It's
1368 # still faster to do the computations over QQ even if
1369 # the numbers in the boxes stay the same.
1372 # Otherwise, convert the coordinate variables back to the
1373 # deorthonormalized ones.
1374 R
= self
.coordinate_polynomial_ring()
1375 from sage
.modules
.free_module_element
import vector
1376 X
= vector(R
, R
.gens())
1377 BX
= self
._deortho
_matrix
*X
1379 subs_dict
= { X[i]: BX[i] for i in range(len(X)) }
1380 return tuple( a_i
.subs(subs_dict
) for a_i
in a
)
1382 class ConcreteEJA(RationalBasisEJA
):
1384 A class for the Euclidean Jordan algebras that we know by name.
1386 These are the Jordan algebras whose basis, multiplication table,
1387 rank, and so on are known a priori. More to the point, they are
1388 the Euclidean Jordan algebras for which we are able to conjure up
1389 a "random instance."
1393 sage: from mjo.eja.eja_algebra import ConcreteEJA
1397 Our basis is normalized with respect to the algebra's inner
1398 product, unless we specify otherwise::
1400 sage: set_random_seed()
1401 sage: J = ConcreteEJA.random_instance()
1402 sage: all( b.norm() == 1 for b in J.gens() )
1405 Since our basis is orthonormal with respect to the algebra's inner
1406 product, and since we know that this algebra is an EJA, any
1407 left-multiplication operator's matrix will be symmetric because
1408 natural->EJA basis representation is an isometry and within the
1409 EJA the operator is self-adjoint by the Jordan axiom::
1411 sage: set_random_seed()
1412 sage: J = ConcreteEJA.random_instance()
1413 sage: x = J.random_element()
1414 sage: x.operator().is_self_adjoint()
1419 def _max_random_instance_size():
1421 Return an integer "size" that is an upper bound on the size of
1422 this algebra when it is used in a random test
1423 case. Unfortunately, the term "size" is ambiguous -- when
1424 dealing with `R^n` under either the Hadamard or Jordan spin
1425 product, the "size" refers to the dimension `n`. When dealing
1426 with a matrix algebra (real symmetric or complex/quaternion
1427 Hermitian), it refers to the size of the matrix, which is far
1428 less than the dimension of the underlying vector space.
1430 This method must be implemented in each subclass.
1432 raise NotImplementedError
1435 def random_instance(cls
, *args
, **kwargs
):
1437 Return a random instance of this type of algebra.
1439 This method should be implemented in each subclass.
1441 from sage
.misc
.prandom
import choice
1442 eja_class
= choice(cls
.__subclasses
__())
1444 # These all bubble up to the RationalBasisEJA superclass
1445 # constructor, so any (kw)args valid there are also valid
1447 return eja_class
.random_instance(*args
, **kwargs
)
1452 def dimension_over_reals():
1454 The dimension of this matrix's base ring over the reals.
1456 The reals are dimension one over themselves, obviously; that's
1457 just `\mathbb{R}^{1}`. Likewise, the complex numbers `a + bi`
1458 have dimension two. Finally, the quaternions have dimension
1459 four over the reals.
1461 This is used to determine the size of the matrix returned from
1462 :meth:`real_embed`, among other things.
1464 raise NotImplementedError
1467 def real_embed(cls
,M
):
1469 Embed the matrix ``M`` into a space of real matrices.
1471 The matrix ``M`` can have entries in any field at the moment:
1472 the real numbers, complex numbers, or quaternions. And although
1473 they are not a field, we can probably support octonions at some
1474 point, too. This function returns a real matrix that "acts like"
1475 the original with respect to matrix multiplication; i.e.
1477 real_embed(M*N) = real_embed(M)*real_embed(N)
1480 if M
.ncols() != M
.nrows():
1481 raise ValueError("the matrix 'M' must be square")
1486 def real_unembed(cls
,M
):
1488 The inverse of :meth:`real_embed`.
1490 if M
.ncols() != M
.nrows():
1491 raise ValueError("the matrix 'M' must be square")
1492 if not ZZ(M
.nrows()).mod(cls
.dimension_over_reals()).is_zero():
1493 raise ValueError("the matrix 'M' must be a real embedding")
1497 def jordan_product(X
,Y
):
1498 return (X
*Y
+ Y
*X
)/2
1501 def trace_inner_product(cls
,X
,Y
):
1503 Compute the trace inner-product of two real-embeddings.
1507 sage: from mjo.eja.eja_algebra import (RealSymmetricEJA,
1508 ....: ComplexHermitianEJA,
1509 ....: QuaternionHermitianEJA)
1513 This gives the same answer as it would if we computed the trace
1514 from the unembedded (original) matrices::
1516 sage: set_random_seed()
1517 sage: J = RealSymmetricEJA.random_instance()
1518 sage: x,y = J.random_elements(2)
1519 sage: Xe = x.to_matrix()
1520 sage: Ye = y.to_matrix()
1521 sage: X = J.real_unembed(Xe)
1522 sage: Y = J.real_unembed(Ye)
1523 sage: expected = (X*Y).trace()
1524 sage: actual = J.trace_inner_product(Xe,Ye)
1525 sage: actual == expected
1530 sage: set_random_seed()
1531 sage: J = ComplexHermitianEJA.random_instance()
1532 sage: x,y = J.random_elements(2)
1533 sage: Xe = x.to_matrix()
1534 sage: Ye = y.to_matrix()
1535 sage: X = J.real_unembed(Xe)
1536 sage: Y = J.real_unembed(Ye)
1537 sage: expected = (X*Y).trace().real()
1538 sage: actual = J.trace_inner_product(Xe,Ye)
1539 sage: actual == expected
1544 sage: set_random_seed()
1545 sage: J = QuaternionHermitianEJA.random_instance()
1546 sage: x,y = J.random_elements(2)
1547 sage: Xe = x.to_matrix()
1548 sage: Ye = y.to_matrix()
1549 sage: X = J.real_unembed(Xe)
1550 sage: Y = J.real_unembed(Ye)
1551 sage: expected = (X*Y).trace().coefficient_tuple()[0]
1552 sage: actual = J.trace_inner_product(Xe,Ye)
1553 sage: actual == expected
1557 Xu
= cls
.real_unembed(X
)
1558 Yu
= cls
.real_unembed(Y
)
1559 tr
= (Xu
*Yu
).trace()
1562 # Works in QQ, AA, RDF, et cetera.
1564 except AttributeError:
1565 # A quaternion doesn't have a real() method, but does
1566 # have coefficient_tuple() method that returns the
1567 # coefficients of 1, i, j, and k -- in that order.
1568 return tr
.coefficient_tuple()[0]
1571 class RealMatrixEJA(MatrixEJA
):
1573 def dimension_over_reals():
1577 class RealSymmetricEJA(ConcreteEJA
, RealMatrixEJA
):
1579 The rank-n simple EJA consisting of real symmetric n-by-n
1580 matrices, the usual symmetric Jordan product, and the trace inner
1581 product. It has dimension `(n^2 + n)/2` over the reals.
1585 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1589 sage: J = RealSymmetricEJA(2)
1590 sage: e0, e1, e2 = J.gens()
1598 In theory, our "field" can be any subfield of the reals::
1600 sage: RealSymmetricEJA(2, field=RDF)
1601 Euclidean Jordan algebra of dimension 3 over Real Double Field
1602 sage: RealSymmetricEJA(2, field=RR)
1603 Euclidean Jordan algebra of dimension 3 over Real Field with
1604 53 bits of precision
1608 The dimension of this algebra is `(n^2 + n) / 2`::
1610 sage: set_random_seed()
1611 sage: n_max = RealSymmetricEJA._max_random_instance_size()
1612 sage: n = ZZ.random_element(1, n_max)
1613 sage: J = RealSymmetricEJA(n)
1614 sage: J.dimension() == (n^2 + n)/2
1617 The Jordan multiplication is what we think it is::
1619 sage: set_random_seed()
1620 sage: J = RealSymmetricEJA.random_instance()
1621 sage: x,y = J.random_elements(2)
1622 sage: actual = (x*y).to_matrix()
1623 sage: X = x.to_matrix()
1624 sage: Y = y.to_matrix()
1625 sage: expected = (X*Y + Y*X)/2
1626 sage: actual == expected
1628 sage: J(expected) == x*y
1631 We can change the generator prefix::
1633 sage: RealSymmetricEJA(3, prefix='q').gens()
1634 (q0, q1, q2, q3, q4, q5)
1636 We can construct the (trivial) algebra of rank zero::
1638 sage: RealSymmetricEJA(0)
1639 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1643 def _denormalized_basis(cls
, n
):
1645 Return a basis for the space of real symmetric n-by-n matrices.
1649 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1653 sage: set_random_seed()
1654 sage: n = ZZ.random_element(1,5)
1655 sage: B = RealSymmetricEJA._denormalized_basis(n)
1656 sage: all( M.is_symmetric() for M in B)
1660 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1664 for j
in range(i
+1):
1665 Eij
= matrix(ZZ
, n
, lambda k
,l
: k
==i
and l
==j
)
1669 Sij
= Eij
+ Eij
.transpose()
1675 def _max_random_instance_size():
1676 return 4 # Dimension 10
1679 def random_instance(cls
, **kwargs
):
1681 Return a random instance of this type of algebra.
1683 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1684 return cls(n
, **kwargs
)
1686 def __init__(self
, n
, **kwargs
):
1687 # We know this is a valid EJA, but will double-check
1688 # if the user passes check_axioms=True.
1689 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
1691 super(RealSymmetricEJA
, self
).__init
__(self
._denormalized
_basis
(n
),
1692 self
.jordan_product
,
1693 self
.trace_inner_product
,
1696 # TODO: this could be factored out somehow, but is left here
1697 # because the MatrixEJA is not presently a subclass of the
1698 # FDEJA class that defines rank() and one().
1699 self
.rank
.set_cache(n
)
1700 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
1701 self
.one
.set_cache(self(idV
))
1705 class ComplexMatrixEJA(MatrixEJA
):
1706 # A manual dictionary-cache for the complex_extension() method,
1707 # since apparently @classmethods can't also be @cached_methods.
1708 _complex_extension
= {}
1711 def complex_extension(cls
,field
):
1713 The complex field that we embed/unembed, as an extension
1714 of the given ``field``.
1716 if field
in cls
._complex
_extension
:
1717 return cls
._complex
_extension
[field
]
1719 # Sage doesn't know how to adjoin the complex "i" (the root of
1720 # x^2 + 1) to a field in a general way. Here, we just enumerate
1721 # all of the cases that I have cared to support so far.
1723 # Sage doesn't know how to embed AA into QQbar, i.e. how
1724 # to adjoin sqrt(-1) to AA.
1726 elif not field
.is_exact():
1728 F
= field
.complex_field()
1730 # Works for QQ and... maybe some other fields.
1731 R
= PolynomialRing(field
, 'z')
1733 F
= field
.extension(z
**2 + 1, 'I', embedding
=CLF(-1).sqrt())
1735 cls
._complex
_extension
[field
] = F
1739 def dimension_over_reals():
1743 def real_embed(cls
,M
):
1745 Embed the n-by-n complex matrix ``M`` into the space of real
1746 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1747 bi` to the block matrix ``[[a,b],[-b,a]]``.
1751 sage: from mjo.eja.eja_algebra import ComplexMatrixEJA
1755 sage: F = QuadraticField(-1, 'I')
1756 sage: x1 = F(4 - 2*i)
1757 sage: x2 = F(1 + 2*i)
1760 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1761 sage: ComplexMatrixEJA.real_embed(M)
1770 Embedding is a homomorphism (isomorphism, in fact)::
1772 sage: set_random_seed()
1773 sage: n = ZZ.random_element(3)
1774 sage: F = QuadraticField(-1, 'I')
1775 sage: X = random_matrix(F, n)
1776 sage: Y = random_matrix(F, n)
1777 sage: Xe = ComplexMatrixEJA.real_embed(X)
1778 sage: Ye = ComplexMatrixEJA.real_embed(Y)
1779 sage: XYe = ComplexMatrixEJA.real_embed(X*Y)
1784 super(ComplexMatrixEJA
,cls
).real_embed(M
)
1787 # We don't need any adjoined elements...
1788 field
= M
.base_ring().base_ring()
1794 blocks
.append(matrix(field
, 2, [ [ a
, b
],
1797 return matrix
.block(field
, n
, blocks
)
1801 def real_unembed(cls
,M
):
1803 The inverse of _embed_complex_matrix().
1807 sage: from mjo.eja.eja_algebra import ComplexMatrixEJA
1811 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1812 ....: [-2, 1, -4, 3],
1813 ....: [ 9, 10, 11, 12],
1814 ....: [-10, 9, -12, 11] ])
1815 sage: ComplexMatrixEJA.real_unembed(A)
1817 [ 10*I + 9 12*I + 11]
1821 Unembedding is the inverse of embedding::
1823 sage: set_random_seed()
1824 sage: F = QuadraticField(-1, 'I')
1825 sage: M = random_matrix(F, 3)
1826 sage: Me = ComplexMatrixEJA.real_embed(M)
1827 sage: ComplexMatrixEJA.real_unembed(Me) == M
1831 super(ComplexMatrixEJA
,cls
).real_unembed(M
)
1833 d
= cls
.dimension_over_reals()
1834 F
= cls
.complex_extension(M
.base_ring())
1837 # Go top-left to bottom-right (reading order), converting every
1838 # 2-by-2 block we see to a single complex element.
1840 for k
in range(n
/d
):
1841 for j
in range(n
/d
):
1842 submat
= M
[d
*k
:d
*k
+d
,d
*j
:d
*j
+d
]
1843 if submat
[0,0] != submat
[1,1]:
1844 raise ValueError('bad on-diagonal submatrix')
1845 if submat
[0,1] != -submat
[1,0]:
1846 raise ValueError('bad off-diagonal submatrix')
1847 z
= submat
[0,0] + submat
[0,1]*i
1850 return matrix(F
, n
/d
, elements
)
1853 class ComplexHermitianEJA(ConcreteEJA
, ComplexMatrixEJA
):
1855 The rank-n simple EJA consisting of complex Hermitian n-by-n
1856 matrices over the real numbers, the usual symmetric Jordan product,
1857 and the real-part-of-trace inner product. It has dimension `n^2` over
1862 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1866 In theory, our "field" can be any subfield of the reals::
1868 sage: ComplexHermitianEJA(2, field=RDF)
1869 Euclidean Jordan algebra of dimension 4 over Real Double Field
1870 sage: ComplexHermitianEJA(2, field=RR)
1871 Euclidean Jordan algebra of dimension 4 over Real Field with
1872 53 bits of precision
1876 The dimension of this algebra is `n^2`::
1878 sage: set_random_seed()
1879 sage: n_max = ComplexHermitianEJA._max_random_instance_size()
1880 sage: n = ZZ.random_element(1, n_max)
1881 sage: J = ComplexHermitianEJA(n)
1882 sage: J.dimension() == n^2
1885 The Jordan multiplication is what we think it is::
1887 sage: set_random_seed()
1888 sage: J = ComplexHermitianEJA.random_instance()
1889 sage: x,y = J.random_elements(2)
1890 sage: actual = (x*y).to_matrix()
1891 sage: X = x.to_matrix()
1892 sage: Y = y.to_matrix()
1893 sage: expected = (X*Y + Y*X)/2
1894 sage: actual == expected
1896 sage: J(expected) == x*y
1899 We can change the generator prefix::
1901 sage: ComplexHermitianEJA(2, prefix='z').gens()
1904 We can construct the (trivial) algebra of rank zero::
1906 sage: ComplexHermitianEJA(0)
1907 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1912 def _denormalized_basis(cls
, n
):
1914 Returns a basis for the space of complex Hermitian n-by-n matrices.
1916 Why do we embed these? Basically, because all of numerical linear
1917 algebra assumes that you're working with vectors consisting of `n`
1918 entries from a field and scalars from the same field. There's no way
1919 to tell SageMath that (for example) the vectors contain complex
1920 numbers, while the scalar field is real.
1924 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1928 sage: set_random_seed()
1929 sage: n = ZZ.random_element(1,5)
1930 sage: B = ComplexHermitianEJA._denormalized_basis(n)
1931 sage: all( M.is_symmetric() for M in B)
1936 R
= PolynomialRing(field
, 'z')
1938 F
= field
.extension(z
**2 + 1, 'I')
1941 # This is like the symmetric case, but we need to be careful:
1943 # * We want conjugate-symmetry, not just symmetry.
1944 # * The diagonal will (as a result) be real.
1947 Eij
= matrix
.zero(F
,n
)
1949 for j
in range(i
+1):
1953 Sij
= cls
.real_embed(Eij
)
1956 # The second one has a minus because it's conjugated.
1957 Eij
[j
,i
] = 1 # Eij = Eij + Eij.transpose()
1958 Sij_real
= cls
.real_embed(Eij
)
1960 # Eij = I*Eij - I*Eij.transpose()
1963 Sij_imag
= cls
.real_embed(Eij
)
1969 # Since we embedded these, we can drop back to the "field" that we
1970 # started with instead of the complex extension "F".
1971 return tuple( s
.change_ring(field
) for s
in S
)
1974 def __init__(self
, n
, **kwargs
):
1975 # We know this is a valid EJA, but will double-check
1976 # if the user passes check_axioms=True.
1977 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
1979 super(ComplexHermitianEJA
, self
).__init
__(self
._denormalized
_basis
(n
),
1980 self
.jordan_product
,
1981 self
.trace_inner_product
,
1983 # TODO: this could be factored out somehow, but is left here
1984 # because the MatrixEJA is not presently a subclass of the
1985 # FDEJA class that defines rank() and one().
1986 self
.rank
.set_cache(n
)
1987 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
1988 self
.one
.set_cache(self(idV
))
1991 def _max_random_instance_size():
1992 return 3 # Dimension 9
1995 def random_instance(cls
, **kwargs
):
1997 Return a random instance of this type of algebra.
1999 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2000 return cls(n
, **kwargs
)
2002 class QuaternionMatrixEJA(MatrixEJA
):
2004 # A manual dictionary-cache for the quaternion_extension() method,
2005 # since apparently @classmethods can't also be @cached_methods.
2006 _quaternion_extension
= {}
2009 def quaternion_extension(cls
,field
):
2011 The quaternion field that we embed/unembed, as an extension
2012 of the given ``field``.
2014 if field
in cls
._quaternion
_extension
:
2015 return cls
._quaternion
_extension
[field
]
2017 Q
= QuaternionAlgebra(field
,-1,-1)
2019 cls
._quaternion
_extension
[field
] = Q
2023 def dimension_over_reals():
2027 def real_embed(cls
,M
):
2029 Embed the n-by-n quaternion matrix ``M`` into the space of real
2030 matrices of size 4n-by-4n by first sending each quaternion entry `z
2031 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
2032 c+di],[-c + di, a-bi]]`, and then embedding those into a real
2037 sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA
2041 sage: Q = QuaternionAlgebra(QQ,-1,-1)
2042 sage: i,j,k = Q.gens()
2043 sage: x = 1 + 2*i + 3*j + 4*k
2044 sage: M = matrix(Q, 1, [[x]])
2045 sage: QuaternionMatrixEJA.real_embed(M)
2051 Embedding is a homomorphism (isomorphism, in fact)::
2053 sage: set_random_seed()
2054 sage: n = ZZ.random_element(2)
2055 sage: Q = QuaternionAlgebra(QQ,-1,-1)
2056 sage: X = random_matrix(Q, n)
2057 sage: Y = random_matrix(Q, n)
2058 sage: Xe = QuaternionMatrixEJA.real_embed(X)
2059 sage: Ye = QuaternionMatrixEJA.real_embed(Y)
2060 sage: XYe = QuaternionMatrixEJA.real_embed(X*Y)
2065 super(QuaternionMatrixEJA
,cls
).real_embed(M
)
2066 quaternions
= M
.base_ring()
2069 F
= QuadraticField(-1, 'I')
2074 t
= z
.coefficient_tuple()
2079 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
2080 [-c
+ d
*i
, a
- b
*i
]])
2081 realM
= ComplexMatrixEJA
.real_embed(cplxM
)
2082 blocks
.append(realM
)
2084 # We should have real entries by now, so use the realest field
2085 # we've got for the return value.
2086 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
2091 def real_unembed(cls
,M
):
2093 The inverse of _embed_quaternion_matrix().
2097 sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA
2101 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
2102 ....: [-2, 1, -4, 3],
2103 ....: [-3, 4, 1, -2],
2104 ....: [-4, -3, 2, 1]])
2105 sage: QuaternionMatrixEJA.real_unembed(M)
2106 [1 + 2*i + 3*j + 4*k]
2110 Unembedding is the inverse of embedding::
2112 sage: set_random_seed()
2113 sage: Q = QuaternionAlgebra(QQ, -1, -1)
2114 sage: M = random_matrix(Q, 3)
2115 sage: Me = QuaternionMatrixEJA.real_embed(M)
2116 sage: QuaternionMatrixEJA.real_unembed(Me) == M
2120 super(QuaternionMatrixEJA
,cls
).real_unembed(M
)
2122 d
= cls
.dimension_over_reals()
2124 # Use the base ring of the matrix to ensure that its entries can be
2125 # multiplied by elements of the quaternion algebra.
2126 Q
= cls
.quaternion_extension(M
.base_ring())
2129 # Go top-left to bottom-right (reading order), converting every
2130 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
2133 for l
in range(n
/d
):
2134 for m
in range(n
/d
):
2135 submat
= ComplexMatrixEJA
.real_unembed(
2136 M
[d
*l
:d
*l
+d
,d
*m
:d
*m
+d
] )
2137 if submat
[0,0] != submat
[1,1].conjugate():
2138 raise ValueError('bad on-diagonal submatrix')
2139 if submat
[0,1] != -submat
[1,0].conjugate():
2140 raise ValueError('bad off-diagonal submatrix')
2141 z
= submat
[0,0].real()
2142 z
+= submat
[0,0].imag()*i
2143 z
+= submat
[0,1].real()*j
2144 z
+= submat
[0,1].imag()*k
2147 return matrix(Q
, n
/d
, elements
)
2150 class QuaternionHermitianEJA(ConcreteEJA
, QuaternionMatrixEJA
):
2152 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
2153 matrices, the usual symmetric Jordan product, and the
2154 real-part-of-trace inner product. It has dimension `2n^2 - n` over
2159 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2163 In theory, our "field" can be any subfield of the reals::
2165 sage: QuaternionHermitianEJA(2, field=RDF)
2166 Euclidean Jordan algebra of dimension 6 over Real Double Field
2167 sage: QuaternionHermitianEJA(2, field=RR)
2168 Euclidean Jordan algebra of dimension 6 over Real Field with
2169 53 bits of precision
2173 The dimension of this algebra is `2*n^2 - n`::
2175 sage: set_random_seed()
2176 sage: n_max = QuaternionHermitianEJA._max_random_instance_size()
2177 sage: n = ZZ.random_element(1, n_max)
2178 sage: J = QuaternionHermitianEJA(n)
2179 sage: J.dimension() == 2*(n^2) - n
2182 The Jordan multiplication is what we think it is::
2184 sage: set_random_seed()
2185 sage: J = QuaternionHermitianEJA.random_instance()
2186 sage: x,y = J.random_elements(2)
2187 sage: actual = (x*y).to_matrix()
2188 sage: X = x.to_matrix()
2189 sage: Y = y.to_matrix()
2190 sage: expected = (X*Y + Y*X)/2
2191 sage: actual == expected
2193 sage: J(expected) == x*y
2196 We can change the generator prefix::
2198 sage: QuaternionHermitianEJA(2, prefix='a').gens()
2199 (a0, a1, a2, a3, a4, a5)
2201 We can construct the (trivial) algebra of rank zero::
2203 sage: QuaternionHermitianEJA(0)
2204 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2208 def _denormalized_basis(cls
, n
):
2210 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
2212 Why do we embed these? Basically, because all of numerical
2213 linear algebra assumes that you're working with vectors consisting
2214 of `n` entries from a field and scalars from the same field. There's
2215 no way to tell SageMath that (for example) the vectors contain
2216 complex numbers, while the scalar field is real.
2220 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2224 sage: set_random_seed()
2225 sage: n = ZZ.random_element(1,5)
2226 sage: B = QuaternionHermitianEJA._denormalized_basis(n)
2227 sage: all( M.is_symmetric() for M in B )
2232 Q
= QuaternionAlgebra(QQ
,-1,-1)
2235 # This is like the symmetric case, but we need to be careful:
2237 # * We want conjugate-symmetry, not just symmetry.
2238 # * The diagonal will (as a result) be real.
2241 Eij
= matrix
.zero(Q
,n
)
2243 for j
in range(i
+1):
2247 Sij
= cls
.real_embed(Eij
)
2250 # The second, third, and fourth ones have a minus
2251 # because they're conjugated.
2252 # Eij = Eij + Eij.transpose()
2254 Sij_real
= cls
.real_embed(Eij
)
2256 # Eij = I*(Eij - Eij.transpose())
2259 Sij_I
= cls
.real_embed(Eij
)
2261 # Eij = J*(Eij - Eij.transpose())
2264 Sij_J
= cls
.real_embed(Eij
)
2266 # Eij = K*(Eij - Eij.transpose())
2269 Sij_K
= cls
.real_embed(Eij
)
2275 # Since we embedded these, we can drop back to the "field" that we
2276 # started with instead of the quaternion algebra "Q".
2277 return tuple( s
.change_ring(field
) for s
in S
)
2280 def __init__(self
, n
, **kwargs
):
2281 # We know this is a valid EJA, but will double-check
2282 # if the user passes check_axioms=True.
2283 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2285 super(QuaternionHermitianEJA
, self
).__init
__(self
._denormalized
_basis
(n
),
2286 self
.jordan_product
,
2287 self
.trace_inner_product
,
2289 # TODO: this could be factored out somehow, but is left here
2290 # because the MatrixEJA is not presently a subclass of the
2291 # FDEJA class that defines rank() and one().
2292 self
.rank
.set_cache(n
)
2293 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
2294 self
.one
.set_cache(self(idV
))
2298 def _max_random_instance_size():
2300 The maximum rank of a random QuaternionHermitianEJA.
2302 return 2 # Dimension 6
2305 def random_instance(cls
, **kwargs
):
2307 Return a random instance of this type of algebra.
2309 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2310 return cls(n
, **kwargs
)
2313 class HadamardEJA(ConcreteEJA
):
2315 Return the Euclidean Jordan Algebra corresponding to the set
2316 `R^n` under the Hadamard product.
2318 Note: this is nothing more than the Cartesian product of ``n``
2319 copies of the spin algebra. Once Cartesian product algebras
2320 are implemented, this can go.
2324 sage: from mjo.eja.eja_algebra import HadamardEJA
2328 This multiplication table can be verified by hand::
2330 sage: J = HadamardEJA(3)
2331 sage: e0,e1,e2 = J.gens()
2347 We can change the generator prefix::
2349 sage: HadamardEJA(3, prefix='r').gens()
2353 def __init__(self
, n
, **kwargs
):
2355 jordan_product
= lambda x
,y
: x
2356 inner_product
= lambda x
,y
: x
2358 def jordan_product(x
,y
):
2360 return P( xi
*yi
for (xi
,yi
) in zip(x
,y
) )
2362 def inner_product(x
,y
):
2365 # New defaults for keyword arguments. Don't orthonormalize
2366 # because our basis is already orthonormal with respect to our
2367 # inner-product. Don't check the axioms, because we know this
2368 # is a valid EJA... but do double-check if the user passes
2369 # check_axioms=True. Note: we DON'T override the "check_field"
2370 # default here, because the user can pass in a field!
2371 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2372 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2374 column_basis
= tuple( b
.column() for b
in FreeModule(ZZ
, n
).basis() )
2375 super().__init
__(column_basis
, jordan_product
, inner_product
, **kwargs
)
2376 self
.rank
.set_cache(n
)
2379 self
.one
.set_cache( self
.zero() )
2381 self
.one
.set_cache( sum(self
.gens()) )
2384 def _max_random_instance_size():
2386 The maximum dimension of a random HadamardEJA.
2391 def random_instance(cls
, **kwargs
):
2393 Return a random instance of this type of algebra.
2395 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2396 return cls(n
, **kwargs
)
2399 class BilinearFormEJA(ConcreteEJA
):
2401 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2402 with the half-trace inner product and jordan product ``x*y =
2403 (<Bx,y>,y_bar>, x0*y_bar + y0*x_bar)`` where `B = 1 \times B22` is
2404 a symmetric positive-definite "bilinear form" matrix. Its
2405 dimension is the size of `B`, and it has rank two in dimensions
2406 larger than two. It reduces to the ``JordanSpinEJA`` when `B` is
2407 the identity matrix of order ``n``.
2409 We insist that the one-by-one upper-left identity block of `B` be
2410 passed in as well so that we can be passed a matrix of size zero
2411 to construct a trivial algebra.
2415 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
2416 ....: JordanSpinEJA)
2420 When no bilinear form is specified, the identity matrix is used,
2421 and the resulting algebra is the Jordan spin algebra::
2423 sage: B = matrix.identity(AA,3)
2424 sage: J0 = BilinearFormEJA(B)
2425 sage: J1 = JordanSpinEJA(3)
2426 sage: J0.multiplication_table() == J0.multiplication_table()
2429 An error is raised if the matrix `B` does not correspond to a
2430 positive-definite bilinear form::
2432 sage: B = matrix.random(QQ,2,3)
2433 sage: J = BilinearFormEJA(B)
2434 Traceback (most recent call last):
2436 ValueError: bilinear form is not positive-definite
2437 sage: B = matrix.zero(QQ,3)
2438 sage: J = BilinearFormEJA(B)
2439 Traceback (most recent call last):
2441 ValueError: bilinear form is not positive-definite
2445 We can create a zero-dimensional algebra::
2447 sage: B = matrix.identity(AA,0)
2448 sage: J = BilinearFormEJA(B)
2452 We can check the multiplication condition given in the Jordan, von
2453 Neumann, and Wigner paper (and also discussed on my "On the
2454 symmetry..." paper). Note that this relies heavily on the standard
2455 choice of basis, as does anything utilizing the bilinear form
2456 matrix. We opt not to orthonormalize the basis, because if we
2457 did, we would have to normalize the `s_{i}` in a similar manner::
2459 sage: set_random_seed()
2460 sage: n = ZZ.random_element(5)
2461 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2462 sage: B11 = matrix.identity(QQ,1)
2463 sage: B22 = M.transpose()*M
2464 sage: B = block_matrix(2,2,[ [B11,0 ],
2466 sage: J = BilinearFormEJA(B, orthonormalize=False)
2467 sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
2468 sage: V = J.vector_space()
2469 sage: sis = [ J( V([0] + (M.inverse()*ei).list()).column() )
2470 ....: for ei in eis ]
2471 sage: actual = [ sis[i]*sis[j]
2472 ....: for i in range(n-1)
2473 ....: for j in range(n-1) ]
2474 sage: expected = [ J.one() if i == j else J.zero()
2475 ....: for i in range(n-1)
2476 ....: for j in range(n-1) ]
2477 sage: actual == expected
2481 def __init__(self
, B
, **kwargs
):
2482 # The matrix "B" is supplied by the user in most cases,
2483 # so it makes sense to check whether or not its positive-
2484 # definite unless we are specifically asked not to...
2485 if ("check_axioms" not in kwargs
) or kwargs
["check_axioms"]:
2486 if not B
.is_positive_definite():
2487 raise ValueError("bilinear form is not positive-definite")
2489 # However, all of the other data for this EJA is computed
2490 # by us in manner that guarantees the axioms are
2491 # satisfied. So, again, unless we are specifically asked to
2492 # verify things, we'll skip the rest of the checks.
2493 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2495 def inner_product(x
,y
):
2496 return (y
.T
*B
*x
)[0,0]
2498 def jordan_product(x
,y
):
2504 z0
= inner_product(y
,x
)
2505 zbar
= y0
*xbar
+ x0
*ybar
2506 return P([z0
] + zbar
.list())
2509 column_basis
= tuple( b
.column() for b
in FreeModule(ZZ
, n
).basis() )
2510 super(BilinearFormEJA
, self
).__init
__(column_basis
,
2515 # The rank of this algebra is two, unless we're in a
2516 # one-dimensional ambient space (because the rank is bounded
2517 # by the ambient dimension).
2518 self
.rank
.set_cache(min(n
,2))
2521 self
.one
.set_cache( self
.zero() )
2523 self
.one
.set_cache( self
.monomial(0) )
2526 def _max_random_instance_size():
2528 The maximum dimension of a random BilinearFormEJA.
2533 def random_instance(cls
, **kwargs
):
2535 Return a random instance of this algebra.
2537 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2539 B
= matrix
.identity(ZZ
, n
)
2540 return cls(B
, **kwargs
)
2542 B11
= matrix
.identity(ZZ
, 1)
2543 M
= matrix
.random(ZZ
, n
-1)
2544 I
= matrix
.identity(ZZ
, n
-1)
2546 while alpha
.is_zero():
2547 alpha
= ZZ
.random_element().abs()
2548 B22
= M
.transpose()*M
+ alpha
*I
2550 from sage
.matrix
.special
import block_matrix
2551 B
= block_matrix(2,2, [ [B11
, ZZ(0) ],
2554 return cls(B
, **kwargs
)
2557 class JordanSpinEJA(BilinearFormEJA
):
2559 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2560 with the usual inner product and jordan product ``x*y =
2561 (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2566 sage: from mjo.eja.eja_algebra import JordanSpinEJA
2570 This multiplication table can be verified by hand::
2572 sage: J = JordanSpinEJA(4)
2573 sage: e0,e1,e2,e3 = J.gens()
2589 We can change the generator prefix::
2591 sage: JordanSpinEJA(2, prefix='B').gens()
2596 Ensure that we have the usual inner product on `R^n`::
2598 sage: set_random_seed()
2599 sage: J = JordanSpinEJA.random_instance()
2600 sage: x,y = J.random_elements(2)
2601 sage: actual = x.inner_product(y)
2602 sage: expected = x.to_vector().inner_product(y.to_vector())
2603 sage: actual == expected
2607 def __init__(self
, n
, **kwargs
):
2608 # This is a special case of the BilinearFormEJA with the
2609 # identity matrix as its bilinear form.
2610 B
= matrix
.identity(ZZ
, n
)
2612 # Don't orthonormalize because our basis is already
2613 # orthonormal with respect to our inner-product.
2614 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2616 # But also don't pass check_field=False here, because the user
2617 # can pass in a field!
2618 super(JordanSpinEJA
, self
).__init
__(B
, **kwargs
)
2621 def _max_random_instance_size():
2623 The maximum dimension of a random JordanSpinEJA.
2628 def random_instance(cls
, **kwargs
):
2630 Return a random instance of this type of algebra.
2632 Needed here to override the implementation for ``BilinearFormEJA``.
2634 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2635 return cls(n
, **kwargs
)
2638 class TrivialEJA(ConcreteEJA
):
2640 The trivial Euclidean Jordan algebra consisting of only a zero element.
2644 sage: from mjo.eja.eja_algebra import TrivialEJA
2648 sage: J = TrivialEJA()
2655 sage: 7*J.one()*12*J.one()
2657 sage: J.one().inner_product(J.one())
2659 sage: J.one().norm()
2661 sage: J.one().subalgebra_generated_by()
2662 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2667 def __init__(self
, **kwargs
):
2668 jordan_product
= lambda x
,y
: x
2669 inner_product
= lambda x
,y
: 0
2672 # New defaults for keyword arguments
2673 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2674 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2676 super(TrivialEJA
, self
).__init
__(basis
,
2680 # The rank is zero using my definition, namely the dimension of the
2681 # largest subalgebra generated by any element.
2682 self
.rank
.set_cache(0)
2683 self
.one
.set_cache( self
.zero() )
2686 def random_instance(cls
, **kwargs
):
2687 # We don't take a "size" argument so the superclass method is
2688 # inappropriate for us.
2689 return cls(**kwargs
)
2692 class CartesianProductEJA(CombinatorialFreeModule_CartesianProduct
,
2693 FiniteDimensionalEJA
):
2695 The external (orthogonal) direct sum of two or more Euclidean
2696 Jordan algebras. Every Euclidean Jordan algebra decomposes into an
2697 orthogonal direct sum of simple Euclidean Jordan algebras which is
2698 then isometric to a Cartesian product, so no generality is lost by
2699 providing only this construction.
2703 sage: from mjo.eja.eja_algebra import (random_eja,
2704 ....: CartesianProductEJA,
2706 ....: JordanSpinEJA,
2707 ....: RealSymmetricEJA)
2711 The Jordan product is inherited from our factors and implemented by
2712 our CombinatorialFreeModule Cartesian product superclass::
2714 sage: set_random_seed()
2715 sage: J1 = HadamardEJA(2)
2716 sage: J2 = RealSymmetricEJA(2)
2717 sage: J = cartesian_product([J1,J2])
2718 sage: x,y = J.random_elements(2)
2722 The ability to retrieve the original factors is implemented by our
2723 CombinatorialFreeModule Cartesian product superclass::
2725 sage: J1 = HadamardEJA(2, field=QQ)
2726 sage: J2 = JordanSpinEJA(3, field=QQ)
2727 sage: J = cartesian_product([J1,J2])
2728 sage: J.cartesian_factors()
2729 (Euclidean Jordan algebra of dimension 2 over Rational Field,
2730 Euclidean Jordan algebra of dimension 3 over Rational Field)
2732 You can provide more than two factors::
2734 sage: J1 = HadamardEJA(2)
2735 sage: J2 = JordanSpinEJA(3)
2736 sage: J3 = RealSymmetricEJA(3)
2737 sage: cartesian_product([J1,J2,J3])
2738 Euclidean Jordan algebra of dimension 2 over Algebraic Real
2739 Field (+) Euclidean Jordan algebra of dimension 3 over Algebraic
2740 Real Field (+) Euclidean Jordan algebra of dimension 6 over
2741 Algebraic Real Field
2743 Rank is additive on a Cartesian product::
2745 sage: J1 = HadamardEJA(1)
2746 sage: J2 = RealSymmetricEJA(2)
2747 sage: J = cartesian_product([J1,J2])
2748 sage: J1.rank.clear_cache()
2749 sage: J2.rank.clear_cache()
2750 sage: J.rank.clear_cache()
2753 sage: J.rank() == J1.rank() + J2.rank()
2756 The same rank computation works over the rationals, with whatever
2759 sage: J1 = HadamardEJA(1, field=QQ, orthonormalize=False)
2760 sage: J2 = RealSymmetricEJA(2, field=QQ, orthonormalize=False)
2761 sage: J = cartesian_product([J1,J2])
2762 sage: J1.rank.clear_cache()
2763 sage: J2.rank.clear_cache()
2764 sage: J.rank.clear_cache()
2767 sage: J.rank() == J1.rank() + J2.rank()
2772 All factors must share the same base field::
2774 sage: J1 = HadamardEJA(2, field=QQ)
2775 sage: J2 = RealSymmetricEJA(2)
2776 sage: CartesianProductEJA((J1,J2))
2777 Traceback (most recent call last):
2779 ValueError: all factors must share the same base field
2781 The "cached" Jordan and inner products are the componentwise
2784 sage: set_random_seed()
2785 sage: J1 = random_eja()
2786 sage: J2 = random_eja()
2787 sage: J = cartesian_product([J1,J2])
2788 sage: x,y = J.random_elements(2)
2789 sage: x*y == J.cartesian_jordan_product(x,y)
2791 sage: x.inner_product(y) == J.cartesian_inner_product(x,y)
2794 The cached unit element is the same one that would be computed::
2796 sage: set_random_seed() # long time
2797 sage: J1 = random_eja() # long time
2798 sage: J2 = random_eja() # long time
2799 sage: J = cartesian_product([J1,J2]) # long time
2800 sage: actual = J.one() # long time
2801 sage: J.one.clear_cache() # long time
2802 sage: expected = J.one() # long time
2803 sage: actual == expected # long time
2807 def __init__(self
, modules
, **kwargs
):
2808 CombinatorialFreeModule_CartesianProduct
.__init
__(self
,
2811 field
= modules
[0].base_ring()
2812 if not all( J
.base_ring() == field
for J
in modules
):
2813 raise ValueError("all factors must share the same base field")
2815 basis
= tuple( b
.to_vector().column() for b
in self
.basis() )
2817 # Define jordan/inner products that operate on the basis.
2818 def jordan_product(x_mat
,y_mat
):
2819 x
= self
.from_vector(_mat2vec(x_mat
))
2820 y
= self
.from_vector(_mat2vec(y_mat
))
2821 return self
.cartesian_jordan_product(x
,y
).to_vector().column()
2823 def inner_product(x_mat
, y_mat
):
2824 x
= self
.from_vector(_mat2vec(x_mat
))
2825 y
= self
.from_vector(_mat2vec(y_mat
))
2826 return self
.cartesian_inner_product(x
,y
)
2828 # Use whatever category the superclass came up with. Usually
2829 # some join of the EJA and Cartesian product
2830 # categories. There's no need to check the field since it
2831 # already came from an EJA.
2833 # If you want the basis to be orthonormalized, orthonormalize
2835 FiniteDimensionalEJA
.__init
__(self
,
2840 orthonormalize
=False,
2841 cartesian_product
=True,
2845 ones
= tuple(J
.one() for J
in modules
)
2846 self
.one
.set_cache(self
._cartesian
_product
_of
_elements
(ones
))
2847 self
.rank
.set_cache(sum(J
.rank() for J
in modules
))
2849 # Now that everything else is ready, we clobber our computed
2850 # matrix basis with the "correct" one consisting of ordered
2851 # tuples. Since we didn't orthonormalize our basis, we can
2852 # create these from the basis that was handed to us; that is,
2853 # we don't need to use the one that the earlier __init__()
2854 # method came up with.
2855 m
= len(self
.cartesian_factors())
2856 MS
= self
.matrix_space()
2857 self
._matrix
_basis
= tuple(
2858 MS(tuple( self
.cartesian_projection(i
)(b
).to_matrix()
2859 for i
in range(m
) ))
2860 for b
in self
.basis()
2863 def matrix_space(self
):
2865 Return the space that our matrix basis lives in as a Cartesian
2870 sage: from mjo.eja.eja_algebra import (HadamardEJA,
2871 ....: RealSymmetricEJA)
2875 sage: J1 = HadamardEJA(1)
2876 sage: J2 = RealSymmetricEJA(2)
2877 sage: J = cartesian_product([J1,J2])
2878 sage: J.matrix_space()
2879 The Cartesian product of (Full MatrixSpace of 1 by 1 dense
2880 matrices over Algebraic Real Field, Full MatrixSpace of 2
2881 by 2 dense matrices over Algebraic Real Field)
2884 from sage
.categories
.cartesian_product
import cartesian_product
2885 return cartesian_product( [J
.matrix_space()
2886 for J
in self
.cartesian_factors()] )
2889 def cartesian_projection(self
, i
):
2893 sage: from mjo.eja.eja_algebra import (random_eja,
2894 ....: JordanSpinEJA,
2896 ....: RealSymmetricEJA,
2897 ....: ComplexHermitianEJA)
2901 The projection morphisms are Euclidean Jordan algebra
2904 sage: J1 = HadamardEJA(2)
2905 sage: J2 = RealSymmetricEJA(2)
2906 sage: J = cartesian_product([J1,J2])
2907 sage: J.cartesian_projection(0)
2908 Linear operator between finite-dimensional Euclidean Jordan
2909 algebras represented by the matrix:
2912 Domain: Euclidean Jordan algebra of dimension 2 over Algebraic
2913 Real Field (+) Euclidean Jordan algebra of dimension 3 over
2914 Algebraic Real Field
2915 Codomain: Euclidean Jordan algebra of dimension 2 over Algebraic
2917 sage: J.cartesian_projection(1)
2918 Linear operator between finite-dimensional Euclidean Jordan
2919 algebras represented by the matrix:
2923 Domain: Euclidean Jordan algebra of dimension 2 over Algebraic
2924 Real Field (+) Euclidean Jordan algebra of dimension 3 over
2925 Algebraic Real Field
2926 Codomain: Euclidean Jordan algebra of dimension 3 over Algebraic
2929 The projections work the way you'd expect on the vector
2930 representation of an element::
2932 sage: J1 = JordanSpinEJA(2)
2933 sage: J2 = ComplexHermitianEJA(2)
2934 sage: J = cartesian_product([J1,J2])
2935 sage: pi_left = J.cartesian_projection(0)
2936 sage: pi_right = J.cartesian_projection(1)
2937 sage: pi_left(J.one()).to_vector()
2939 sage: pi_right(J.one()).to_vector()
2941 sage: J.one().to_vector()
2946 The answer never changes::
2948 sage: set_random_seed()
2949 sage: J1 = random_eja()
2950 sage: J2 = random_eja()
2951 sage: J = cartesian_product([J1,J2])
2952 sage: P0 = J.cartesian_projection(0)
2953 sage: P1 = J.cartesian_projection(0)
2958 Ji
= self
.cartesian_factors()[i
]
2959 # Requires the fix on Trac 31421/31422 to work!
2960 Pi
= super().cartesian_projection(i
)
2961 return FiniteDimensionalEJAOperator(self
,Ji
,Pi
.matrix())
2964 def cartesian_embedding(self
, i
):
2968 sage: from mjo.eja.eja_algebra import (random_eja,
2969 ....: JordanSpinEJA,
2971 ....: RealSymmetricEJA)
2975 The embedding morphisms are Euclidean Jordan algebra
2978 sage: J1 = HadamardEJA(2)
2979 sage: J2 = RealSymmetricEJA(2)
2980 sage: J = cartesian_product([J1,J2])
2981 sage: J.cartesian_embedding(0)
2982 Linear operator between finite-dimensional Euclidean Jordan
2983 algebras represented by the matrix:
2989 Domain: Euclidean Jordan algebra of dimension 2 over
2990 Algebraic Real Field
2991 Codomain: Euclidean Jordan algebra of dimension 2 over
2992 Algebraic Real Field (+) Euclidean Jordan algebra of
2993 dimension 3 over Algebraic Real Field
2994 sage: J.cartesian_embedding(1)
2995 Linear operator between finite-dimensional Euclidean Jordan
2996 algebras represented by the matrix:
3002 Domain: Euclidean Jordan algebra of dimension 3 over
3003 Algebraic Real Field
3004 Codomain: Euclidean Jordan algebra of dimension 2 over
3005 Algebraic Real Field (+) Euclidean Jordan algebra of
3006 dimension 3 over Algebraic Real Field
3008 The embeddings work the way you'd expect on the vector
3009 representation of an element::
3011 sage: J1 = JordanSpinEJA(3)
3012 sage: J2 = RealSymmetricEJA(2)
3013 sage: J = cartesian_product([J1,J2])
3014 sage: iota_left = J.cartesian_embedding(0)
3015 sage: iota_right = J.cartesian_embedding(1)
3016 sage: iota_left(J1.zero()) == J.zero()
3018 sage: iota_right(J2.zero()) == J.zero()
3020 sage: J1.one().to_vector()
3022 sage: iota_left(J1.one()).to_vector()
3024 sage: J2.one().to_vector()
3026 sage: iota_right(J2.one()).to_vector()
3028 sage: J.one().to_vector()
3033 The answer never changes::
3035 sage: set_random_seed()
3036 sage: J1 = random_eja()
3037 sage: J2 = random_eja()
3038 sage: J = cartesian_product([J1,J2])
3039 sage: E0 = J.cartesian_embedding(0)
3040 sage: E1 = J.cartesian_embedding(0)
3044 Composing a projection with the corresponding inclusion should
3045 produce the identity map, and mismatching them should produce
3048 sage: set_random_seed()
3049 sage: J1 = random_eja()
3050 sage: J2 = random_eja()
3051 sage: J = cartesian_product([J1,J2])
3052 sage: iota_left = J.cartesian_embedding(0)
3053 sage: iota_right = J.cartesian_embedding(1)
3054 sage: pi_left = J.cartesian_projection(0)
3055 sage: pi_right = J.cartesian_projection(1)
3056 sage: pi_left*iota_left == J1.one().operator()
3058 sage: pi_right*iota_right == J2.one().operator()
3060 sage: (pi_left*iota_right).is_zero()
3062 sage: (pi_right*iota_left).is_zero()
3066 Ji
= self
.cartesian_factors()[i
]
3067 # Requires the fix on Trac 31421/31422 to work!
3068 Ei
= super().cartesian_embedding(i
)
3069 return FiniteDimensionalEJAOperator(Ji
,self
,Ei
.matrix())
3072 def cartesian_jordan_product(self
, x
, y
):
3074 The componentwise Jordan product.
3076 We project ``x`` and ``y`` onto our factors, and add up the
3077 Jordan products from the subalgebras. This may still be useful
3078 after (if) the default Jordan product in the Cartesian product
3079 algebra is overridden.
3083 sage: from mjo.eja.eja_algebra import (HadamardEJA,
3084 ....: JordanSpinEJA)
3088 sage: J1 = HadamardEJA(3)
3089 sage: J2 = JordanSpinEJA(3)
3090 sage: J = cartesian_product([J1,J2])
3091 sage: x1 = J1.from_vector(vector(QQ,(1,2,1)))
3092 sage: y1 = J1.from_vector(vector(QQ,(1,0,2)))
3093 sage: x2 = J2.from_vector(vector(QQ,(1,2,3)))
3094 sage: y2 = J2.from_vector(vector(QQ,(1,1,1)))
3095 sage: z1 = J.from_vector(vector(QQ,(1,2,1,1,2,3)))
3096 sage: z2 = J.from_vector(vector(QQ,(1,0,2,1,1,1)))
3097 sage: (x1*y1).to_vector()
3099 sage: (x2*y2).to_vector()
3101 sage: J.cartesian_jordan_product(z1,z2).to_vector()
3105 m
= len(self
.cartesian_factors())
3106 projections
= ( self
.cartesian_projection(i
) for i
in range(m
) )
3107 products
= ( P(x
)*P(y
) for P
in projections
)
3108 return self
._cartesian
_product
_of
_elements
(tuple(products
))
3110 def cartesian_inner_product(self
, x
, y
):
3112 The standard componentwise Cartesian inner-product.
3114 We project ``x`` and ``y`` onto our factors, and add up the
3115 inner-products from the subalgebras. This may still be useful
3116 after (if) the default inner product in the Cartesian product
3117 algebra is overridden.
3121 sage: from mjo.eja.eja_algebra import (HadamardEJA,
3122 ....: QuaternionHermitianEJA)
3126 sage: J1 = HadamardEJA(3,field=QQ)
3127 sage: J2 = QuaternionHermitianEJA(2,field=QQ,orthonormalize=False)
3128 sage: J = cartesian_product([J1,J2])
3133 sage: x1.inner_product(x2)
3135 sage: y1.inner_product(y2)
3137 sage: z1 = J._cartesian_product_of_elements((x1,y1))
3138 sage: z2 = J._cartesian_product_of_elements((x2,y2))
3139 sage: J.cartesian_inner_product(z1,z2)
3143 m
= len(self
.cartesian_factors())
3144 projections
= ( self
.cartesian_projection(i
) for i
in range(m
) )
3145 return sum( P(x
).inner_product(P(y
)) for P
in projections
)
3148 def _element_constructor_(self
, elt
):
3150 Construct an element of this algebra from an ordered tuple.
3152 We just apply the element constructor from each of our factors
3153 to the corresponding component of the tuple, and package up
3158 sage: from mjo.eja.eja_algebra import (HadamardEJA,
3159 ....: RealSymmetricEJA)
3163 sage: J1 = HadamardEJA(3)
3164 sage: J2 = RealSymmetricEJA(2)
3165 sage: J = cartesian_product([J1,J2])
3166 sage: J( (J1.matrix_basis()[1], J2.matrix_basis()[2]) )
3169 m
= len(self
.cartesian_factors())
3171 z
= tuple( self
.cartesian_factors()[i
](elt
[i
]) for i
in range(m
) )
3172 return self
._cartesian
_product
_of
_elements
(z
)
3174 raise ValueError("not an element of this algebra")
3176 Element
= CartesianProductEJAElement
3179 FiniteDimensionalEJA
.CartesianProduct
= CartesianProductEJA
3180 random_eja
= ConcreteEJA
.random_instance