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
= 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
.zero(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 elt
= jordan_product(q_i
, q_j
)
175 ip
= inner_product(q_i
, q_j
)
177 # The jordan product returns a matrixy answer, so we
178 # have to convert it to the algebra coordinates.
179 elt
= W
.coordinate_vector(V(elt
.list()))
180 self
._multiplication
_table
[i
][j
] = self
.from_vector(elt
)
181 self
._inner
_product
_matrix
[i
,j
] = ip
182 self
._inner
_product
_matrix
[j
,i
] = ip
184 self
._inner
_product
_matrix
._cache
= {'hermitian': True}
185 self
._inner
_product
_matrix
.set_immutable()
188 if not self
._is
_jordanian
():
189 raise ValueError("Jordan identity does not hold")
190 if not self
._inner
_product
_is
_associative
():
191 raise ValueError("inner product is not associative")
194 def _coerce_map_from_base_ring(self
):
196 Disable the map from the base ring into the algebra.
198 Performing a nonsense conversion like this automatically
199 is counterpedagogical. The fallback is to try the usual
200 element constructor, which should also fail.
204 sage: from mjo.eja.eja_algebra import random_eja
208 sage: set_random_seed()
209 sage: J = random_eja()
211 Traceback (most recent call last):
213 ValueError: not an element of this algebra
219 def product_on_basis(self
, i
, j
):
220 # We only stored the lower-triangular portion of the
221 # multiplication table.
223 return self
._multiplication
_table
[i
][j
]
225 return self
._multiplication
_table
[j
][i
]
227 def inner_product(self
, x
, y
):
229 The inner product associated with this Euclidean Jordan algebra.
231 Defaults to the trace inner product, but can be overridden by
232 subclasses if they are sure that the necessary properties are
237 sage: from mjo.eja.eja_algebra import (random_eja,
239 ....: BilinearFormEJA)
243 Our inner product is "associative," which means the following for
244 a symmetric bilinear form::
246 sage: set_random_seed()
247 sage: J = random_eja()
248 sage: x,y,z = J.random_elements(3)
249 sage: (x*y).inner_product(z) == y.inner_product(x*z)
254 Ensure that this is the usual inner product for the algebras
257 sage: set_random_seed()
258 sage: J = HadamardEJA.random_instance()
259 sage: x,y = J.random_elements(2)
260 sage: actual = x.inner_product(y)
261 sage: expected = x.to_vector().inner_product(y.to_vector())
262 sage: actual == expected
265 Ensure that this is one-half of the trace inner-product in a
266 BilinearFormEJA that isn't just the reals (when ``n`` isn't
267 one). This is in Faraut and Koranyi, and also my "On the
270 sage: set_random_seed()
271 sage: J = BilinearFormEJA.random_instance()
272 sage: n = J.dimension()
273 sage: x = J.random_element()
274 sage: y = J.random_element()
275 sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
278 B
= self
._inner
_product
_matrix
279 return (B
*x
.to_vector()).inner_product(y
.to_vector())
282 def _is_commutative(self
):
284 Whether or not this algebra's multiplication table is commutative.
286 This method should of course always return ``True``, unless
287 this algebra was constructed with ``check_axioms=False`` and
288 passed an invalid multiplication table.
290 return all( self
.product_on_basis(i
,j
) == self
.product_on_basis(i
,j
)
291 for i
in range(self
.dimension())
292 for j
in range(self
.dimension()) )
294 def _is_jordanian(self
):
296 Whether or not this algebra's multiplication table respects the
297 Jordan identity `(x^{2})(xy) = x(x^{2}y)`.
299 We only check one arrangement of `x` and `y`, so for a
300 ``True`` result to be truly true, you should also check
301 :meth:`_is_commutative`. This method should of course always
302 return ``True``, unless this algebra was constructed with
303 ``check_axioms=False`` and passed an invalid multiplication table.
305 return all( (self
.monomial(i
)**2)*(self
.monomial(i
)*self
.monomial(j
))
307 (self
.monomial(i
))*((self
.monomial(i
)**2)*self
.monomial(j
))
308 for i
in range(self
.dimension())
309 for j
in range(self
.dimension()) )
311 def _inner_product_is_associative(self
):
313 Return whether or not this algebra's inner product `B` is
314 associative; that is, whether or not `B(xy,z) = B(x,yz)`.
316 This method should of course always return ``True``, unless
317 this algebra was constructed with ``check_axioms=False`` and
318 passed an invalid multiplication table.
321 # Used to check whether or not something is zero in an inexact
322 # ring. This number is sufficient to allow the construction of
323 # QuaternionHermitianEJA(2, field=RDF) with check_axioms=True.
326 for i
in range(self
.dimension()):
327 for j
in range(self
.dimension()):
328 for k
in range(self
.dimension()):
332 diff
= (x
*y
).inner_product(z
) - x
.inner_product(y
*z
)
334 if self
.base_ring().is_exact():
338 if diff
.abs() > epsilon
:
343 def _element_constructor_(self
, elt
):
345 Construct an element of this algebra from its vector or matrix
348 This gets called only after the parent element _call_ method
349 fails to find a coercion for the argument.
353 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
355 ....: RealSymmetricEJA)
359 The identity in `S^n` is converted to the identity in the EJA::
361 sage: J = RealSymmetricEJA(3)
362 sage: I = matrix.identity(QQ,3)
363 sage: J(I) == J.one()
366 This skew-symmetric matrix can't be represented in the EJA::
368 sage: J = RealSymmetricEJA(3)
369 sage: A = matrix(QQ,3, lambda i,j: i-j)
371 Traceback (most recent call last):
373 ValueError: not an element of this algebra
377 Ensure that we can convert any element of the two non-matrix
378 simple algebras (whose matrix representations are columns)
379 back and forth faithfully::
381 sage: set_random_seed()
382 sage: J = HadamardEJA.random_instance()
383 sage: x = J.random_element()
384 sage: J(x.to_vector().column()) == x
386 sage: J = JordanSpinEJA.random_instance()
387 sage: x = J.random_element()
388 sage: J(x.to_vector().column()) == x
392 msg
= "not an element of this algebra"
394 # The superclass implementation of random_element()
395 # needs to be able to coerce "0" into the algebra.
397 elif elt
in self
.base_ring():
398 # Ensure that no base ring -> algebra coercion is performed
399 # by this method. There's some stupidity in sage that would
400 # otherwise propagate to this method; for example, sage thinks
401 # that the integer 3 belongs to the space of 2-by-2 matrices.
402 raise ValueError(msg
)
406 except (AttributeError, TypeError):
407 # Try to convert a vector into a column-matrix
410 if elt
not in self
.matrix_space():
411 raise ValueError(msg
)
413 # Thanks for nothing! Matrix spaces aren't vector spaces in
414 # Sage, so we have to figure out its matrix-basis coordinates
415 # ourselves. We use the basis space's ring instead of the
416 # element's ring because the basis space might be an algebraic
417 # closure whereas the base ring of the 3-by-3 identity matrix
418 # could be QQ instead of QQbar.
420 # We pass check=False because the matrix basis is "guaranteed"
421 # to be linearly independent... right? Ha ha.
422 V
= VectorSpace(self
.base_ring(), elt
.nrows()*elt
.ncols())
423 W
= V
.span_of_basis( (_mat2vec(s
) for s
in self
.matrix_basis()),
427 coords
= W
.coordinate_vector(_mat2vec(elt
))
428 except ArithmeticError: # vector is not in free module
429 raise ValueError(msg
)
431 return self
.from_vector(coords
)
435 Return a string representation of ``self``.
439 sage: from mjo.eja.eja_algebra import JordanSpinEJA
443 Ensure that it says what we think it says::
445 sage: JordanSpinEJA(2, field=AA)
446 Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
447 sage: JordanSpinEJA(3, field=RDF)
448 Euclidean Jordan algebra of dimension 3 over Real Double Field
451 fmt
= "Euclidean Jordan algebra of dimension {} over {}"
452 return fmt
.format(self
.dimension(), self
.base_ring())
456 def characteristic_polynomial_of(self
):
458 Return the algebra's "characteristic polynomial of" function,
459 which is itself a multivariate polynomial that, when evaluated
460 at the coordinates of some algebra element, returns that
461 element's characteristic polynomial.
463 The resulting polynomial has `n+1` variables, where `n` is the
464 dimension of this algebra. The first `n` variables correspond to
465 the coordinates of an algebra element: when evaluated at the
466 coordinates of an algebra element with respect to a certain
467 basis, the result is a univariate polynomial (in the one
468 remaining variable ``t``), namely the characteristic polynomial
473 sage: from mjo.eja.eja_algebra import JordanSpinEJA, TrivialEJA
477 The characteristic polynomial in the spin algebra is given in
478 Alizadeh, Example 11.11::
480 sage: J = JordanSpinEJA(3)
481 sage: p = J.characteristic_polynomial_of(); p
482 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
483 sage: xvec = J.one().to_vector()
487 By definition, the characteristic polynomial is a monic
488 degree-zero polynomial in a rank-zero algebra. Note that
489 Cayley-Hamilton is indeed satisfied since the polynomial
490 ``1`` evaluates to the identity element of the algebra on
493 sage: J = TrivialEJA()
494 sage: J.characteristic_polynomial_of()
501 # The list of coefficient polynomials a_0, a_1, a_2, ..., a_(r-1).
502 a
= self
._charpoly
_coefficients
()
504 # We go to a bit of trouble here to reorder the
505 # indeterminates, so that it's easier to evaluate the
506 # characteristic polynomial at x's coordinates and get back
507 # something in terms of t, which is what we want.
508 S
= PolynomialRing(self
.base_ring(),'t')
512 S
= PolynomialRing(S
, R
.variable_names())
515 return (t
**r
+ sum( a
[k
]*(t
**k
) for k
in range(r
) ))
517 def coordinate_polynomial_ring(self
):
519 The multivariate polynomial ring in which this algebra's
520 :meth:`characteristic_polynomial_of` lives.
524 sage: from mjo.eja.eja_algebra import (HadamardEJA,
525 ....: RealSymmetricEJA)
529 sage: J = HadamardEJA(2)
530 sage: J.coordinate_polynomial_ring()
531 Multivariate Polynomial Ring in X1, X2...
532 sage: J = RealSymmetricEJA(3,field=QQ,orthonormalize=False)
533 sage: J.coordinate_polynomial_ring()
534 Multivariate Polynomial Ring in X1, X2, X3, X4, X5, X6...
537 var_names
= tuple( "X%d" % z
for z
in range(1, self
.dimension()+1) )
538 return PolynomialRing(self
.base_ring(), var_names
)
540 def inner_product(self
, x
, y
):
542 The inner product associated with this Euclidean Jordan algebra.
544 Defaults to the trace inner product, but can be overridden by
545 subclasses if they are sure that the necessary properties are
550 sage: from mjo.eja.eja_algebra import (random_eja,
552 ....: BilinearFormEJA)
556 Our inner product is "associative," which means the following for
557 a symmetric bilinear form::
559 sage: set_random_seed()
560 sage: J = random_eja()
561 sage: x,y,z = J.random_elements(3)
562 sage: (x*y).inner_product(z) == y.inner_product(x*z)
567 Ensure that this is the usual inner product for the algebras
570 sage: set_random_seed()
571 sage: J = HadamardEJA.random_instance()
572 sage: x,y = J.random_elements(2)
573 sage: actual = x.inner_product(y)
574 sage: expected = x.to_vector().inner_product(y.to_vector())
575 sage: actual == expected
578 Ensure that this is one-half of the trace inner-product in a
579 BilinearFormEJA that isn't just the reals (when ``n`` isn't
580 one). This is in Faraut and Koranyi, and also my "On the
583 sage: set_random_seed()
584 sage: J = BilinearFormEJA.random_instance()
585 sage: n = J.dimension()
586 sage: x = J.random_element()
587 sage: y = J.random_element()
588 sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
591 B
= self
._inner
_product
_matrix
592 return (B
*x
.to_vector()).inner_product(y
.to_vector())
595 def is_trivial(self
):
597 Return whether or not this algebra is trivial.
599 A trivial algebra contains only the zero element.
603 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
608 sage: J = ComplexHermitianEJA(3)
614 sage: J = TrivialEJA()
619 return self
.dimension() == 0
622 def multiplication_table(self
):
624 Return a visual representation of this algebra's multiplication
625 table (on basis elements).
629 sage: from mjo.eja.eja_algebra import JordanSpinEJA
633 sage: J = JordanSpinEJA(4)
634 sage: J.multiplication_table()
635 +----++----+----+----+----+
636 | * || e0 | e1 | e2 | e3 |
637 +====++====+====+====+====+
638 | e0 || e0 | e1 | e2 | e3 |
639 +----++----+----+----+----+
640 | e1 || e1 | e0 | 0 | 0 |
641 +----++----+----+----+----+
642 | e2 || e2 | 0 | e0 | 0 |
643 +----++----+----+----+----+
644 | e3 || e3 | 0 | 0 | e0 |
645 +----++----+----+----+----+
649 # Prepend the header row.
650 M
= [["*"] + list(self
.gens())]
652 # And to each subsequent row, prepend an entry that belongs to
653 # the left-side "header column."
654 M
+= [ [self
.monomial(i
)] + [ self
.product_on_basis(i
,j
)
658 return table(M
, header_row
=True, header_column
=True, frame
=True)
661 def matrix_basis(self
):
663 Return an (often more natural) representation of this algebras
664 basis as an ordered tuple of matrices.
666 Every finite-dimensional Euclidean Jordan Algebra is a, up to
667 Jordan isomorphism, a direct sum of five simple
668 algebras---four of which comprise Hermitian matrices. And the
669 last type of algebra can of course be thought of as `n`-by-`1`
670 column matrices (ambiguusly called column vectors) to avoid
671 special cases. As a result, matrices (and column vectors) are
672 a natural representation format for Euclidean Jordan algebra
675 But, when we construct an algebra from a basis of matrices,
676 those matrix representations are lost in favor of coordinate
677 vectors *with respect to* that basis. We could eventually
678 convert back if we tried hard enough, but having the original
679 representations handy is valuable enough that we simply store
680 them and return them from this method.
682 Why implement this for non-matrix algebras? Avoiding special
683 cases for the :class:`BilinearFormEJA` pays with simplicity in
684 its own right. But mainly, we would like to be able to assume
685 that elements of a :class:`DirectSumEJA` can be displayed
686 nicely, without having to have special classes for direct sums
687 one of whose components was a matrix algebra.
691 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
692 ....: RealSymmetricEJA)
696 sage: J = RealSymmetricEJA(2)
698 Finite family {0: e0, 1: e1, 2: e2}
699 sage: J.matrix_basis()
701 [1 0] [ 0 0.7071067811865475?] [0 0]
702 [0 0], [0.7071067811865475? 0], [0 1]
707 sage: J = JordanSpinEJA(2)
709 Finite family {0: e0, 1: e1}
710 sage: J.matrix_basis()
716 return self
._matrix
_basis
719 def matrix_space(self
):
721 Return the matrix space in which this algebra's elements live, if
722 we think of them as matrices (including column vectors of the
725 Generally this will be an `n`-by-`1` column-vector space,
726 except when the algebra is trivial. There it's `n`-by-`n`
727 (where `n` is zero), to ensure that two elements of the matrix
728 space (empty matrices) can be multiplied.
730 Matrix algebras override this with something more useful.
732 if self
.is_trivial():
733 return MatrixSpace(self
.base_ring(), 0)
735 return self
._matrix
_basis
[0].matrix_space()
741 Return the unit element of this algebra.
745 sage: from mjo.eja.eja_algebra import (HadamardEJA,
750 sage: J = HadamardEJA(5)
752 e0 + e1 + e2 + e3 + e4
756 The identity element acts like the identity::
758 sage: set_random_seed()
759 sage: J = random_eja()
760 sage: x = J.random_element()
761 sage: J.one()*x == x and x*J.one() == x
764 The matrix of the unit element's operator is the identity::
766 sage: set_random_seed()
767 sage: J = random_eja()
768 sage: actual = J.one().operator().matrix()
769 sage: expected = matrix.identity(J.base_ring(), J.dimension())
770 sage: actual == expected
773 Ensure that the cached unit element (often precomputed by
774 hand) agrees with the computed one::
776 sage: set_random_seed()
777 sage: J = random_eja()
778 sage: cached = J.one()
779 sage: J.one.clear_cache()
780 sage: J.one() == cached
784 # We can brute-force compute the matrices of the operators
785 # that correspond to the basis elements of this algebra.
786 # If some linear combination of those basis elements is the
787 # algebra identity, then the same linear combination of
788 # their matrices has to be the identity matrix.
790 # Of course, matrices aren't vectors in sage, so we have to
791 # appeal to the "long vectors" isometry.
792 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
794 # Now we use basic linear algebra to find the coefficients,
795 # of the matrices-as-vectors-linear-combination, which should
796 # work for the original algebra basis too.
797 A
= matrix(self
.base_ring(), oper_vecs
)
799 # We used the isometry on the left-hand side already, but we
800 # still need to do it for the right-hand side. Recall that we
801 # wanted something that summed to the identity matrix.
802 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
804 # Now if there's an identity element in the algebra, this
805 # should work. We solve on the left to avoid having to
806 # transpose the matrix "A".
807 return self
.from_vector(A
.solve_left(b
))
810 def peirce_decomposition(self
, c
):
812 The Peirce decomposition of this algebra relative to the
815 In the future, this can be extended to a complete system of
816 orthogonal idempotents.
820 - ``c`` -- an idempotent of this algebra.
824 A triple (J0, J5, J1) containing two subalgebras and one subspace
827 - ``J0`` -- the algebra on the eigenspace of ``c.operator()``
828 corresponding to the eigenvalue zero.
830 - ``J5`` -- the eigenspace (NOT a subalgebra) of ``c.operator()``
831 corresponding to the eigenvalue one-half.
833 - ``J1`` -- the algebra on the eigenspace of ``c.operator()``
834 corresponding to the eigenvalue one.
836 These are the only possible eigenspaces for that operator, and this
837 algebra is a direct sum of them. The spaces ``J0`` and ``J1`` are
838 orthogonal, and are subalgebras of this algebra with the appropriate
843 sage: from mjo.eja.eja_algebra import random_eja, RealSymmetricEJA
847 The canonical example comes from the symmetric matrices, which
848 decompose into diagonal and off-diagonal parts::
850 sage: J = RealSymmetricEJA(3)
851 sage: C = matrix(QQ, [ [1,0,0],
855 sage: J0,J5,J1 = J.peirce_decomposition(c)
857 Euclidean Jordan algebra of dimension 1...
859 Vector space of degree 6 and dimension 2...
861 Euclidean Jordan algebra of dimension 3...
862 sage: J0.one().to_matrix()
866 sage: orig_df = AA.options.display_format
867 sage: AA.options.display_format = 'radical'
868 sage: J.from_vector(J5.basis()[0]).to_matrix()
872 sage: J.from_vector(J5.basis()[1]).to_matrix()
876 sage: AA.options.display_format = orig_df
877 sage: J1.one().to_matrix()
884 Every algebra decomposes trivially with respect to its identity
887 sage: set_random_seed()
888 sage: J = random_eja()
889 sage: J0,J5,J1 = J.peirce_decomposition(J.one())
890 sage: J0.dimension() == 0 and J5.dimension() == 0
892 sage: J1.superalgebra() == J and J1.dimension() == J.dimension()
895 The decomposition is into eigenspaces, and its components are
896 therefore necessarily orthogonal. Moreover, the identity
897 elements in the two subalgebras are the projections onto their
898 respective subspaces of the superalgebra's identity element::
900 sage: set_random_seed()
901 sage: J = random_eja()
902 sage: x = J.random_element()
903 sage: if not J.is_trivial():
904 ....: while x.is_nilpotent():
905 ....: x = J.random_element()
906 sage: c = x.subalgebra_idempotent()
907 sage: J0,J5,J1 = J.peirce_decomposition(c)
909 sage: for (w,y,z) in zip(J0.basis(), J5.basis(), J1.basis()):
910 ....: w = w.superalgebra_element()
911 ....: y = J.from_vector(y)
912 ....: z = z.superalgebra_element()
913 ....: ipsum += w.inner_product(y).abs()
914 ....: ipsum += w.inner_product(z).abs()
915 ....: ipsum += y.inner_product(z).abs()
918 sage: J1(c) == J1.one()
920 sage: J0(J.one() - c) == J0.one()
924 if not c
.is_idempotent():
925 raise ValueError("element is not idempotent: %s" % c
)
927 from mjo
.eja
.eja_subalgebra
import FiniteDimensionalEJASubalgebra
929 # Default these to what they should be if they turn out to be
930 # trivial, because eigenspaces_left() won't return eigenvalues
931 # corresponding to trivial spaces (e.g. it returns only the
932 # eigenspace corresponding to lambda=1 if you take the
933 # decomposition relative to the identity element).
934 trivial
= FiniteDimensionalEJASubalgebra(self
, ())
935 J0
= trivial
# eigenvalue zero
936 J5
= VectorSpace(self
.base_ring(), 0) # eigenvalue one-half
937 J1
= trivial
# eigenvalue one
939 for (eigval
, eigspace
) in c
.operator().matrix().right_eigenspaces():
940 if eigval
== ~
(self
.base_ring()(2)):
943 gens
= tuple( self
.from_vector(b
) for b
in eigspace
.basis() )
944 subalg
= FiniteDimensionalEJASubalgebra(self
,
952 raise ValueError("unexpected eigenvalue: %s" % eigval
)
957 def random_element(self
, thorough
=False):
959 Return a random element of this algebra.
961 Our algebra superclass method only returns a linear
962 combination of at most two basis elements. We instead
963 want the vector space "random element" method that
964 returns a more diverse selection.
968 - ``thorough`` -- (boolean; default False) whether or not we
969 should generate irrational coefficients for the random
970 element when our base ring is irrational; this slows the
971 algebra operations to a crawl, but any truly random method
975 # For a general base ring... maybe we can trust this to do the
976 # right thing? Unlikely, but.
977 V
= self
.vector_space()
978 v
= V
.random_element()
980 if self
.base_ring() is AA
:
981 # The "random element" method of the algebraic reals is
982 # stupid at the moment, and only returns integers between
983 # -2 and 2, inclusive:
985 # https://trac.sagemath.org/ticket/30875
987 # Instead, we implement our own "random vector" method,
988 # and then coerce that into the algebra. We use the vector
989 # space degree here instead of the dimension because a
990 # subalgebra could (for example) be spanned by only two
991 # vectors, each with five coordinates. We need to
992 # generate all five coordinates.
994 v
*= QQbar
.random_element().real()
996 v
*= QQ
.random_element()
998 return self
.from_vector(V
.coordinate_vector(v
))
1000 def random_elements(self
, count
, thorough
=False):
1002 Return ``count`` random elements as a tuple.
1006 - ``thorough`` -- (boolean; default False) whether or not we
1007 should generate irrational coefficients for the random
1008 elements when our base ring is irrational; this slows the
1009 algebra operations to a crawl, but any truly random method
1014 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1018 sage: J = JordanSpinEJA(3)
1019 sage: x,y,z = J.random_elements(3)
1020 sage: all( [ x in J, y in J, z in J ])
1022 sage: len( J.random_elements(10) ) == 10
1026 return tuple( self
.random_element(thorough
)
1027 for idx
in range(count
) )
1031 def _charpoly_coefficients(self
):
1033 The `r` polynomial coefficients of the "characteristic polynomial
1036 n
= self
.dimension()
1037 R
= self
.coordinate_polynomial_ring()
1039 F
= R
.fraction_field()
1042 # From a result in my book, these are the entries of the
1043 # basis representation of L_x.
1044 return sum( vars[k
]*self
.monomial(k
).operator().matrix()[i
,j
]
1047 L_x
= matrix(F
, n
, n
, L_x_i_j
)
1050 if self
.rank
.is_in_cache():
1052 # There's no need to pad the system with redundant
1053 # columns if we *know* they'll be redundant.
1056 # Compute an extra power in case the rank is equal to
1057 # the dimension (otherwise, we would stop at x^(r-1)).
1058 x_powers
= [ (L_x
**k
)*self
.one().to_vector()
1059 for k
in range(n
+1) ]
1060 A
= matrix
.column(F
, x_powers
[:n
])
1061 AE
= A
.extended_echelon_form()
1068 # The theory says that only the first "r" coefficients are
1069 # nonzero, and they actually live in the original polynomial
1070 # ring and not the fraction field. We negate them because
1071 # in the actual characteristic polynomial, they get moved
1072 # to the other side where x^r lives.
1073 return -A_rref
.solve_right(E
*b
).change_ring(R
)[:r
]
1078 Return the rank of this EJA.
1080 This is a cached method because we know the rank a priori for
1081 all of the algebras we can construct. Thus we can avoid the
1082 expensive ``_charpoly_coefficients()`` call unless we truly
1083 need to compute the whole characteristic polynomial.
1087 sage: from mjo.eja.eja_algebra import (HadamardEJA,
1088 ....: JordanSpinEJA,
1089 ....: RealSymmetricEJA,
1090 ....: ComplexHermitianEJA,
1091 ....: QuaternionHermitianEJA,
1096 The rank of the Jordan spin algebra is always two::
1098 sage: JordanSpinEJA(2).rank()
1100 sage: JordanSpinEJA(3).rank()
1102 sage: JordanSpinEJA(4).rank()
1105 The rank of the `n`-by-`n` Hermitian real, complex, or
1106 quaternion matrices is `n`::
1108 sage: RealSymmetricEJA(4).rank()
1110 sage: ComplexHermitianEJA(3).rank()
1112 sage: QuaternionHermitianEJA(2).rank()
1117 Ensure that every EJA that we know how to construct has a
1118 positive integer rank, unless the algebra is trivial in
1119 which case its rank will be zero::
1121 sage: set_random_seed()
1122 sage: J = random_eja()
1126 sage: r > 0 or (r == 0 and J.is_trivial())
1129 Ensure that computing the rank actually works, since the ranks
1130 of all simple algebras are known and will be cached by default::
1132 sage: set_random_seed() # long time
1133 sage: J = random_eja() # long time
1134 sage: caches = J.rank() # long time
1135 sage: J.rank.clear_cache() # long time
1136 sage: J.rank() == cached # long time
1140 return len(self
._charpoly
_coefficients
())
1143 def vector_space(self
):
1145 Return the vector space that underlies this algebra.
1149 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1153 sage: J = RealSymmetricEJA(2)
1154 sage: J.vector_space()
1155 Vector space of dimension 3 over...
1158 return self
.zero().to_vector().parent().ambient_vector_space()
1161 Element
= FiniteDimensionalEJAElement
1163 class RationalBasisEJA(FiniteDimensionalEJA
):
1165 New class for algebras whose supplied basis elements have all rational entries.
1169 sage: from mjo.eja.eja_algebra import BilinearFormEJA
1173 The supplied basis is orthonormalized by default::
1175 sage: B = matrix(QQ, [[1, 0, 0], [0, 25, -32], [0, -32, 41]])
1176 sage: J = BilinearFormEJA(B)
1177 sage: J.matrix_basis()
1190 orthonormalize
=True,
1196 # Abuse the check_field parameter to check that the entries of
1197 # out basis (in ambient coordinates) are in the field QQ.
1198 if not all( all(b_i
in QQ
for b_i
in b
.list()) for b
in basis
):
1199 raise TypeError("basis not rational")
1202 # There's no point in constructing the extra algebra if this
1203 # one is already rational.
1205 # Note: the same Jordan and inner-products work here,
1206 # because they are necessarily defined with respect to
1207 # ambient coordinates and not any particular basis.
1208 self
._rational
_algebra
= FiniteDimensionalEJA(
1213 orthonormalize
=False,
1218 super().__init
__(basis
,
1222 check_field
=check_field
,
1223 check_axioms
=check_axioms
,
1227 def _charpoly_coefficients(self
):
1231 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
1232 ....: JordanSpinEJA)
1236 The base ring of the resulting polynomial coefficients is what
1237 it should be, and not the rationals (unless the algebra was
1238 already over the rationals)::
1240 sage: J = JordanSpinEJA(3)
1241 sage: J._charpoly_coefficients()
1242 (X1^2 - X2^2 - X3^2, -2*X1)
1243 sage: a0 = J._charpoly_coefficients()[0]
1245 Algebraic Real Field
1246 sage: a0.base_ring()
1247 Algebraic Real Field
1250 if self
._rational
_algebra
is None:
1251 # There's no need to construct *another* algebra over the
1252 # rationals if this one is already over the
1253 # rationals. Likewise, if we never orthonormalized our
1254 # basis, we might as well just use the given one.
1255 return super()._charpoly
_coefficients
()
1257 # Do the computation over the rationals. The answer will be
1258 # the same, because all we've done is a change of basis.
1259 # Then, change back from QQ to our real base ring
1260 a
= ( a_i
.change_ring(self
.base_ring())
1261 for a_i
in self
._rational
_algebra
._charpoly
_coefficients
() )
1263 # Now convert the coordinate variables back to the
1264 # deorthonormalized ones.
1265 R
= self
.coordinate_polynomial_ring()
1266 from sage
.modules
.free_module_element
import vector
1267 X
= vector(R
, R
.gens())
1268 BX
= self
._deortho
_matrix
*X
1270 subs_dict
= { X[i]: BX[i] for i in range(len(X)) }
1271 return tuple( a_i
.subs(subs_dict
) for a_i
in a
)
1273 class ConcreteEJA(RationalBasisEJA
):
1275 A class for the Euclidean Jordan algebras that we know by name.
1277 These are the Jordan algebras whose basis, multiplication table,
1278 rank, and so on are known a priori. More to the point, they are
1279 the Euclidean Jordan algebras for which we are able to conjure up
1280 a "random instance."
1284 sage: from mjo.eja.eja_algebra import ConcreteEJA
1288 Our basis is normalized with respect to the algebra's inner
1289 product, unless we specify otherwise::
1291 sage: set_random_seed()
1292 sage: J = ConcreteEJA.random_instance()
1293 sage: all( b.norm() == 1 for b in J.gens() )
1296 Since our basis is orthonormal with respect to the algebra's inner
1297 product, and since we know that this algebra is an EJA, any
1298 left-multiplication operator's matrix will be symmetric because
1299 natural->EJA basis representation is an isometry and within the
1300 EJA the operator is self-adjoint by the Jordan axiom::
1302 sage: set_random_seed()
1303 sage: J = ConcreteEJA.random_instance()
1304 sage: x = J.random_element()
1305 sage: x.operator().is_self_adjoint()
1310 def _max_random_instance_size():
1312 Return an integer "size" that is an upper bound on the size of
1313 this algebra when it is used in a random test
1314 case. Unfortunately, the term "size" is ambiguous -- when
1315 dealing with `R^n` under either the Hadamard or Jordan spin
1316 product, the "size" refers to the dimension `n`. When dealing
1317 with a matrix algebra (real symmetric or complex/quaternion
1318 Hermitian), it refers to the size of the matrix, which is far
1319 less than the dimension of the underlying vector space.
1321 This method must be implemented in each subclass.
1323 raise NotImplementedError
1326 def random_instance(cls
, *args
, **kwargs
):
1328 Return a random instance of this type of algebra.
1330 This method should be implemented in each subclass.
1332 from sage
.misc
.prandom
import choice
1333 eja_class
= choice(cls
.__subclasses
__())
1335 # These all bubble up to the RationalBasisEJA superclass
1336 # constructor, so any (kw)args valid there are also valid
1338 return eja_class
.random_instance(*args
, **kwargs
)
1343 def dimension_over_reals():
1345 The dimension of this matrix's base ring over the reals.
1347 The reals are dimension one over themselves, obviously; that's
1348 just `\mathbb{R}^{1}`. Likewise, the complex numbers `a + bi`
1349 have dimension two. Finally, the quaternions have dimension
1350 four over the reals.
1352 This is used to determine the size of the matrix returned from
1353 :meth:`real_embed`, among other things.
1355 raise NotImplementedError
1358 def real_embed(cls
,M
):
1360 Embed the matrix ``M`` into a space of real matrices.
1362 The matrix ``M`` can have entries in any field at the moment:
1363 the real numbers, complex numbers, or quaternions. And although
1364 they are not a field, we can probably support octonions at some
1365 point, too. This function returns a real matrix that "acts like"
1366 the original with respect to matrix multiplication; i.e.
1368 real_embed(M*N) = real_embed(M)*real_embed(N)
1371 if M
.ncols() != M
.nrows():
1372 raise ValueError("the matrix 'M' must be square")
1377 def real_unembed(cls
,M
):
1379 The inverse of :meth:`real_embed`.
1381 if M
.ncols() != M
.nrows():
1382 raise ValueError("the matrix 'M' must be square")
1383 if not ZZ(M
.nrows()).mod(cls
.dimension_over_reals()).is_zero():
1384 raise ValueError("the matrix 'M' must be a real embedding")
1388 def jordan_product(X
,Y
):
1389 return (X
*Y
+ Y
*X
)/2
1392 def trace_inner_product(cls
,X
,Y
):
1394 Compute the trace inner-product of two real-embeddings.
1398 sage: from mjo.eja.eja_algebra import (RealSymmetricEJA,
1399 ....: ComplexHermitianEJA,
1400 ....: QuaternionHermitianEJA)
1404 This gives the same answer as it would if we computed the trace
1405 from the unembedded (original) matrices::
1407 sage: set_random_seed()
1408 sage: J = RealSymmetricEJA.random_instance()
1409 sage: x,y = J.random_elements(2)
1410 sage: Xe = x.to_matrix()
1411 sage: Ye = y.to_matrix()
1412 sage: X = J.real_unembed(Xe)
1413 sage: Y = J.real_unembed(Ye)
1414 sage: expected = (X*Y).trace()
1415 sage: actual = J.trace_inner_product(Xe,Ye)
1416 sage: actual == expected
1421 sage: set_random_seed()
1422 sage: J = ComplexHermitianEJA.random_instance()
1423 sage: x,y = J.random_elements(2)
1424 sage: Xe = x.to_matrix()
1425 sage: Ye = y.to_matrix()
1426 sage: X = J.real_unembed(Xe)
1427 sage: Y = J.real_unembed(Ye)
1428 sage: expected = (X*Y).trace().real()
1429 sage: actual = J.trace_inner_product(Xe,Ye)
1430 sage: actual == expected
1435 sage: set_random_seed()
1436 sage: J = QuaternionHermitianEJA.random_instance()
1437 sage: x,y = J.random_elements(2)
1438 sage: Xe = x.to_matrix()
1439 sage: Ye = y.to_matrix()
1440 sage: X = J.real_unembed(Xe)
1441 sage: Y = J.real_unembed(Ye)
1442 sage: expected = (X*Y).trace().coefficient_tuple()[0]
1443 sage: actual = J.trace_inner_product(Xe,Ye)
1444 sage: actual == expected
1448 Xu
= cls
.real_unembed(X
)
1449 Yu
= cls
.real_unembed(Y
)
1450 tr
= (Xu
*Yu
).trace()
1453 # Works in QQ, AA, RDF, et cetera.
1455 except AttributeError:
1456 # A quaternion doesn't have a real() method, but does
1457 # have coefficient_tuple() method that returns the
1458 # coefficients of 1, i, j, and k -- in that order.
1459 return tr
.coefficient_tuple()[0]
1462 class RealMatrixEJA(MatrixEJA
):
1464 def dimension_over_reals():
1468 class RealSymmetricEJA(ConcreteEJA
, RealMatrixEJA
):
1470 The rank-n simple EJA consisting of real symmetric n-by-n
1471 matrices, the usual symmetric Jordan product, and the trace inner
1472 product. It has dimension `(n^2 + n)/2` over the reals.
1476 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1480 sage: J = RealSymmetricEJA(2)
1481 sage: e0, e1, e2 = J.gens()
1489 In theory, our "field" can be any subfield of the reals::
1491 sage: RealSymmetricEJA(2, field=RDF)
1492 Euclidean Jordan algebra of dimension 3 over Real Double Field
1493 sage: RealSymmetricEJA(2, field=RR)
1494 Euclidean Jordan algebra of dimension 3 over Real Field with
1495 53 bits of precision
1499 The dimension of this algebra is `(n^2 + n) / 2`::
1501 sage: set_random_seed()
1502 sage: n_max = RealSymmetricEJA._max_random_instance_size()
1503 sage: n = ZZ.random_element(1, n_max)
1504 sage: J = RealSymmetricEJA(n)
1505 sage: J.dimension() == (n^2 + n)/2
1508 The Jordan multiplication is what we think it is::
1510 sage: set_random_seed()
1511 sage: J = RealSymmetricEJA.random_instance()
1512 sage: x,y = J.random_elements(2)
1513 sage: actual = (x*y).to_matrix()
1514 sage: X = x.to_matrix()
1515 sage: Y = y.to_matrix()
1516 sage: expected = (X*Y + Y*X)/2
1517 sage: actual == expected
1519 sage: J(expected) == x*y
1522 We can change the generator prefix::
1524 sage: RealSymmetricEJA(3, prefix='q').gens()
1525 (q0, q1, q2, q3, q4, q5)
1527 We can construct the (trivial) algebra of rank zero::
1529 sage: RealSymmetricEJA(0)
1530 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1534 def _denormalized_basis(cls
, n
):
1536 Return a basis for the space of real symmetric n-by-n matrices.
1540 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1544 sage: set_random_seed()
1545 sage: n = ZZ.random_element(1,5)
1546 sage: B = RealSymmetricEJA._denormalized_basis(n)
1547 sage: all( M.is_symmetric() for M in B)
1551 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1555 for j
in range(i
+1):
1556 Eij
= matrix(ZZ
, n
, lambda k
,l
: k
==i
and l
==j
)
1560 Sij
= Eij
+ Eij
.transpose()
1566 def _max_random_instance_size():
1567 return 4 # Dimension 10
1570 def random_instance(cls
, **kwargs
):
1572 Return a random instance of this type of algebra.
1574 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1575 return cls(n
, **kwargs
)
1577 def __init__(self
, n
, **kwargs
):
1578 # We know this is a valid EJA, but will double-check
1579 # if the user passes check_axioms=True.
1580 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
1582 super(RealSymmetricEJA
, self
).__init
__(self
._denormalized
_basis
(n
),
1583 self
.jordan_product
,
1584 self
.trace_inner_product
,
1587 # TODO: this could be factored out somehow, but is left here
1588 # because the MatrixEJA is not presently a subclass of the
1589 # FDEJA class that defines rank() and one().
1590 self
.rank
.set_cache(n
)
1591 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
1592 self
.one
.set_cache(self(idV
))
1596 class ComplexMatrixEJA(MatrixEJA
):
1598 def dimension_over_reals():
1602 def real_embed(cls
,M
):
1604 Embed the n-by-n complex matrix ``M`` into the space of real
1605 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1606 bi` to the block matrix ``[[a,b],[-b,a]]``.
1610 sage: from mjo.eja.eja_algebra import ComplexMatrixEJA
1614 sage: F = QuadraticField(-1, 'I')
1615 sage: x1 = F(4 - 2*i)
1616 sage: x2 = F(1 + 2*i)
1619 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1620 sage: ComplexMatrixEJA.real_embed(M)
1629 Embedding is a homomorphism (isomorphism, in fact)::
1631 sage: set_random_seed()
1632 sage: n = ZZ.random_element(3)
1633 sage: F = QuadraticField(-1, 'I')
1634 sage: X = random_matrix(F, n)
1635 sage: Y = random_matrix(F, n)
1636 sage: Xe = ComplexMatrixEJA.real_embed(X)
1637 sage: Ye = ComplexMatrixEJA.real_embed(Y)
1638 sage: XYe = ComplexMatrixEJA.real_embed(X*Y)
1643 super(ComplexMatrixEJA
,cls
).real_embed(M
)
1646 # We don't need any adjoined elements...
1647 field
= M
.base_ring().base_ring()
1651 a
= z
.list()[0] # real part, I guess
1652 b
= z
.list()[1] # imag part, I guess
1653 blocks
.append(matrix(field
, 2, [[a
,b
],[-b
,a
]]))
1655 return matrix
.block(field
, n
, blocks
)
1659 def real_unembed(cls
,M
):
1661 The inverse of _embed_complex_matrix().
1665 sage: from mjo.eja.eja_algebra import ComplexMatrixEJA
1669 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1670 ....: [-2, 1, -4, 3],
1671 ....: [ 9, 10, 11, 12],
1672 ....: [-10, 9, -12, 11] ])
1673 sage: ComplexMatrixEJA.real_unembed(A)
1675 [ 10*I + 9 12*I + 11]
1679 Unembedding is the inverse of embedding::
1681 sage: set_random_seed()
1682 sage: F = QuadraticField(-1, 'I')
1683 sage: M = random_matrix(F, 3)
1684 sage: Me = ComplexMatrixEJA.real_embed(M)
1685 sage: ComplexMatrixEJA.real_unembed(Me) == M
1689 super(ComplexMatrixEJA
,cls
).real_unembed(M
)
1691 d
= cls
.dimension_over_reals()
1693 # If "M" was normalized, its base ring might have roots
1694 # adjoined and they can stick around after unembedding.
1695 field
= M
.base_ring()
1696 R
= PolynomialRing(field
, 'z')
1699 # Sage doesn't know how to adjoin the complex "i" (the root of
1700 # x^2 + 1) to a field in a general way. Here, we just enumerate
1701 # all of the cases that I have cared to support so far.
1703 # Sage doesn't know how to embed AA into QQbar, i.e. how
1704 # to adjoin sqrt(-1) to AA.
1706 elif not field
.is_exact():
1708 F
= field
.complex_field()
1710 # Works for QQ and... maybe some other fields.
1711 F
= field
.extension(z
**2 + 1, 'I', embedding
=CLF(-1).sqrt())
1714 # Go top-left to bottom-right (reading order), converting every
1715 # 2-by-2 block we see to a single complex element.
1717 for k
in range(n
/d
):
1718 for j
in range(n
/d
):
1719 submat
= M
[d
*k
:d
*k
+d
,d
*j
:d
*j
+d
]
1720 if submat
[0,0] != submat
[1,1]:
1721 raise ValueError('bad on-diagonal submatrix')
1722 if submat
[0,1] != -submat
[1,0]:
1723 raise ValueError('bad off-diagonal submatrix')
1724 z
= submat
[0,0] + submat
[0,1]*i
1727 return matrix(F
, n
/d
, elements
)
1730 class ComplexHermitianEJA(ConcreteEJA
, ComplexMatrixEJA
):
1732 The rank-n simple EJA consisting of complex Hermitian n-by-n
1733 matrices over the real numbers, the usual symmetric Jordan product,
1734 and the real-part-of-trace inner product. It has dimension `n^2` over
1739 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1743 In theory, our "field" can be any subfield of the reals::
1745 sage: ComplexHermitianEJA(2, field=RDF)
1746 Euclidean Jordan algebra of dimension 4 over Real Double Field
1747 sage: ComplexHermitianEJA(2, field=RR)
1748 Euclidean Jordan algebra of dimension 4 over Real Field with
1749 53 bits of precision
1753 The dimension of this algebra is `n^2`::
1755 sage: set_random_seed()
1756 sage: n_max = ComplexHermitianEJA._max_random_instance_size()
1757 sage: n = ZZ.random_element(1, n_max)
1758 sage: J = ComplexHermitianEJA(n)
1759 sage: J.dimension() == n^2
1762 The Jordan multiplication is what we think it is::
1764 sage: set_random_seed()
1765 sage: J = ComplexHermitianEJA.random_instance()
1766 sage: x,y = J.random_elements(2)
1767 sage: actual = (x*y).to_matrix()
1768 sage: X = x.to_matrix()
1769 sage: Y = y.to_matrix()
1770 sage: expected = (X*Y + Y*X)/2
1771 sage: actual == expected
1773 sage: J(expected) == x*y
1776 We can change the generator prefix::
1778 sage: ComplexHermitianEJA(2, prefix='z').gens()
1781 We can construct the (trivial) algebra of rank zero::
1783 sage: ComplexHermitianEJA(0)
1784 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1789 def _denormalized_basis(cls
, n
):
1791 Returns a basis for the space of complex Hermitian n-by-n matrices.
1793 Why do we embed these? Basically, because all of numerical linear
1794 algebra assumes that you're working with vectors consisting of `n`
1795 entries from a field and scalars from the same field. There's no way
1796 to tell SageMath that (for example) the vectors contain complex
1797 numbers, while the scalar field is real.
1801 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1805 sage: set_random_seed()
1806 sage: n = ZZ.random_element(1,5)
1807 sage: field = QuadraticField(2, 'sqrt2')
1808 sage: B = ComplexHermitianEJA._denormalized_basis(n)
1809 sage: all( M.is_symmetric() for M in B)
1814 R
= PolynomialRing(field
, 'z')
1816 F
= field
.extension(z
**2 + 1, 'I')
1819 # This is like the symmetric case, but we need to be careful:
1821 # * We want conjugate-symmetry, not just symmetry.
1822 # * The diagonal will (as a result) be real.
1826 for j
in range(i
+1):
1827 Eij
= matrix(F
, n
, lambda k
,l
: k
==i
and l
==j
)
1829 Sij
= cls
.real_embed(Eij
)
1832 # The second one has a minus because it's conjugated.
1833 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
1835 Sij_imag
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
1838 # Since we embedded these, we can drop back to the "field" that we
1839 # started with instead of the complex extension "F".
1840 return tuple( s
.change_ring(field
) for s
in S
)
1843 def __init__(self
, n
, **kwargs
):
1844 # We know this is a valid EJA, but will double-check
1845 # if the user passes check_axioms=True.
1846 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
1848 super(ComplexHermitianEJA
, self
).__init
__(self
._denormalized
_basis
(n
),
1849 self
.jordan_product
,
1850 self
.trace_inner_product
,
1852 # TODO: this could be factored out somehow, but is left here
1853 # because the MatrixEJA is not presently a subclass of the
1854 # FDEJA class that defines rank() and one().
1855 self
.rank
.set_cache(n
)
1856 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
1857 self
.one
.set_cache(self(idV
))
1860 def _max_random_instance_size():
1861 return 3 # Dimension 9
1864 def random_instance(cls
, **kwargs
):
1866 Return a random instance of this type of algebra.
1868 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1869 return cls(n
, **kwargs
)
1871 class QuaternionMatrixEJA(MatrixEJA
):
1873 def dimension_over_reals():
1877 def real_embed(cls
,M
):
1879 Embed the n-by-n quaternion matrix ``M`` into the space of real
1880 matrices of size 4n-by-4n by first sending each quaternion entry `z
1881 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
1882 c+di],[-c + di, a-bi]]`, and then embedding those into a real
1887 sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA
1891 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1892 sage: i,j,k = Q.gens()
1893 sage: x = 1 + 2*i + 3*j + 4*k
1894 sage: M = matrix(Q, 1, [[x]])
1895 sage: QuaternionMatrixEJA.real_embed(M)
1901 Embedding is a homomorphism (isomorphism, in fact)::
1903 sage: set_random_seed()
1904 sage: n = ZZ.random_element(2)
1905 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1906 sage: X = random_matrix(Q, n)
1907 sage: Y = random_matrix(Q, n)
1908 sage: Xe = QuaternionMatrixEJA.real_embed(X)
1909 sage: Ye = QuaternionMatrixEJA.real_embed(Y)
1910 sage: XYe = QuaternionMatrixEJA.real_embed(X*Y)
1915 super(QuaternionMatrixEJA
,cls
).real_embed(M
)
1916 quaternions
= M
.base_ring()
1919 F
= QuadraticField(-1, 'I')
1924 t
= z
.coefficient_tuple()
1929 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
1930 [-c
+ d
*i
, a
- b
*i
]])
1931 realM
= ComplexMatrixEJA
.real_embed(cplxM
)
1932 blocks
.append(realM
)
1934 # We should have real entries by now, so use the realest field
1935 # we've got for the return value.
1936 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
1941 def real_unembed(cls
,M
):
1943 The inverse of _embed_quaternion_matrix().
1947 sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA
1951 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
1952 ....: [-2, 1, -4, 3],
1953 ....: [-3, 4, 1, -2],
1954 ....: [-4, -3, 2, 1]])
1955 sage: QuaternionMatrixEJA.real_unembed(M)
1956 [1 + 2*i + 3*j + 4*k]
1960 Unembedding is the inverse of embedding::
1962 sage: set_random_seed()
1963 sage: Q = QuaternionAlgebra(QQ, -1, -1)
1964 sage: M = random_matrix(Q, 3)
1965 sage: Me = QuaternionMatrixEJA.real_embed(M)
1966 sage: QuaternionMatrixEJA.real_unembed(Me) == M
1970 super(QuaternionMatrixEJA
,cls
).real_unembed(M
)
1972 d
= cls
.dimension_over_reals()
1974 # Use the base ring of the matrix to ensure that its entries can be
1975 # multiplied by elements of the quaternion algebra.
1976 field
= M
.base_ring()
1977 Q
= QuaternionAlgebra(field
,-1,-1)
1980 # Go top-left to bottom-right (reading order), converting every
1981 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
1984 for l
in range(n
/d
):
1985 for m
in range(n
/d
):
1986 submat
= ComplexMatrixEJA
.real_unembed(
1987 M
[d
*l
:d
*l
+d
,d
*m
:d
*m
+d
] )
1988 if submat
[0,0] != submat
[1,1].conjugate():
1989 raise ValueError('bad on-diagonal submatrix')
1990 if submat
[0,1] != -submat
[1,0].conjugate():
1991 raise ValueError('bad off-diagonal submatrix')
1992 z
= submat
[0,0].real()
1993 z
+= submat
[0,0].imag()*i
1994 z
+= submat
[0,1].real()*j
1995 z
+= submat
[0,1].imag()*k
1998 return matrix(Q
, n
/d
, elements
)
2001 class QuaternionHermitianEJA(ConcreteEJA
, QuaternionMatrixEJA
):
2003 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
2004 matrices, the usual symmetric Jordan product, and the
2005 real-part-of-trace inner product. It has dimension `2n^2 - n` over
2010 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2014 In theory, our "field" can be any subfield of the reals::
2016 sage: QuaternionHermitianEJA(2, field=RDF)
2017 Euclidean Jordan algebra of dimension 6 over Real Double Field
2018 sage: QuaternionHermitianEJA(2, field=RR)
2019 Euclidean Jordan algebra of dimension 6 over Real Field with
2020 53 bits of precision
2024 The dimension of this algebra is `2*n^2 - n`::
2026 sage: set_random_seed()
2027 sage: n_max = QuaternionHermitianEJA._max_random_instance_size()
2028 sage: n = ZZ.random_element(1, n_max)
2029 sage: J = QuaternionHermitianEJA(n)
2030 sage: J.dimension() == 2*(n^2) - n
2033 The Jordan multiplication is what we think it is::
2035 sage: set_random_seed()
2036 sage: J = QuaternionHermitianEJA.random_instance()
2037 sage: x,y = J.random_elements(2)
2038 sage: actual = (x*y).to_matrix()
2039 sage: X = x.to_matrix()
2040 sage: Y = y.to_matrix()
2041 sage: expected = (X*Y + Y*X)/2
2042 sage: actual == expected
2044 sage: J(expected) == x*y
2047 We can change the generator prefix::
2049 sage: QuaternionHermitianEJA(2, prefix='a').gens()
2050 (a0, a1, a2, a3, a4, a5)
2052 We can construct the (trivial) algebra of rank zero::
2054 sage: QuaternionHermitianEJA(0)
2055 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2059 def _denormalized_basis(cls
, n
):
2061 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
2063 Why do we embed these? Basically, because all of numerical
2064 linear algebra assumes that you're working with vectors consisting
2065 of `n` entries from a field and scalars from the same field. There's
2066 no way to tell SageMath that (for example) the vectors contain
2067 complex numbers, while the scalar field is real.
2071 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2075 sage: set_random_seed()
2076 sage: n = ZZ.random_element(1,5)
2077 sage: B = QuaternionHermitianEJA._denormalized_basis(n)
2078 sage: all( M.is_symmetric() for M in B )
2083 Q
= QuaternionAlgebra(QQ
,-1,-1)
2086 # This is like the symmetric case, but we need to be careful:
2088 # * We want conjugate-symmetry, not just symmetry.
2089 # * The diagonal will (as a result) be real.
2093 for j
in range(i
+1):
2094 Eij
= matrix(Q
, n
, lambda k
,l
: k
==i
and l
==j
)
2096 Sij
= cls
.real_embed(Eij
)
2099 # The second, third, and fourth ones have a minus
2100 # because they're conjugated.
2101 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
2103 Sij_I
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
2105 Sij_J
= cls
.real_embed(J
*Eij
- J
*Eij
.transpose())
2107 Sij_K
= cls
.real_embed(K
*Eij
- K
*Eij
.transpose())
2110 # Since we embedded these, we can drop back to the "field" that we
2111 # started with instead of the quaternion algebra "Q".
2112 return tuple( s
.change_ring(field
) for s
in S
)
2115 def __init__(self
, n
, **kwargs
):
2116 # We know this is a valid EJA, but will double-check
2117 # if the user passes check_axioms=True.
2118 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2120 super(QuaternionHermitianEJA
, self
).__init
__(self
._denormalized
_basis
(n
),
2121 self
.jordan_product
,
2122 self
.trace_inner_product
,
2124 # TODO: this could be factored out somehow, but is left here
2125 # because the MatrixEJA is not presently a subclass of the
2126 # FDEJA class that defines rank() and one().
2127 self
.rank
.set_cache(n
)
2128 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
2129 self
.one
.set_cache(self(idV
))
2133 def _max_random_instance_size():
2135 The maximum rank of a random QuaternionHermitianEJA.
2137 return 2 # Dimension 6
2140 def random_instance(cls
, **kwargs
):
2142 Return a random instance of this type of algebra.
2144 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2145 return cls(n
, **kwargs
)
2148 class HadamardEJA(ConcreteEJA
):
2150 Return the Euclidean Jordan Algebra corresponding to the set
2151 `R^n` under the Hadamard product.
2153 Note: this is nothing more than the Cartesian product of ``n``
2154 copies of the spin algebra. Once Cartesian product algebras
2155 are implemented, this can go.
2159 sage: from mjo.eja.eja_algebra import HadamardEJA
2163 This multiplication table can be verified by hand::
2165 sage: J = HadamardEJA(3)
2166 sage: e0,e1,e2 = J.gens()
2182 We can change the generator prefix::
2184 sage: HadamardEJA(3, prefix='r').gens()
2188 def __init__(self
, n
, **kwargs
):
2190 jordan_product
= lambda x
,y
: x
2191 inner_product
= lambda x
,y
: x
2193 def jordan_product(x
,y
):
2195 return P( xi
*yi
for (xi
,yi
) in zip(x
,y
) )
2197 def inner_product(x
,y
):
2200 # New defaults for keyword arguments. Don't orthonormalize
2201 # because our basis is already orthonormal with respect to our
2202 # inner-product. Don't check the axioms, because we know this
2203 # is a valid EJA... but do double-check if the user passes
2204 # check_axioms=True. Note: we DON'T override the "check_field"
2205 # default here, because the user can pass in a field!
2206 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2207 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2209 column_basis
= tuple( b
.column() for b
in FreeModule(ZZ
, n
).basis() )
2210 super().__init
__(column_basis
, jordan_product
, inner_product
, **kwargs
)
2211 self
.rank
.set_cache(n
)
2214 self
.one
.set_cache( self
.zero() )
2216 self
.one
.set_cache( sum(self
.gens()) )
2219 def _max_random_instance_size():
2221 The maximum dimension of a random HadamardEJA.
2226 def random_instance(cls
, **kwargs
):
2228 Return a random instance of this type of algebra.
2230 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2231 return cls(n
, **kwargs
)
2234 class BilinearFormEJA(ConcreteEJA
):
2236 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2237 with the half-trace inner product and jordan product ``x*y =
2238 (<Bx,y>,y_bar>, x0*y_bar + y0*x_bar)`` where `B = 1 \times B22` is
2239 a symmetric positive-definite "bilinear form" matrix. Its
2240 dimension is the size of `B`, and it has rank two in dimensions
2241 larger than two. It reduces to the ``JordanSpinEJA`` when `B` is
2242 the identity matrix of order ``n``.
2244 We insist that the one-by-one upper-left identity block of `B` be
2245 passed in as well so that we can be passed a matrix of size zero
2246 to construct a trivial algebra.
2250 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
2251 ....: JordanSpinEJA)
2255 When no bilinear form is specified, the identity matrix is used,
2256 and the resulting algebra is the Jordan spin algebra::
2258 sage: B = matrix.identity(AA,3)
2259 sage: J0 = BilinearFormEJA(B)
2260 sage: J1 = JordanSpinEJA(3)
2261 sage: J0.multiplication_table() == J0.multiplication_table()
2264 An error is raised if the matrix `B` does not correspond to a
2265 positive-definite bilinear form::
2267 sage: B = matrix.random(QQ,2,3)
2268 sage: J = BilinearFormEJA(B)
2269 Traceback (most recent call last):
2271 ValueError: bilinear form is not positive-definite
2272 sage: B = matrix.zero(QQ,3)
2273 sage: J = BilinearFormEJA(B)
2274 Traceback (most recent call last):
2276 ValueError: bilinear form is not positive-definite
2280 We can create a zero-dimensional algebra::
2282 sage: B = matrix.identity(AA,0)
2283 sage: J = BilinearFormEJA(B)
2287 We can check the multiplication condition given in the Jordan, von
2288 Neumann, and Wigner paper (and also discussed on my "On the
2289 symmetry..." paper). Note that this relies heavily on the standard
2290 choice of basis, as does anything utilizing the bilinear form
2291 matrix. We opt not to orthonormalize the basis, because if we
2292 did, we would have to normalize the `s_{i}` in a similar manner::
2294 sage: set_random_seed()
2295 sage: n = ZZ.random_element(5)
2296 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2297 sage: B11 = matrix.identity(QQ,1)
2298 sage: B22 = M.transpose()*M
2299 sage: B = block_matrix(2,2,[ [B11,0 ],
2301 sage: J = BilinearFormEJA(B, orthonormalize=False)
2302 sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
2303 sage: V = J.vector_space()
2304 sage: sis = [ J( V([0] + (M.inverse()*ei).list()).column() )
2305 ....: for ei in eis ]
2306 sage: actual = [ sis[i]*sis[j]
2307 ....: for i in range(n-1)
2308 ....: for j in range(n-1) ]
2309 sage: expected = [ J.one() if i == j else J.zero()
2310 ....: for i in range(n-1)
2311 ....: for j in range(n-1) ]
2312 sage: actual == expected
2315 def __init__(self
, B
, **kwargs
):
2316 if not B
.is_positive_definite():
2317 raise ValueError("bilinear form is not positive-definite")
2319 def inner_product(x
,y
):
2320 return (y
.T
*B
*x
)[0,0]
2322 def jordan_product(x
,y
):
2329 zbar
= y0
*xbar
+ x0
*ybar
2330 return P([0] + zbar
.list())
2332 # We know this is a valid EJA, but will double-check
2333 # if the user passes check_axioms=True.
2334 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2337 column_basis
= tuple( b
.column() for b
in FreeModule(ZZ
, n
).basis() )
2338 super(BilinearFormEJA
, self
).__init
__(column_basis
,
2343 # The rank of this algebra is two, unless we're in a
2344 # one-dimensional ambient space (because the rank is bounded
2345 # by the ambient dimension).
2346 self
.rank
.set_cache(min(n
,2))
2349 self
.one
.set_cache( self
.zero() )
2351 self
.one
.set_cache( self
.monomial(0) )
2354 def _max_random_instance_size():
2356 The maximum dimension of a random BilinearFormEJA.
2361 def random_instance(cls
, **kwargs
):
2363 Return a random instance of this algebra.
2365 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2367 B
= matrix
.identity(ZZ
, n
)
2368 return cls(B
, **kwargs
)
2370 B11
= matrix
.identity(ZZ
, 1)
2371 M
= matrix
.random(ZZ
, n
-1)
2372 I
= matrix
.identity(ZZ
, n
-1)
2374 while alpha
.is_zero():
2375 alpha
= ZZ
.random_element().abs()
2376 B22
= M
.transpose()*M
+ alpha
*I
2378 from sage
.matrix
.special
import block_matrix
2379 B
= block_matrix(2,2, [ [B11
, ZZ(0) ],
2382 return cls(B
, **kwargs
)
2385 class JordanSpinEJA(BilinearFormEJA
):
2387 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2388 with the usual inner product and jordan product ``x*y =
2389 (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2394 sage: from mjo.eja.eja_algebra import JordanSpinEJA
2398 This multiplication table can be verified by hand::
2400 sage: J = JordanSpinEJA(4)
2401 sage: e0,e1,e2,e3 = J.gens()
2417 We can change the generator prefix::
2419 sage: JordanSpinEJA(2, prefix='B').gens()
2424 Ensure that we have the usual inner product on `R^n`::
2426 sage: set_random_seed()
2427 sage: J = JordanSpinEJA.random_instance()
2428 sage: x,y = J.random_elements(2)
2429 sage: actual = x.inner_product(y)
2430 sage: expected = x.to_vector().inner_product(y.to_vector())
2431 sage: actual == expected
2435 def __init__(self
, n
, **kwargs
):
2436 # This is a special case of the BilinearFormEJA with the
2437 # identity matrix as its bilinear form.
2438 B
= matrix
.identity(ZZ
, n
)
2440 # Don't orthonormalize because our basis is already
2441 # orthonormal with respect to our inner-product.
2442 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2444 # But also don't pass check_field=False here, because the user
2445 # can pass in a field!
2446 super(JordanSpinEJA
, self
).__init
__(B
, **kwargs
)
2449 def _max_random_instance_size():
2451 The maximum dimension of a random JordanSpinEJA.
2456 def random_instance(cls
, **kwargs
):
2458 Return a random instance of this type of algebra.
2460 Needed here to override the implementation for ``BilinearFormEJA``.
2462 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2463 return cls(n
, **kwargs
)
2466 class TrivialEJA(ConcreteEJA
):
2468 The trivial Euclidean Jordan algebra consisting of only a zero element.
2472 sage: from mjo.eja.eja_algebra import TrivialEJA
2476 sage: J = TrivialEJA()
2483 sage: 7*J.one()*12*J.one()
2485 sage: J.one().inner_product(J.one())
2487 sage: J.one().norm()
2489 sage: J.one().subalgebra_generated_by()
2490 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2495 def __init__(self
, **kwargs
):
2496 jordan_product
= lambda x
,y
: x
2497 inner_product
= lambda x
,y
: 0
2500 # New defaults for keyword arguments
2501 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2502 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2504 super(TrivialEJA
, self
).__init
__(basis
,
2508 # The rank is zero using my definition, namely the dimension of the
2509 # largest subalgebra generated by any element.
2510 self
.rank
.set_cache(0)
2511 self
.one
.set_cache( self
.zero() )
2514 def random_instance(cls
, **kwargs
):
2515 # We don't take a "size" argument so the superclass method is
2516 # inappropriate for us.
2517 return cls(**kwargs
)
2519 # class DirectSumEJA(ConcreteEJA):
2521 # The external (orthogonal) direct sum of two other Euclidean Jordan
2522 # algebras. Essentially the Cartesian product of its two factors.
2523 # Every Euclidean Jordan algebra decomposes into an orthogonal
2524 # direct sum of simple Euclidean Jordan algebras, so no generality
2525 # is lost by providing only this construction.
2529 # sage: from mjo.eja.eja_algebra import (random_eja,
2530 # ....: HadamardEJA,
2531 # ....: RealSymmetricEJA,
2532 # ....: DirectSumEJA)
2536 # sage: J1 = HadamardEJA(2)
2537 # sage: J2 = RealSymmetricEJA(3)
2538 # sage: J = DirectSumEJA(J1,J2)
2539 # sage: J.dimension()
2546 # The external direct sum construction is only valid when the two factors
2547 # have the same base ring; an error is raised otherwise::
2549 # sage: set_random_seed()
2550 # sage: J1 = random_eja(field=AA)
2551 # sage: J2 = random_eja(field=QQ,orthonormalize=False)
2552 # sage: J = DirectSumEJA(J1,J2)
2553 # Traceback (most recent call last):
2555 # ValueError: algebras must share the same base field
2558 # def __init__(self, J1, J2, **kwargs):
2559 # if J1.base_ring() != J2.base_ring():
2560 # raise ValueError("algebras must share the same base field")
2561 # field = J1.base_ring()
2563 # self._factors = (J1, J2)
2564 # n1 = J1.dimension()
2565 # n2 = J2.dimension()
2567 # V = VectorSpace(field, n)
2568 # mult_table = [ [ V.zero() for j in range(i+1) ]
2569 # for i in range(n) ]
2570 # for i in range(n1):
2571 # for j in range(i+1):
2572 # p = (J1.monomial(i)*J1.monomial(j)).to_vector()
2573 # mult_table[i][j] = V(p.list() + [field.zero()]*n2)
2575 # for i in range(n2):
2576 # for j in range(i+1):
2577 # p = (J2.monomial(i)*J2.monomial(j)).to_vector()
2578 # mult_table[n1+i][n1+j] = V([field.zero()]*n1 + p.list())
2580 # # TODO: build the IP table here from the two constituent IP
2581 # # matrices (it'll be block diagonal, I think).
2582 # ip_table = [ [ field.zero() for j in range(i+1) ]
2583 # for i in range(n) ]
2584 # super(DirectSumEJA, self).__init__(field,
2587 # check_axioms=False,
2589 # self.rank.set_cache(J1.rank() + J2.rank())
2592 # def factors(self):
2594 # Return the pair of this algebra's factors.
2598 # sage: from mjo.eja.eja_algebra import (HadamardEJA,
2599 # ....: JordanSpinEJA,
2600 # ....: DirectSumEJA)
2604 # sage: J1 = HadamardEJA(2, field=QQ)
2605 # sage: J2 = JordanSpinEJA(3, field=QQ)
2606 # sage: J = DirectSumEJA(J1,J2)
2608 # (Euclidean Jordan algebra of dimension 2 over Rational Field,
2609 # Euclidean Jordan algebra of dimension 3 over Rational Field)
2612 # return self._factors
2614 # def projections(self):
2616 # Return a pair of projections onto this algebra's factors.
2620 # sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
2621 # ....: ComplexHermitianEJA,
2622 # ....: DirectSumEJA)
2626 # sage: J1 = JordanSpinEJA(2)
2627 # sage: J2 = ComplexHermitianEJA(2)
2628 # sage: J = DirectSumEJA(J1,J2)
2629 # sage: (pi_left, pi_right) = J.projections()
2630 # sage: J.one().to_vector()
2631 # (1, 0, 1, 0, 0, 1)
2632 # sage: pi_left(J.one()).to_vector()
2634 # sage: pi_right(J.one()).to_vector()
2638 # (J1,J2) = self.factors()
2639 # m = J1.dimension()
2640 # n = J2.dimension()
2641 # V_basis = self.vector_space().basis()
2642 # # Need to specify the dimensions explicitly so that we don't
2643 # # wind up with a zero-by-zero matrix when we want e.g. a
2644 # # zero-by-two matrix (important for composing things).
2645 # P1 = matrix(self.base_ring(), m, m+n, V_basis[:m])
2646 # P2 = matrix(self.base_ring(), n, m+n, V_basis[m:])
2647 # pi_left = FiniteDimensionalEJAOperator(self,J1,P1)
2648 # pi_right = FiniteDimensionalEJAOperator(self,J2,P2)
2649 # return (pi_left, pi_right)
2651 # def inclusions(self):
2653 # Return the pair of inclusion maps from our factors into us.
2657 # sage: from mjo.eja.eja_algebra import (random_eja,
2658 # ....: JordanSpinEJA,
2659 # ....: RealSymmetricEJA,
2660 # ....: DirectSumEJA)
2664 # sage: J1 = JordanSpinEJA(3)
2665 # sage: J2 = RealSymmetricEJA(2)
2666 # sage: J = DirectSumEJA(J1,J2)
2667 # sage: (iota_left, iota_right) = J.inclusions()
2668 # sage: iota_left(J1.zero()) == J.zero()
2670 # sage: iota_right(J2.zero()) == J.zero()
2672 # sage: J1.one().to_vector()
2674 # sage: iota_left(J1.one()).to_vector()
2675 # (1, 0, 0, 0, 0, 0)
2676 # sage: J2.one().to_vector()
2678 # sage: iota_right(J2.one()).to_vector()
2679 # (0, 0, 0, 1, 0, 1)
2680 # sage: J.one().to_vector()
2681 # (1, 0, 0, 1, 0, 1)
2685 # Composing a projection with the corresponding inclusion should
2686 # produce the identity map, and mismatching them should produce
2689 # sage: set_random_seed()
2690 # sage: J1 = random_eja()
2691 # sage: J2 = random_eja()
2692 # sage: J = DirectSumEJA(J1,J2)
2693 # sage: (iota_left, iota_right) = J.inclusions()
2694 # sage: (pi_left, pi_right) = J.projections()
2695 # sage: pi_left*iota_left == J1.one().operator()
2697 # sage: pi_right*iota_right == J2.one().operator()
2699 # sage: (pi_left*iota_right).is_zero()
2701 # sage: (pi_right*iota_left).is_zero()
2705 # (J1,J2) = self.factors()
2706 # m = J1.dimension()
2707 # n = J2.dimension()
2708 # V_basis = self.vector_space().basis()
2709 # # Need to specify the dimensions explicitly so that we don't
2710 # # wind up with a zero-by-zero matrix when we want e.g. a
2711 # # two-by-zero matrix (important for composing things).
2712 # I1 = matrix.column(self.base_ring(), m, m+n, V_basis[:m])
2713 # I2 = matrix.column(self.base_ring(), n, m+n, V_basis[m:])
2714 # iota_left = FiniteDimensionalEJAOperator(J1,self,I1)
2715 # iota_right = FiniteDimensionalEJAOperator(J2,self,I2)
2716 # return (iota_left, iota_right)
2718 # def inner_product(self, x, y):
2720 # The standard Cartesian inner-product.
2722 # We project ``x`` and ``y`` onto our factors, and add up the
2723 # inner-products from the subalgebras.
2728 # sage: from mjo.eja.eja_algebra import (HadamardEJA,
2729 # ....: QuaternionHermitianEJA,
2730 # ....: DirectSumEJA)
2734 # sage: J1 = HadamardEJA(3,field=QQ)
2735 # sage: J2 = QuaternionHermitianEJA(2,field=QQ,orthonormalize=False)
2736 # sage: J = DirectSumEJA(J1,J2)
2737 # sage: x1 = J1.one()
2739 # sage: y1 = J2.one()
2741 # sage: x1.inner_product(x2)
2743 # sage: y1.inner_product(y2)
2745 # sage: J.one().inner_product(J.one())
2749 # (pi_left, pi_right) = self.projections()
2755 # return (x1.inner_product(y1) + x2.inner_product(y2))
2759 random_eja
= ConcreteEJA
.random_instance