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 sage: J = HadamardEJA(5)
757 e0 + e1 + e2 + e3 + e4
761 The identity element acts like the identity::
763 sage: set_random_seed()
764 sage: J = random_eja()
765 sage: x = J.random_element()
766 sage: J.one()*x == x and x*J.one() == x
769 The matrix of the unit element's operator is the identity::
771 sage: set_random_seed()
772 sage: J = random_eja()
773 sage: actual = J.one().operator().matrix()
774 sage: expected = matrix.identity(J.base_ring(), J.dimension())
775 sage: actual == expected
778 Ensure that the cached unit element (often precomputed by
779 hand) agrees with the computed one::
781 sage: set_random_seed()
782 sage: J = random_eja()
783 sage: cached = J.one()
784 sage: J.one.clear_cache()
785 sage: J.one() == cached
789 # We can brute-force compute the matrices of the operators
790 # that correspond to the basis elements of this algebra.
791 # If some linear combination of those basis elements is the
792 # algebra identity, then the same linear combination of
793 # their matrices has to be the identity matrix.
795 # Of course, matrices aren't vectors in sage, so we have to
796 # appeal to the "long vectors" isometry.
797 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
799 # Now we use basic linear algebra to find the coefficients,
800 # of the matrices-as-vectors-linear-combination, which should
801 # work for the original algebra basis too.
802 A
= matrix(self
.base_ring(), oper_vecs
)
804 # We used the isometry on the left-hand side already, but we
805 # still need to do it for the right-hand side. Recall that we
806 # wanted something that summed to the identity matrix.
807 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
809 # Now if there's an identity element in the algebra, this
810 # should work. We solve on the left to avoid having to
811 # transpose the matrix "A".
812 return self
.from_vector(A
.solve_left(b
))
815 def peirce_decomposition(self
, c
):
817 The Peirce decomposition of this algebra relative to the
820 In the future, this can be extended to a complete system of
821 orthogonal idempotents.
825 - ``c`` -- an idempotent of this algebra.
829 A triple (J0, J5, J1) containing two subalgebras and one subspace
832 - ``J0`` -- the algebra on the eigenspace of ``c.operator()``
833 corresponding to the eigenvalue zero.
835 - ``J5`` -- the eigenspace (NOT a subalgebra) of ``c.operator()``
836 corresponding to the eigenvalue one-half.
838 - ``J1`` -- the algebra on the eigenspace of ``c.operator()``
839 corresponding to the eigenvalue one.
841 These are the only possible eigenspaces for that operator, and this
842 algebra is a direct sum of them. The spaces ``J0`` and ``J1`` are
843 orthogonal, and are subalgebras of this algebra with the appropriate
848 sage: from mjo.eja.eja_algebra import random_eja, RealSymmetricEJA
852 The canonical example comes from the symmetric matrices, which
853 decompose into diagonal and off-diagonal parts::
855 sage: J = RealSymmetricEJA(3)
856 sage: C = matrix(QQ, [ [1,0,0],
860 sage: J0,J5,J1 = J.peirce_decomposition(c)
862 Euclidean Jordan algebra of dimension 1...
864 Vector space of degree 6 and dimension 2...
866 Euclidean Jordan algebra of dimension 3...
867 sage: J0.one().to_matrix()
871 sage: orig_df = AA.options.display_format
872 sage: AA.options.display_format = 'radical'
873 sage: J.from_vector(J5.basis()[0]).to_matrix()
877 sage: J.from_vector(J5.basis()[1]).to_matrix()
881 sage: AA.options.display_format = orig_df
882 sage: J1.one().to_matrix()
889 Every algebra decomposes trivially with respect to its identity
892 sage: set_random_seed()
893 sage: J = random_eja()
894 sage: J0,J5,J1 = J.peirce_decomposition(J.one())
895 sage: J0.dimension() == 0 and J5.dimension() == 0
897 sage: J1.superalgebra() == J and J1.dimension() == J.dimension()
900 The decomposition is into eigenspaces, and its components are
901 therefore necessarily orthogonal. Moreover, the identity
902 elements in the two subalgebras are the projections onto their
903 respective subspaces of the superalgebra's identity element::
905 sage: set_random_seed()
906 sage: J = random_eja()
907 sage: x = J.random_element()
908 sage: if not J.is_trivial():
909 ....: while x.is_nilpotent():
910 ....: x = J.random_element()
911 sage: c = x.subalgebra_idempotent()
912 sage: J0,J5,J1 = J.peirce_decomposition(c)
914 sage: for (w,y,z) in zip(J0.basis(), J5.basis(), J1.basis()):
915 ....: w = w.superalgebra_element()
916 ....: y = J.from_vector(y)
917 ....: z = z.superalgebra_element()
918 ....: ipsum += w.inner_product(y).abs()
919 ....: ipsum += w.inner_product(z).abs()
920 ....: ipsum += y.inner_product(z).abs()
923 sage: J1(c) == J1.one()
925 sage: J0(J.one() - c) == J0.one()
929 if not c
.is_idempotent():
930 raise ValueError("element is not idempotent: %s" % c
)
932 from mjo
.eja
.eja_subalgebra
import FiniteDimensionalEJASubalgebra
934 # Default these to what they should be if they turn out to be
935 # trivial, because eigenspaces_left() won't return eigenvalues
936 # corresponding to trivial spaces (e.g. it returns only the
937 # eigenspace corresponding to lambda=1 if you take the
938 # decomposition relative to the identity element).
939 trivial
= FiniteDimensionalEJASubalgebra(self
, ())
940 J0
= trivial
# eigenvalue zero
941 J5
= VectorSpace(self
.base_ring(), 0) # eigenvalue one-half
942 J1
= trivial
# eigenvalue one
944 for (eigval
, eigspace
) in c
.operator().matrix().right_eigenspaces():
945 if eigval
== ~
(self
.base_ring()(2)):
948 gens
= tuple( self
.from_vector(b
) for b
in eigspace
.basis() )
949 subalg
= FiniteDimensionalEJASubalgebra(self
,
957 raise ValueError("unexpected eigenvalue: %s" % eigval
)
962 def random_element(self
, thorough
=False):
964 Return a random element of this algebra.
966 Our algebra superclass method only returns a linear
967 combination of at most two basis elements. We instead
968 want the vector space "random element" method that
969 returns a more diverse selection.
973 - ``thorough`` -- (boolean; default False) whether or not we
974 should generate irrational coefficients for the random
975 element when our base ring is irrational; this slows the
976 algebra operations to a crawl, but any truly random method
980 # For a general base ring... maybe we can trust this to do the
981 # right thing? Unlikely, but.
982 V
= self
.vector_space()
983 v
= V
.random_element()
985 if self
.base_ring() is AA
:
986 # The "random element" method of the algebraic reals is
987 # stupid at the moment, and only returns integers between
988 # -2 and 2, inclusive:
990 # https://trac.sagemath.org/ticket/30875
992 # Instead, we implement our own "random vector" method,
993 # and then coerce that into the algebra. We use the vector
994 # space degree here instead of the dimension because a
995 # subalgebra could (for example) be spanned by only two
996 # vectors, each with five coordinates. We need to
997 # generate all five coordinates.
999 v
*= QQbar
.random_element().real()
1001 v
*= QQ
.random_element()
1003 return self
.from_vector(V
.coordinate_vector(v
))
1005 def random_elements(self
, count
, thorough
=False):
1007 Return ``count`` random elements as a tuple.
1011 - ``thorough`` -- (boolean; default False) whether or not we
1012 should generate irrational coefficients for the random
1013 elements when our base ring is irrational; this slows the
1014 algebra operations to a crawl, but any truly random method
1019 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1023 sage: J = JordanSpinEJA(3)
1024 sage: x,y,z = J.random_elements(3)
1025 sage: all( [ x in J, y in J, z in J ])
1027 sage: len( J.random_elements(10) ) == 10
1031 return tuple( self
.random_element(thorough
)
1032 for idx
in range(count
) )
1036 def _charpoly_coefficients(self
):
1038 The `r` polynomial coefficients of the "characteristic polynomial
1041 n
= self
.dimension()
1042 R
= self
.coordinate_polynomial_ring()
1044 F
= R
.fraction_field()
1047 # From a result in my book, these are the entries of the
1048 # basis representation of L_x.
1049 return sum( vars[k
]*self
.monomial(k
).operator().matrix()[i
,j
]
1052 L_x
= matrix(F
, n
, n
, L_x_i_j
)
1055 if self
.rank
.is_in_cache():
1057 # There's no need to pad the system with redundant
1058 # columns if we *know* they'll be redundant.
1061 # Compute an extra power in case the rank is equal to
1062 # the dimension (otherwise, we would stop at x^(r-1)).
1063 x_powers
= [ (L_x
**k
)*self
.one().to_vector()
1064 for k
in range(n
+1) ]
1065 A
= matrix
.column(F
, x_powers
[:n
])
1066 AE
= A
.extended_echelon_form()
1073 # The theory says that only the first "r" coefficients are
1074 # nonzero, and they actually live in the original polynomial
1075 # ring and not the fraction field. We negate them because
1076 # in the actual characteristic polynomial, they get moved
1077 # to the other side where x^r lives.
1078 return -A_rref
.solve_right(E
*b
).change_ring(R
)[:r
]
1083 Return the rank of this EJA.
1085 This is a cached method because we know the rank a priori for
1086 all of the algebras we can construct. Thus we can avoid the
1087 expensive ``_charpoly_coefficients()`` call unless we truly
1088 need to compute the whole characteristic polynomial.
1092 sage: from mjo.eja.eja_algebra import (HadamardEJA,
1093 ....: JordanSpinEJA,
1094 ....: RealSymmetricEJA,
1095 ....: ComplexHermitianEJA,
1096 ....: QuaternionHermitianEJA,
1101 The rank of the Jordan spin algebra is always two::
1103 sage: JordanSpinEJA(2).rank()
1105 sage: JordanSpinEJA(3).rank()
1107 sage: JordanSpinEJA(4).rank()
1110 The rank of the `n`-by-`n` Hermitian real, complex, or
1111 quaternion matrices is `n`::
1113 sage: RealSymmetricEJA(4).rank()
1115 sage: ComplexHermitianEJA(3).rank()
1117 sage: QuaternionHermitianEJA(2).rank()
1122 Ensure that every EJA that we know how to construct has a
1123 positive integer rank, unless the algebra is trivial in
1124 which case its rank will be zero::
1126 sage: set_random_seed()
1127 sage: J = random_eja()
1131 sage: r > 0 or (r == 0 and J.is_trivial())
1134 Ensure that computing the rank actually works, since the ranks
1135 of all simple algebras are known and will be cached by default::
1137 sage: set_random_seed() # long time
1138 sage: J = random_eja() # long time
1139 sage: caches = J.rank() # long time
1140 sage: J.rank.clear_cache() # long time
1141 sage: J.rank() == cached # long time
1145 return len(self
._charpoly
_coefficients
())
1148 def vector_space(self
):
1150 Return the vector space that underlies this algebra.
1154 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1158 sage: J = RealSymmetricEJA(2)
1159 sage: J.vector_space()
1160 Vector space of dimension 3 over...
1163 return self
.zero().to_vector().parent().ambient_vector_space()
1166 Element
= FiniteDimensionalEJAElement
1168 class RationalBasisEJA(FiniteDimensionalEJA
):
1170 New class for algebras whose supplied basis elements have all rational entries.
1174 sage: from mjo.eja.eja_algebra import BilinearFormEJA
1178 The supplied basis is orthonormalized by default::
1180 sage: B = matrix(QQ, [[1, 0, 0], [0, 25, -32], [0, -32, 41]])
1181 sage: J = BilinearFormEJA(B)
1182 sage: J.matrix_basis()
1199 # Abuse the check_field parameter to check that the entries of
1200 # out basis (in ambient coordinates) are in the field QQ.
1201 if not all( all(b_i
in QQ
for b_i
in b
.list()) for b
in basis
):
1202 raise TypeError("basis not rational")
1204 self
._rational
_algebra
= None
1206 # There's no point in constructing the extra algebra if this
1207 # one is already rational.
1209 # Note: the same Jordan and inner-products work here,
1210 # because they are necessarily defined with respect to
1211 # ambient coordinates and not any particular basis.
1212 self
._rational
_algebra
= FiniteDimensionalEJA(
1217 orthonormalize
=False,
1221 super().__init
__(basis
,
1225 check_field
=check_field
,
1229 def _charpoly_coefficients(self
):
1233 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
1234 ....: JordanSpinEJA)
1238 The base ring of the resulting polynomial coefficients is what
1239 it should be, and not the rationals (unless the algebra was
1240 already over the rationals)::
1242 sage: J = JordanSpinEJA(3)
1243 sage: J._charpoly_coefficients()
1244 (X1^2 - X2^2 - X3^2, -2*X1)
1245 sage: a0 = J._charpoly_coefficients()[0]
1247 Algebraic Real Field
1248 sage: a0.base_ring()
1249 Algebraic Real Field
1252 if self
._rational
_algebra
is None:
1253 # There's no need to construct *another* algebra over the
1254 # rationals if this one is already over the
1255 # rationals. Likewise, if we never orthonormalized our
1256 # basis, we might as well just use the given one.
1257 return super()._charpoly
_coefficients
()
1259 # Do the computation over the rationals. The answer will be
1260 # the same, because all we've done is a change of basis.
1261 # Then, change back from QQ to our real base ring
1262 a
= ( a_i
.change_ring(self
.base_ring())
1263 for a_i
in self
._rational
_algebra
._charpoly
_coefficients
() )
1265 if self
._deortho
_matrix
is None:
1266 # This can happen if our base ring was, say, AA and we
1267 # chose not to (or didn't need to) orthonormalize. It's
1268 # still faster to do the computations over QQ even if
1269 # the numbers in the boxes stay the same.
1272 # Otherwise, convert the coordinate variables back to the
1273 # deorthonormalized ones.
1274 R
= self
.coordinate_polynomial_ring()
1275 from sage
.modules
.free_module_element
import vector
1276 X
= vector(R
, R
.gens())
1277 BX
= self
._deortho
_matrix
*X
1279 subs_dict
= { X[i]: BX[i] for i in range(len(X)) }
1280 return tuple( a_i
.subs(subs_dict
) for a_i
in a
)
1282 class ConcreteEJA(RationalBasisEJA
):
1284 A class for the Euclidean Jordan algebras that we know by name.
1286 These are the Jordan algebras whose basis, multiplication table,
1287 rank, and so on are known a priori. More to the point, they are
1288 the Euclidean Jordan algebras for which we are able to conjure up
1289 a "random instance."
1293 sage: from mjo.eja.eja_algebra import ConcreteEJA
1297 Our basis is normalized with respect to the algebra's inner
1298 product, unless we specify otherwise::
1300 sage: set_random_seed()
1301 sage: J = ConcreteEJA.random_instance()
1302 sage: all( b.norm() == 1 for b in J.gens() )
1305 Since our basis is orthonormal with respect to the algebra's inner
1306 product, and since we know that this algebra is an EJA, any
1307 left-multiplication operator's matrix will be symmetric because
1308 natural->EJA basis representation is an isometry and within the
1309 EJA the operator is self-adjoint by the Jordan axiom::
1311 sage: set_random_seed()
1312 sage: J = ConcreteEJA.random_instance()
1313 sage: x = J.random_element()
1314 sage: x.operator().is_self_adjoint()
1319 def _max_random_instance_size():
1321 Return an integer "size" that is an upper bound on the size of
1322 this algebra when it is used in a random test
1323 case. Unfortunately, the term "size" is ambiguous -- when
1324 dealing with `R^n` under either the Hadamard or Jordan spin
1325 product, the "size" refers to the dimension `n`. When dealing
1326 with a matrix algebra (real symmetric or complex/quaternion
1327 Hermitian), it refers to the size of the matrix, which is far
1328 less than the dimension of the underlying vector space.
1330 This method must be implemented in each subclass.
1332 raise NotImplementedError
1335 def random_instance(cls
, *args
, **kwargs
):
1337 Return a random instance of this type of algebra.
1339 This method should be implemented in each subclass.
1341 from sage
.misc
.prandom
import choice
1342 eja_class
= choice(cls
.__subclasses
__())
1344 # These all bubble up to the RationalBasisEJA superclass
1345 # constructor, so any (kw)args valid there are also valid
1347 return eja_class
.random_instance(*args
, **kwargs
)
1352 def dimension_over_reals():
1354 The dimension of this matrix's base ring over the reals.
1356 The reals are dimension one over themselves, obviously; that's
1357 just `\mathbb{R}^{1}`. Likewise, the complex numbers `a + bi`
1358 have dimension two. Finally, the quaternions have dimension
1359 four over the reals.
1361 This is used to determine the size of the matrix returned from
1362 :meth:`real_embed`, among other things.
1364 raise NotImplementedError
1367 def real_embed(cls
,M
):
1369 Embed the matrix ``M`` into a space of real matrices.
1371 The matrix ``M`` can have entries in any field at the moment:
1372 the real numbers, complex numbers, or quaternions. And although
1373 they are not a field, we can probably support octonions at some
1374 point, too. This function returns a real matrix that "acts like"
1375 the original with respect to matrix multiplication; i.e.
1377 real_embed(M*N) = real_embed(M)*real_embed(N)
1380 if M
.ncols() != M
.nrows():
1381 raise ValueError("the matrix 'M' must be square")
1386 def real_unembed(cls
,M
):
1388 The inverse of :meth:`real_embed`.
1390 if M
.ncols() != M
.nrows():
1391 raise ValueError("the matrix 'M' must be square")
1392 if not ZZ(M
.nrows()).mod(cls
.dimension_over_reals()).is_zero():
1393 raise ValueError("the matrix 'M' must be a real embedding")
1397 def jordan_product(X
,Y
):
1398 return (X
*Y
+ Y
*X
)/2
1401 def trace_inner_product(cls
,X
,Y
):
1403 Compute the trace inner-product of two real-embeddings.
1407 sage: from mjo.eja.eja_algebra import (RealSymmetricEJA,
1408 ....: ComplexHermitianEJA,
1409 ....: QuaternionHermitianEJA)
1413 This gives the same answer as it would if we computed the trace
1414 from the unembedded (original) matrices::
1416 sage: set_random_seed()
1417 sage: J = RealSymmetricEJA.random_instance()
1418 sage: x,y = J.random_elements(2)
1419 sage: Xe = x.to_matrix()
1420 sage: Ye = y.to_matrix()
1421 sage: X = J.real_unembed(Xe)
1422 sage: Y = J.real_unembed(Ye)
1423 sage: expected = (X*Y).trace()
1424 sage: actual = J.trace_inner_product(Xe,Ye)
1425 sage: actual == expected
1430 sage: set_random_seed()
1431 sage: J = ComplexHermitianEJA.random_instance()
1432 sage: x,y = J.random_elements(2)
1433 sage: Xe = x.to_matrix()
1434 sage: Ye = y.to_matrix()
1435 sage: X = J.real_unembed(Xe)
1436 sage: Y = J.real_unembed(Ye)
1437 sage: expected = (X*Y).trace().real()
1438 sage: actual = J.trace_inner_product(Xe,Ye)
1439 sage: actual == expected
1444 sage: set_random_seed()
1445 sage: J = QuaternionHermitianEJA.random_instance()
1446 sage: x,y = J.random_elements(2)
1447 sage: Xe = x.to_matrix()
1448 sage: Ye = y.to_matrix()
1449 sage: X = J.real_unembed(Xe)
1450 sage: Y = J.real_unembed(Ye)
1451 sage: expected = (X*Y).trace().coefficient_tuple()[0]
1452 sage: actual = J.trace_inner_product(Xe,Ye)
1453 sage: actual == expected
1457 Xu
= cls
.real_unembed(X
)
1458 Yu
= cls
.real_unembed(Y
)
1459 tr
= (Xu
*Yu
).trace()
1462 # Works in QQ, AA, RDF, et cetera.
1464 except AttributeError:
1465 # A quaternion doesn't have a real() method, but does
1466 # have coefficient_tuple() method that returns the
1467 # coefficients of 1, i, j, and k -- in that order.
1468 return tr
.coefficient_tuple()[0]
1471 class RealMatrixEJA(MatrixEJA
):
1473 def dimension_over_reals():
1477 class RealSymmetricEJA(ConcreteEJA
, RealMatrixEJA
):
1479 The rank-n simple EJA consisting of real symmetric n-by-n
1480 matrices, the usual symmetric Jordan product, and the trace inner
1481 product. It has dimension `(n^2 + n)/2` over the reals.
1485 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1489 sage: J = RealSymmetricEJA(2)
1490 sage: e0, e1, e2 = J.gens()
1498 In theory, our "field" can be any subfield of the reals::
1500 sage: RealSymmetricEJA(2, field=RDF)
1501 Euclidean Jordan algebra of dimension 3 over Real Double Field
1502 sage: RealSymmetricEJA(2, field=RR)
1503 Euclidean Jordan algebra of dimension 3 over Real Field with
1504 53 bits of precision
1508 The dimension of this algebra is `(n^2 + n) / 2`::
1510 sage: set_random_seed()
1511 sage: n_max = RealSymmetricEJA._max_random_instance_size()
1512 sage: n = ZZ.random_element(1, n_max)
1513 sage: J = RealSymmetricEJA(n)
1514 sage: J.dimension() == (n^2 + n)/2
1517 The Jordan multiplication is what we think it is::
1519 sage: set_random_seed()
1520 sage: J = RealSymmetricEJA.random_instance()
1521 sage: x,y = J.random_elements(2)
1522 sage: actual = (x*y).to_matrix()
1523 sage: X = x.to_matrix()
1524 sage: Y = y.to_matrix()
1525 sage: expected = (X*Y + Y*X)/2
1526 sage: actual == expected
1528 sage: J(expected) == x*y
1531 We can change the generator prefix::
1533 sage: RealSymmetricEJA(3, prefix='q').gens()
1534 (q0, q1, q2, q3, q4, q5)
1536 We can construct the (trivial) algebra of rank zero::
1538 sage: RealSymmetricEJA(0)
1539 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1543 def _denormalized_basis(cls
, n
):
1545 Return a basis for the space of real symmetric n-by-n matrices.
1549 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1553 sage: set_random_seed()
1554 sage: n = ZZ.random_element(1,5)
1555 sage: B = RealSymmetricEJA._denormalized_basis(n)
1556 sage: all( M.is_symmetric() for M in B)
1560 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1564 for j
in range(i
+1):
1565 Eij
= matrix(ZZ
, n
, lambda k
,l
: k
==i
and l
==j
)
1569 Sij
= Eij
+ Eij
.transpose()
1575 def _max_random_instance_size():
1576 return 4 # Dimension 10
1579 def random_instance(cls
, **kwargs
):
1581 Return a random instance of this type of algebra.
1583 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1584 return cls(n
, **kwargs
)
1586 def __init__(self
, n
, **kwargs
):
1587 # We know this is a valid EJA, but will double-check
1588 # if the user passes check_axioms=True.
1589 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
1591 super(RealSymmetricEJA
, self
).__init
__(self
._denormalized
_basis
(n
),
1592 self
.jordan_product
,
1593 self
.trace_inner_product
,
1596 # TODO: this could be factored out somehow, but is left here
1597 # because the MatrixEJA is not presently a subclass of the
1598 # FDEJA class that defines rank() and one().
1599 self
.rank
.set_cache(n
)
1600 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
1601 self
.one
.set_cache(self(idV
))
1605 class ComplexMatrixEJA(MatrixEJA
):
1607 def dimension_over_reals():
1611 def real_embed(cls
,M
):
1613 Embed the n-by-n complex matrix ``M`` into the space of real
1614 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1615 bi` to the block matrix ``[[a,b],[-b,a]]``.
1619 sage: from mjo.eja.eja_algebra import ComplexMatrixEJA
1623 sage: F = QuadraticField(-1, 'I')
1624 sage: x1 = F(4 - 2*i)
1625 sage: x2 = F(1 + 2*i)
1628 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1629 sage: ComplexMatrixEJA.real_embed(M)
1638 Embedding is a homomorphism (isomorphism, in fact)::
1640 sage: set_random_seed()
1641 sage: n = ZZ.random_element(3)
1642 sage: F = QuadraticField(-1, 'I')
1643 sage: X = random_matrix(F, n)
1644 sage: Y = random_matrix(F, n)
1645 sage: Xe = ComplexMatrixEJA.real_embed(X)
1646 sage: Ye = ComplexMatrixEJA.real_embed(Y)
1647 sage: XYe = ComplexMatrixEJA.real_embed(X*Y)
1652 super(ComplexMatrixEJA
,cls
).real_embed(M
)
1655 # We don't need any adjoined elements...
1656 field
= M
.base_ring().base_ring()
1660 a
= z
.list()[0] # real part, I guess
1661 b
= z
.list()[1] # imag part, I guess
1662 blocks
.append(matrix(field
, 2, [[a
,b
],[-b
,a
]]))
1664 return matrix
.block(field
, n
, blocks
)
1668 def real_unembed(cls
,M
):
1670 The inverse of _embed_complex_matrix().
1674 sage: from mjo.eja.eja_algebra import ComplexMatrixEJA
1678 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1679 ....: [-2, 1, -4, 3],
1680 ....: [ 9, 10, 11, 12],
1681 ....: [-10, 9, -12, 11] ])
1682 sage: ComplexMatrixEJA.real_unembed(A)
1684 [ 10*I + 9 12*I + 11]
1688 Unembedding is the inverse of embedding::
1690 sage: set_random_seed()
1691 sage: F = QuadraticField(-1, 'I')
1692 sage: M = random_matrix(F, 3)
1693 sage: Me = ComplexMatrixEJA.real_embed(M)
1694 sage: ComplexMatrixEJA.real_unembed(Me) == M
1698 super(ComplexMatrixEJA
,cls
).real_unembed(M
)
1700 d
= cls
.dimension_over_reals()
1702 # If "M" was normalized, its base ring might have roots
1703 # adjoined and they can stick around after unembedding.
1704 field
= M
.base_ring()
1705 R
= PolynomialRing(field
, 'z')
1708 # Sage doesn't know how to adjoin the complex "i" (the root of
1709 # x^2 + 1) to a field in a general way. Here, we just enumerate
1710 # all of the cases that I have cared to support so far.
1712 # Sage doesn't know how to embed AA into QQbar, i.e. how
1713 # to adjoin sqrt(-1) to AA.
1715 elif not field
.is_exact():
1717 F
= field
.complex_field()
1719 # Works for QQ and... maybe some other fields.
1720 F
= field
.extension(z
**2 + 1, 'I', embedding
=CLF(-1).sqrt())
1723 # Go top-left to bottom-right (reading order), converting every
1724 # 2-by-2 block we see to a single complex element.
1726 for k
in range(n
/d
):
1727 for j
in range(n
/d
):
1728 submat
= M
[d
*k
:d
*k
+d
,d
*j
:d
*j
+d
]
1729 if submat
[0,0] != submat
[1,1]:
1730 raise ValueError('bad on-diagonal submatrix')
1731 if submat
[0,1] != -submat
[1,0]:
1732 raise ValueError('bad off-diagonal submatrix')
1733 z
= submat
[0,0] + submat
[0,1]*i
1736 return matrix(F
, n
/d
, elements
)
1739 class ComplexHermitianEJA(ConcreteEJA
, ComplexMatrixEJA
):
1741 The rank-n simple EJA consisting of complex Hermitian n-by-n
1742 matrices over the real numbers, the usual symmetric Jordan product,
1743 and the real-part-of-trace inner product. It has dimension `n^2` over
1748 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1752 In theory, our "field" can be any subfield of the reals::
1754 sage: ComplexHermitianEJA(2, field=RDF)
1755 Euclidean Jordan algebra of dimension 4 over Real Double Field
1756 sage: ComplexHermitianEJA(2, field=RR)
1757 Euclidean Jordan algebra of dimension 4 over Real Field with
1758 53 bits of precision
1762 The dimension of this algebra is `n^2`::
1764 sage: set_random_seed()
1765 sage: n_max = ComplexHermitianEJA._max_random_instance_size()
1766 sage: n = ZZ.random_element(1, n_max)
1767 sage: J = ComplexHermitianEJA(n)
1768 sage: J.dimension() == n^2
1771 The Jordan multiplication is what we think it is::
1773 sage: set_random_seed()
1774 sage: J = ComplexHermitianEJA.random_instance()
1775 sage: x,y = J.random_elements(2)
1776 sage: actual = (x*y).to_matrix()
1777 sage: X = x.to_matrix()
1778 sage: Y = y.to_matrix()
1779 sage: expected = (X*Y + Y*X)/2
1780 sage: actual == expected
1782 sage: J(expected) == x*y
1785 We can change the generator prefix::
1787 sage: ComplexHermitianEJA(2, prefix='z').gens()
1790 We can construct the (trivial) algebra of rank zero::
1792 sage: ComplexHermitianEJA(0)
1793 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1798 def _denormalized_basis(cls
, n
):
1800 Returns a basis for the space of complex Hermitian n-by-n matrices.
1802 Why do we embed these? Basically, because all of numerical linear
1803 algebra assumes that you're working with vectors consisting of `n`
1804 entries from a field and scalars from the same field. There's no way
1805 to tell SageMath that (for example) the vectors contain complex
1806 numbers, while the scalar field is real.
1810 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1814 sage: set_random_seed()
1815 sage: n = ZZ.random_element(1,5)
1816 sage: field = QuadraticField(2, 'sqrt2')
1817 sage: B = ComplexHermitianEJA._denormalized_basis(n)
1818 sage: all( M.is_symmetric() for M in B)
1823 R
= PolynomialRing(field
, 'z')
1825 F
= field
.extension(z
**2 + 1, 'I')
1828 # This is like the symmetric case, but we need to be careful:
1830 # * We want conjugate-symmetry, not just symmetry.
1831 # * The diagonal will (as a result) be real.
1835 for j
in range(i
+1):
1836 Eij
= matrix(F
, n
, lambda k
,l
: k
==i
and l
==j
)
1838 Sij
= cls
.real_embed(Eij
)
1841 # The second one has a minus because it's conjugated.
1842 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
1844 Sij_imag
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
1847 # Since we embedded these, we can drop back to the "field" that we
1848 # started with instead of the complex extension "F".
1849 return tuple( s
.change_ring(field
) for s
in S
)
1852 def __init__(self
, n
, **kwargs
):
1853 # We know this is a valid EJA, but will double-check
1854 # if the user passes check_axioms=True.
1855 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
1857 super(ComplexHermitianEJA
, self
).__init
__(self
._denormalized
_basis
(n
),
1858 self
.jordan_product
,
1859 self
.trace_inner_product
,
1861 # TODO: this could be factored out somehow, but is left here
1862 # because the MatrixEJA is not presently a subclass of the
1863 # FDEJA class that defines rank() and one().
1864 self
.rank
.set_cache(n
)
1865 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
1866 self
.one
.set_cache(self(idV
))
1869 def _max_random_instance_size():
1870 return 3 # Dimension 9
1873 def random_instance(cls
, **kwargs
):
1875 Return a random instance of this type of algebra.
1877 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1878 return cls(n
, **kwargs
)
1880 class QuaternionMatrixEJA(MatrixEJA
):
1882 def dimension_over_reals():
1886 def real_embed(cls
,M
):
1888 Embed the n-by-n quaternion matrix ``M`` into the space of real
1889 matrices of size 4n-by-4n by first sending each quaternion entry `z
1890 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
1891 c+di],[-c + di, a-bi]]`, and then embedding those into a real
1896 sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA
1900 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1901 sage: i,j,k = Q.gens()
1902 sage: x = 1 + 2*i + 3*j + 4*k
1903 sage: M = matrix(Q, 1, [[x]])
1904 sage: QuaternionMatrixEJA.real_embed(M)
1910 Embedding is a homomorphism (isomorphism, in fact)::
1912 sage: set_random_seed()
1913 sage: n = ZZ.random_element(2)
1914 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1915 sage: X = random_matrix(Q, n)
1916 sage: Y = random_matrix(Q, n)
1917 sage: Xe = QuaternionMatrixEJA.real_embed(X)
1918 sage: Ye = QuaternionMatrixEJA.real_embed(Y)
1919 sage: XYe = QuaternionMatrixEJA.real_embed(X*Y)
1924 super(QuaternionMatrixEJA
,cls
).real_embed(M
)
1925 quaternions
= M
.base_ring()
1928 F
= QuadraticField(-1, 'I')
1933 t
= z
.coefficient_tuple()
1938 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
1939 [-c
+ d
*i
, a
- b
*i
]])
1940 realM
= ComplexMatrixEJA
.real_embed(cplxM
)
1941 blocks
.append(realM
)
1943 # We should have real entries by now, so use the realest field
1944 # we've got for the return value.
1945 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
1950 def real_unembed(cls
,M
):
1952 The inverse of _embed_quaternion_matrix().
1956 sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA
1960 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
1961 ....: [-2, 1, -4, 3],
1962 ....: [-3, 4, 1, -2],
1963 ....: [-4, -3, 2, 1]])
1964 sage: QuaternionMatrixEJA.real_unembed(M)
1965 [1 + 2*i + 3*j + 4*k]
1969 Unembedding is the inverse of embedding::
1971 sage: set_random_seed()
1972 sage: Q = QuaternionAlgebra(QQ, -1, -1)
1973 sage: M = random_matrix(Q, 3)
1974 sage: Me = QuaternionMatrixEJA.real_embed(M)
1975 sage: QuaternionMatrixEJA.real_unembed(Me) == M
1979 super(QuaternionMatrixEJA
,cls
).real_unembed(M
)
1981 d
= cls
.dimension_over_reals()
1983 # Use the base ring of the matrix to ensure that its entries can be
1984 # multiplied by elements of the quaternion algebra.
1985 field
= M
.base_ring()
1986 Q
= QuaternionAlgebra(field
,-1,-1)
1989 # Go top-left to bottom-right (reading order), converting every
1990 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
1993 for l
in range(n
/d
):
1994 for m
in range(n
/d
):
1995 submat
= ComplexMatrixEJA
.real_unembed(
1996 M
[d
*l
:d
*l
+d
,d
*m
:d
*m
+d
] )
1997 if submat
[0,0] != submat
[1,1].conjugate():
1998 raise ValueError('bad on-diagonal submatrix')
1999 if submat
[0,1] != -submat
[1,0].conjugate():
2000 raise ValueError('bad off-diagonal submatrix')
2001 z
= submat
[0,0].real()
2002 z
+= submat
[0,0].imag()*i
2003 z
+= submat
[0,1].real()*j
2004 z
+= submat
[0,1].imag()*k
2007 return matrix(Q
, n
/d
, elements
)
2010 class QuaternionHermitianEJA(ConcreteEJA
, QuaternionMatrixEJA
):
2012 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
2013 matrices, the usual symmetric Jordan product, and the
2014 real-part-of-trace inner product. It has dimension `2n^2 - n` over
2019 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2023 In theory, our "field" can be any subfield of the reals::
2025 sage: QuaternionHermitianEJA(2, field=RDF)
2026 Euclidean Jordan algebra of dimension 6 over Real Double Field
2027 sage: QuaternionHermitianEJA(2, field=RR)
2028 Euclidean Jordan algebra of dimension 6 over Real Field with
2029 53 bits of precision
2033 The dimension of this algebra is `2*n^2 - n`::
2035 sage: set_random_seed()
2036 sage: n_max = QuaternionHermitianEJA._max_random_instance_size()
2037 sage: n = ZZ.random_element(1, n_max)
2038 sage: J = QuaternionHermitianEJA(n)
2039 sage: J.dimension() == 2*(n^2) - n
2042 The Jordan multiplication is what we think it is::
2044 sage: set_random_seed()
2045 sage: J = QuaternionHermitianEJA.random_instance()
2046 sage: x,y = J.random_elements(2)
2047 sage: actual = (x*y).to_matrix()
2048 sage: X = x.to_matrix()
2049 sage: Y = y.to_matrix()
2050 sage: expected = (X*Y + Y*X)/2
2051 sage: actual == expected
2053 sage: J(expected) == x*y
2056 We can change the generator prefix::
2058 sage: QuaternionHermitianEJA(2, prefix='a').gens()
2059 (a0, a1, a2, a3, a4, a5)
2061 We can construct the (trivial) algebra of rank zero::
2063 sage: QuaternionHermitianEJA(0)
2064 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2068 def _denormalized_basis(cls
, n
):
2070 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
2072 Why do we embed these? Basically, because all of numerical
2073 linear algebra assumes that you're working with vectors consisting
2074 of `n` entries from a field and scalars from the same field. There's
2075 no way to tell SageMath that (for example) the vectors contain
2076 complex numbers, while the scalar field is real.
2080 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2084 sage: set_random_seed()
2085 sage: n = ZZ.random_element(1,5)
2086 sage: B = QuaternionHermitianEJA._denormalized_basis(n)
2087 sage: all( M.is_symmetric() for M in B )
2092 Q
= QuaternionAlgebra(QQ
,-1,-1)
2095 # This is like the symmetric case, but we need to be careful:
2097 # * We want conjugate-symmetry, not just symmetry.
2098 # * The diagonal will (as a result) be real.
2102 for j
in range(i
+1):
2103 Eij
= matrix(Q
, n
, lambda k
,l
: k
==i
and l
==j
)
2105 Sij
= cls
.real_embed(Eij
)
2108 # The second, third, and fourth ones have a minus
2109 # because they're conjugated.
2110 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
2112 Sij_I
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
2114 Sij_J
= cls
.real_embed(J
*Eij
- J
*Eij
.transpose())
2116 Sij_K
= cls
.real_embed(K
*Eij
- K
*Eij
.transpose())
2119 # Since we embedded these, we can drop back to the "field" that we
2120 # started with instead of the quaternion algebra "Q".
2121 return tuple( s
.change_ring(field
) for s
in S
)
2124 def __init__(self
, n
, **kwargs
):
2125 # We know this is a valid EJA, but will double-check
2126 # if the user passes check_axioms=True.
2127 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2129 super(QuaternionHermitianEJA
, self
).__init
__(self
._denormalized
_basis
(n
),
2130 self
.jordan_product
,
2131 self
.trace_inner_product
,
2133 # TODO: this could be factored out somehow, but is left here
2134 # because the MatrixEJA is not presently a subclass of the
2135 # FDEJA class that defines rank() and one().
2136 self
.rank
.set_cache(n
)
2137 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
2138 self
.one
.set_cache(self(idV
))
2142 def _max_random_instance_size():
2144 The maximum rank of a random QuaternionHermitianEJA.
2146 return 2 # Dimension 6
2149 def random_instance(cls
, **kwargs
):
2151 Return a random instance of this type of algebra.
2153 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2154 return cls(n
, **kwargs
)
2157 class HadamardEJA(ConcreteEJA
):
2159 Return the Euclidean Jordan Algebra corresponding to the set
2160 `R^n` under the Hadamard product.
2162 Note: this is nothing more than the Cartesian product of ``n``
2163 copies of the spin algebra. Once Cartesian product algebras
2164 are implemented, this can go.
2168 sage: from mjo.eja.eja_algebra import HadamardEJA
2172 This multiplication table can be verified by hand::
2174 sage: J = HadamardEJA(3)
2175 sage: e0,e1,e2 = J.gens()
2191 We can change the generator prefix::
2193 sage: HadamardEJA(3, prefix='r').gens()
2197 def __init__(self
, n
, **kwargs
):
2199 jordan_product
= lambda x
,y
: x
2200 inner_product
= lambda x
,y
: x
2202 def jordan_product(x
,y
):
2204 return P( xi
*yi
for (xi
,yi
) in zip(x
,y
) )
2206 def inner_product(x
,y
):
2209 # New defaults for keyword arguments. Don't orthonormalize
2210 # because our basis is already orthonormal with respect to our
2211 # inner-product. Don't check the axioms, because we know this
2212 # is a valid EJA... but do double-check if the user passes
2213 # check_axioms=True. Note: we DON'T override the "check_field"
2214 # default here, because the user can pass in a field!
2215 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2216 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2218 column_basis
= tuple( b
.column() for b
in FreeModule(ZZ
, n
).basis() )
2219 super().__init
__(column_basis
, jordan_product
, inner_product
, **kwargs
)
2220 self
.rank
.set_cache(n
)
2223 self
.one
.set_cache( self
.zero() )
2225 self
.one
.set_cache( sum(self
.gens()) )
2228 def _max_random_instance_size():
2230 The maximum dimension of a random HadamardEJA.
2235 def random_instance(cls
, **kwargs
):
2237 Return a random instance of this type of algebra.
2239 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2240 return cls(n
, **kwargs
)
2243 class BilinearFormEJA(ConcreteEJA
):
2245 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2246 with the half-trace inner product and jordan product ``x*y =
2247 (<Bx,y>,y_bar>, x0*y_bar + y0*x_bar)`` where `B = 1 \times B22` is
2248 a symmetric positive-definite "bilinear form" matrix. Its
2249 dimension is the size of `B`, and it has rank two in dimensions
2250 larger than two. It reduces to the ``JordanSpinEJA`` when `B` is
2251 the identity matrix of order ``n``.
2253 We insist that the one-by-one upper-left identity block of `B` be
2254 passed in as well so that we can be passed a matrix of size zero
2255 to construct a trivial algebra.
2259 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
2260 ....: JordanSpinEJA)
2264 When no bilinear form is specified, the identity matrix is used,
2265 and the resulting algebra is the Jordan spin algebra::
2267 sage: B = matrix.identity(AA,3)
2268 sage: J0 = BilinearFormEJA(B)
2269 sage: J1 = JordanSpinEJA(3)
2270 sage: J0.multiplication_table() == J0.multiplication_table()
2273 An error is raised if the matrix `B` does not correspond to a
2274 positive-definite bilinear form::
2276 sage: B = matrix.random(QQ,2,3)
2277 sage: J = BilinearFormEJA(B)
2278 Traceback (most recent call last):
2280 ValueError: bilinear form is not positive-definite
2281 sage: B = matrix.zero(QQ,3)
2282 sage: J = BilinearFormEJA(B)
2283 Traceback (most recent call last):
2285 ValueError: bilinear form is not positive-definite
2289 We can create a zero-dimensional algebra::
2291 sage: B = matrix.identity(AA,0)
2292 sage: J = BilinearFormEJA(B)
2296 We can check the multiplication condition given in the Jordan, von
2297 Neumann, and Wigner paper (and also discussed on my "On the
2298 symmetry..." paper). Note that this relies heavily on the standard
2299 choice of basis, as does anything utilizing the bilinear form
2300 matrix. We opt not to orthonormalize the basis, because if we
2301 did, we would have to normalize the `s_{i}` in a similar manner::
2303 sage: set_random_seed()
2304 sage: n = ZZ.random_element(5)
2305 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2306 sage: B11 = matrix.identity(QQ,1)
2307 sage: B22 = M.transpose()*M
2308 sage: B = block_matrix(2,2,[ [B11,0 ],
2310 sage: J = BilinearFormEJA(B, orthonormalize=False)
2311 sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
2312 sage: V = J.vector_space()
2313 sage: sis = [ J( V([0] + (M.inverse()*ei).list()).column() )
2314 ....: for ei in eis ]
2315 sage: actual = [ sis[i]*sis[j]
2316 ....: for i in range(n-1)
2317 ....: for j in range(n-1) ]
2318 sage: expected = [ J.one() if i == j else J.zero()
2319 ....: for i in range(n-1)
2320 ....: for j in range(n-1) ]
2321 sage: actual == expected
2325 def __init__(self
, B
, **kwargs
):
2326 # The matrix "B" is supplied by the user in most cases,
2327 # so it makes sense to check whether or not its positive-
2328 # definite unless we are specifically asked not to...
2329 if ("check_axioms" not in kwargs
) or kwargs
["check_axioms"]:
2330 if not B
.is_positive_definite():
2331 raise ValueError("bilinear form is not positive-definite")
2333 # However, all of the other data for this EJA is computed
2334 # by us in manner that guarantees the axioms are
2335 # satisfied. So, again, unless we are specifically asked to
2336 # verify things, we'll skip the rest of the checks.
2337 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2339 def inner_product(x
,y
):
2340 return (y
.T
*B
*x
)[0,0]
2342 def jordan_product(x
,y
):
2348 z0
= inner_product(y
,x
)
2349 zbar
= y0
*xbar
+ x0
*ybar
2350 return P([z0
] + zbar
.list())
2353 column_basis
= tuple( b
.column() for b
in FreeModule(ZZ
, n
).basis() )
2354 super(BilinearFormEJA
, self
).__init
__(column_basis
,
2359 # The rank of this algebra is two, unless we're in a
2360 # one-dimensional ambient space (because the rank is bounded
2361 # by the ambient dimension).
2362 self
.rank
.set_cache(min(n
,2))
2365 self
.one
.set_cache( self
.zero() )
2367 self
.one
.set_cache( self
.monomial(0) )
2370 def _max_random_instance_size():
2372 The maximum dimension of a random BilinearFormEJA.
2377 def random_instance(cls
, **kwargs
):
2379 Return a random instance of this algebra.
2381 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2383 B
= matrix
.identity(ZZ
, n
)
2384 return cls(B
, **kwargs
)
2386 B11
= matrix
.identity(ZZ
, 1)
2387 M
= matrix
.random(ZZ
, n
-1)
2388 I
= matrix
.identity(ZZ
, n
-1)
2390 while alpha
.is_zero():
2391 alpha
= ZZ
.random_element().abs()
2392 B22
= M
.transpose()*M
+ alpha
*I
2394 from sage
.matrix
.special
import block_matrix
2395 B
= block_matrix(2,2, [ [B11
, ZZ(0) ],
2398 return cls(B
, **kwargs
)
2401 class JordanSpinEJA(BilinearFormEJA
):
2403 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2404 with the usual inner product and jordan product ``x*y =
2405 (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2410 sage: from mjo.eja.eja_algebra import JordanSpinEJA
2414 This multiplication table can be verified by hand::
2416 sage: J = JordanSpinEJA(4)
2417 sage: e0,e1,e2,e3 = J.gens()
2433 We can change the generator prefix::
2435 sage: JordanSpinEJA(2, prefix='B').gens()
2440 Ensure that we have the usual inner product on `R^n`::
2442 sage: set_random_seed()
2443 sage: J = JordanSpinEJA.random_instance()
2444 sage: x,y = J.random_elements(2)
2445 sage: actual = x.inner_product(y)
2446 sage: expected = x.to_vector().inner_product(y.to_vector())
2447 sage: actual == expected
2451 def __init__(self
, n
, **kwargs
):
2452 # This is a special case of the BilinearFormEJA with the
2453 # identity matrix as its bilinear form.
2454 B
= matrix
.identity(ZZ
, n
)
2456 # Don't orthonormalize because our basis is already
2457 # orthonormal with respect to our inner-product.
2458 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2460 # But also don't pass check_field=False here, because the user
2461 # can pass in a field!
2462 super(JordanSpinEJA
, self
).__init
__(B
, **kwargs
)
2465 def _max_random_instance_size():
2467 The maximum dimension of a random JordanSpinEJA.
2472 def random_instance(cls
, **kwargs
):
2474 Return a random instance of this type of algebra.
2476 Needed here to override the implementation for ``BilinearFormEJA``.
2478 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2479 return cls(n
, **kwargs
)
2482 class TrivialEJA(ConcreteEJA
):
2484 The trivial Euclidean Jordan algebra consisting of only a zero element.
2488 sage: from mjo.eja.eja_algebra import TrivialEJA
2492 sage: J = TrivialEJA()
2499 sage: 7*J.one()*12*J.one()
2501 sage: J.one().inner_product(J.one())
2503 sage: J.one().norm()
2505 sage: J.one().subalgebra_generated_by()
2506 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2511 def __init__(self
, **kwargs
):
2512 jordan_product
= lambda x
,y
: x
2513 inner_product
= lambda x
,y
: 0
2516 # New defaults for keyword arguments
2517 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2518 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2520 super(TrivialEJA
, self
).__init
__(basis
,
2524 # The rank is zero using my definition, namely the dimension of the
2525 # largest subalgebra generated by any element.
2526 self
.rank
.set_cache(0)
2527 self
.one
.set_cache( self
.zero() )
2530 def random_instance(cls
, **kwargs
):
2531 # We don't take a "size" argument so the superclass method is
2532 # inappropriate for us.
2533 return cls(**kwargs
)
2535 # class DirectSumEJA(ConcreteEJA):
2537 # The external (orthogonal) direct sum of two other Euclidean Jordan
2538 # algebras. Essentially the Cartesian product of its two factors.
2539 # Every Euclidean Jordan algebra decomposes into an orthogonal
2540 # direct sum of simple Euclidean Jordan algebras, so no generality
2541 # is lost by providing only this construction.
2545 # sage: from mjo.eja.eja_algebra import (random_eja,
2546 # ....: HadamardEJA,
2547 # ....: RealSymmetricEJA,
2548 # ....: DirectSumEJA)
2552 # sage: J1 = HadamardEJA(2)
2553 # sage: J2 = RealSymmetricEJA(3)
2554 # sage: J = DirectSumEJA(J1,J2)
2555 # sage: J.dimension()
2562 # The external direct sum construction is only valid when the two factors
2563 # have the same base ring; an error is raised otherwise::
2565 # sage: set_random_seed()
2566 # sage: J1 = random_eja(field=AA)
2567 # sage: J2 = random_eja(field=QQ,orthonormalize=False)
2568 # sage: J = DirectSumEJA(J1,J2)
2569 # Traceback (most recent call last):
2571 # ValueError: algebras must share the same base field
2574 # def __init__(self, J1, J2, **kwargs):
2575 # if J1.base_ring() != J2.base_ring():
2576 # raise ValueError("algebras must share the same base field")
2577 # field = J1.base_ring()
2579 # self._factors = (J1, J2)
2580 # n1 = J1.dimension()
2581 # n2 = J2.dimension()
2583 # V = VectorSpace(field, n)
2584 # mult_table = [ [ V.zero() for j in range(i+1) ]
2585 # for i in range(n) ]
2586 # for i in range(n1):
2587 # for j in range(i+1):
2588 # p = (J1.monomial(i)*J1.monomial(j)).to_vector()
2589 # mult_table[i][j] = V(p.list() + [field.zero()]*n2)
2591 # for i in range(n2):
2592 # for j in range(i+1):
2593 # p = (J2.monomial(i)*J2.monomial(j)).to_vector()
2594 # mult_table[n1+i][n1+j] = V([field.zero()]*n1 + p.list())
2596 # # TODO: build the IP table here from the two constituent IP
2597 # # matrices (it'll be block diagonal, I think).
2598 # ip_table = [ [ field.zero() for j in range(i+1) ]
2599 # for i in range(n) ]
2600 # super(DirectSumEJA, self).__init__(field,
2603 # check_axioms=False,
2605 # self.rank.set_cache(J1.rank() + J2.rank())
2608 # def factors(self):
2610 # Return the pair of this algebra's factors.
2614 # sage: from mjo.eja.eja_algebra import (HadamardEJA,
2615 # ....: JordanSpinEJA,
2616 # ....: DirectSumEJA)
2620 # sage: J1 = HadamardEJA(2, field=QQ)
2621 # sage: J2 = JordanSpinEJA(3, field=QQ)
2622 # sage: J = DirectSumEJA(J1,J2)
2624 # (Euclidean Jordan algebra of dimension 2 over Rational Field,
2625 # Euclidean Jordan algebra of dimension 3 over Rational Field)
2628 # return self._factors
2630 # def projections(self):
2632 # Return a pair of projections onto this algebra's factors.
2636 # sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
2637 # ....: ComplexHermitianEJA,
2638 # ....: DirectSumEJA)
2642 # sage: J1 = JordanSpinEJA(2)
2643 # sage: J2 = ComplexHermitianEJA(2)
2644 # sage: J = DirectSumEJA(J1,J2)
2645 # sage: (pi_left, pi_right) = J.projections()
2646 # sage: J.one().to_vector()
2647 # (1, 0, 1, 0, 0, 1)
2648 # sage: pi_left(J.one()).to_vector()
2650 # sage: pi_right(J.one()).to_vector()
2654 # (J1,J2) = self.factors()
2655 # m = J1.dimension()
2656 # n = J2.dimension()
2657 # V_basis = self.vector_space().basis()
2658 # # Need to specify the dimensions explicitly so that we don't
2659 # # wind up with a zero-by-zero matrix when we want e.g. a
2660 # # zero-by-two matrix (important for composing things).
2661 # P1 = matrix(self.base_ring(), m, m+n, V_basis[:m])
2662 # P2 = matrix(self.base_ring(), n, m+n, V_basis[m:])
2663 # pi_left = FiniteDimensionalEJAOperator(self,J1,P1)
2664 # pi_right = FiniteDimensionalEJAOperator(self,J2,P2)
2665 # return (pi_left, pi_right)
2667 # def inclusions(self):
2669 # Return the pair of inclusion maps from our factors into us.
2673 # sage: from mjo.eja.eja_algebra import (random_eja,
2674 # ....: JordanSpinEJA,
2675 # ....: RealSymmetricEJA,
2676 # ....: DirectSumEJA)
2680 # sage: J1 = JordanSpinEJA(3)
2681 # sage: J2 = RealSymmetricEJA(2)
2682 # sage: J = DirectSumEJA(J1,J2)
2683 # sage: (iota_left, iota_right) = J.inclusions()
2684 # sage: iota_left(J1.zero()) == J.zero()
2686 # sage: iota_right(J2.zero()) == J.zero()
2688 # sage: J1.one().to_vector()
2690 # sage: iota_left(J1.one()).to_vector()
2691 # (1, 0, 0, 0, 0, 0)
2692 # sage: J2.one().to_vector()
2694 # sage: iota_right(J2.one()).to_vector()
2695 # (0, 0, 0, 1, 0, 1)
2696 # sage: J.one().to_vector()
2697 # (1, 0, 0, 1, 0, 1)
2701 # Composing a projection with the corresponding inclusion should
2702 # produce the identity map, and mismatching them should produce
2705 # sage: set_random_seed()
2706 # sage: J1 = random_eja()
2707 # sage: J2 = random_eja()
2708 # sage: J = DirectSumEJA(J1,J2)
2709 # sage: (iota_left, iota_right) = J.inclusions()
2710 # sage: (pi_left, pi_right) = J.projections()
2711 # sage: pi_left*iota_left == J1.one().operator()
2713 # sage: pi_right*iota_right == J2.one().operator()
2715 # sage: (pi_left*iota_right).is_zero()
2717 # sage: (pi_right*iota_left).is_zero()
2721 # (J1,J2) = self.factors()
2722 # m = J1.dimension()
2723 # n = J2.dimension()
2724 # V_basis = self.vector_space().basis()
2725 # # Need to specify the dimensions explicitly so that we don't
2726 # # wind up with a zero-by-zero matrix when we want e.g. a
2727 # # two-by-zero matrix (important for composing things).
2728 # I1 = matrix.column(self.base_ring(), m, m+n, V_basis[:m])
2729 # I2 = matrix.column(self.base_ring(), n, m+n, V_basis[m:])
2730 # iota_left = FiniteDimensionalEJAOperator(J1,self,I1)
2731 # iota_right = FiniteDimensionalEJAOperator(J2,self,I2)
2732 # return (iota_left, iota_right)
2734 # def inner_product(self, x, y):
2736 # The standard Cartesian inner-product.
2738 # We project ``x`` and ``y`` onto our factors, and add up the
2739 # inner-products from the subalgebras.
2744 # sage: from mjo.eja.eja_algebra import (HadamardEJA,
2745 # ....: QuaternionHermitianEJA,
2746 # ....: DirectSumEJA)
2750 # sage: J1 = HadamardEJA(3,field=QQ)
2751 # sage: J2 = QuaternionHermitianEJA(2,field=QQ,orthonormalize=False)
2752 # sage: J = DirectSumEJA(J1,J2)
2753 # sage: x1 = J1.one()
2755 # sage: y1 = J2.one()
2757 # sage: x1.inner_product(x2)
2759 # sage: y1.inner_product(y2)
2761 # sage: J.one().inner_product(J.one())
2765 # (pi_left, pi_right) = self.projections()
2771 # return (x1.inner_product(y1) + x2.inner_product(y2))
2775 random_eja
= ConcreteEJA
.random_instance