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 Jordan or inner-product.
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()
1194 # Abuse the check_field parameter to check that the entries of
1195 # out basis (in ambient coordinates) are in the field QQ.
1196 if not all( all(b_i
in QQ
for b_i
in b
.list()) for b
in basis
):
1197 raise TypeError("basis not rational")
1200 # There's no point in constructing the extra algebra if this
1201 # one is already rational.
1203 # Note: the same Jordan and inner-products work here,
1204 # because they are necessarily defined with respect to
1205 # ambient coordinates and not any particular basis.
1206 self
._rational
_algebra
= FiniteDimensionalEJA(
1211 orthonormalize
=False,
1215 super().__init
__(basis
,
1219 check_field
=check_field
,
1223 def _charpoly_coefficients(self
):
1227 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
1228 ....: JordanSpinEJA)
1232 The base ring of the resulting polynomial coefficients is what
1233 it should be, and not the rationals (unless the algebra was
1234 already over the rationals)::
1236 sage: J = JordanSpinEJA(3)
1237 sage: J._charpoly_coefficients()
1238 (X1^2 - X2^2 - X3^2, -2*X1)
1239 sage: a0 = J._charpoly_coefficients()[0]
1241 Algebraic Real Field
1242 sage: a0.base_ring()
1243 Algebraic Real Field
1246 if self
._rational
_algebra
is None:
1247 # There's no need to construct *another* algebra over the
1248 # rationals if this one is already over the
1249 # rationals. Likewise, if we never orthonormalized our
1250 # basis, we might as well just use the given one.
1251 return super()._charpoly
_coefficients
()
1253 # Do the computation over the rationals. The answer will be
1254 # the same, because all we've done is a change of basis.
1255 # Then, change back from QQ to our real base ring
1256 a
= ( a_i
.change_ring(self
.base_ring())
1257 for a_i
in self
._rational
_algebra
._charpoly
_coefficients
() )
1259 # Now convert the coordinate variables back to the
1260 # deorthonormalized ones.
1261 R
= self
.coordinate_polynomial_ring()
1262 from sage
.modules
.free_module_element
import vector
1263 X
= vector(R
, R
.gens())
1264 BX
= self
._deortho
_matrix
*X
1266 subs_dict
= { X[i]: BX[i] for i in range(len(X)) }
1267 return tuple( a_i
.subs(subs_dict
) for a_i
in a
)
1269 class ConcreteEJA(RationalBasisEJA
):
1271 A class for the Euclidean Jordan algebras that we know by name.
1273 These are the Jordan algebras whose basis, multiplication table,
1274 rank, and so on are known a priori. More to the point, they are
1275 the Euclidean Jordan algebras for which we are able to conjure up
1276 a "random instance."
1280 sage: from mjo.eja.eja_algebra import ConcreteEJA
1284 Our basis is normalized with respect to the algebra's inner
1285 product, unless we specify otherwise::
1287 sage: set_random_seed()
1288 sage: J = ConcreteEJA.random_instance()
1289 sage: all( b.norm() == 1 for b in J.gens() )
1292 Since our basis is orthonormal with respect to the algebra's inner
1293 product, and since we know that this algebra is an EJA, any
1294 left-multiplication operator's matrix will be symmetric because
1295 natural->EJA basis representation is an isometry and within the
1296 EJA the operator is self-adjoint by the Jordan axiom::
1298 sage: set_random_seed()
1299 sage: J = ConcreteEJA.random_instance()
1300 sage: x = J.random_element()
1301 sage: x.operator().is_self_adjoint()
1306 def _max_random_instance_size():
1308 Return an integer "size" that is an upper bound on the size of
1309 this algebra when it is used in a random test
1310 case. Unfortunately, the term "size" is ambiguous -- when
1311 dealing with `R^n` under either the Hadamard or Jordan spin
1312 product, the "size" refers to the dimension `n`. When dealing
1313 with a matrix algebra (real symmetric or complex/quaternion
1314 Hermitian), it refers to the size of the matrix, which is far
1315 less than the dimension of the underlying vector space.
1317 This method must be implemented in each subclass.
1319 raise NotImplementedError
1322 def random_instance(cls
, *args
, **kwargs
):
1324 Return a random instance of this type of algebra.
1326 This method should be implemented in each subclass.
1328 from sage
.misc
.prandom
import choice
1329 eja_class
= choice(cls
.__subclasses
__())
1331 # These all bubble up to the RationalBasisEJA superclass
1332 # constructor, so any (kw)args valid there are also valid
1334 return eja_class
.random_instance(*args
, **kwargs
)
1339 def dimension_over_reals():
1341 The dimension of this matrix's base ring over the reals.
1343 The reals are dimension one over themselves, obviously; that's
1344 just `\mathbb{R}^{1}`. Likewise, the complex numbers `a + bi`
1345 have dimension two. Finally, the quaternions have dimension
1346 four over the reals.
1348 This is used to determine the size of the matrix returned from
1349 :meth:`real_embed`, among other things.
1351 raise NotImplementedError
1354 def real_embed(cls
,M
):
1356 Embed the matrix ``M`` into a space of real matrices.
1358 The matrix ``M`` can have entries in any field at the moment:
1359 the real numbers, complex numbers, or quaternions. And although
1360 they are not a field, we can probably support octonions at some
1361 point, too. This function returns a real matrix that "acts like"
1362 the original with respect to matrix multiplication; i.e.
1364 real_embed(M*N) = real_embed(M)*real_embed(N)
1367 if M
.ncols() != M
.nrows():
1368 raise ValueError("the matrix 'M' must be square")
1373 def real_unembed(cls
,M
):
1375 The inverse of :meth:`real_embed`.
1377 if M
.ncols() != M
.nrows():
1378 raise ValueError("the matrix 'M' must be square")
1379 if not ZZ(M
.nrows()).mod(cls
.dimension_over_reals()).is_zero():
1380 raise ValueError("the matrix 'M' must be a real embedding")
1384 def jordan_product(X
,Y
):
1385 return (X
*Y
+ Y
*X
)/2
1388 def trace_inner_product(cls
,X
,Y
):
1390 Compute the trace inner-product of two real-embeddings.
1394 sage: from mjo.eja.eja_algebra import (RealSymmetricEJA,
1395 ....: ComplexHermitianEJA,
1396 ....: QuaternionHermitianEJA)
1400 This gives the same answer as it would if we computed the trace
1401 from the unembedded (original) matrices::
1403 sage: set_random_seed()
1404 sage: J = RealSymmetricEJA.random_instance()
1405 sage: x,y = J.random_elements(2)
1406 sage: Xe = x.to_matrix()
1407 sage: Ye = y.to_matrix()
1408 sage: X = J.real_unembed(Xe)
1409 sage: Y = J.real_unembed(Ye)
1410 sage: expected = (X*Y).trace()
1411 sage: actual = J.trace_inner_product(Xe,Ye)
1412 sage: actual == expected
1417 sage: set_random_seed()
1418 sage: J = ComplexHermitianEJA.random_instance()
1419 sage: x,y = J.random_elements(2)
1420 sage: Xe = x.to_matrix()
1421 sage: Ye = y.to_matrix()
1422 sage: X = J.real_unembed(Xe)
1423 sage: Y = J.real_unembed(Ye)
1424 sage: expected = (X*Y).trace().real()
1425 sage: actual = J.trace_inner_product(Xe,Ye)
1426 sage: actual == expected
1431 sage: set_random_seed()
1432 sage: J = QuaternionHermitianEJA.random_instance()
1433 sage: x,y = J.random_elements(2)
1434 sage: Xe = x.to_matrix()
1435 sage: Ye = y.to_matrix()
1436 sage: X = J.real_unembed(Xe)
1437 sage: Y = J.real_unembed(Ye)
1438 sage: expected = (X*Y).trace().coefficient_tuple()[0]
1439 sage: actual = J.trace_inner_product(Xe,Ye)
1440 sage: actual == expected
1444 Xu
= cls
.real_unembed(X
)
1445 Yu
= cls
.real_unembed(Y
)
1446 tr
= (Xu
*Yu
).trace()
1449 # Works in QQ, AA, RDF, et cetera.
1451 except AttributeError:
1452 # A quaternion doesn't have a real() method, but does
1453 # have coefficient_tuple() method that returns the
1454 # coefficients of 1, i, j, and k -- in that order.
1455 return tr
.coefficient_tuple()[0]
1458 class RealMatrixEJA(MatrixEJA
):
1460 def dimension_over_reals():
1464 class RealSymmetricEJA(ConcreteEJA
, RealMatrixEJA
):
1466 The rank-n simple EJA consisting of real symmetric n-by-n
1467 matrices, the usual symmetric Jordan product, and the trace inner
1468 product. It has dimension `(n^2 + n)/2` over the reals.
1472 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1476 sage: J = RealSymmetricEJA(2)
1477 sage: e0, e1, e2 = J.gens()
1485 In theory, our "field" can be any subfield of the reals::
1487 sage: RealSymmetricEJA(2, field=RDF)
1488 Euclidean Jordan algebra of dimension 3 over Real Double Field
1489 sage: RealSymmetricEJA(2, field=RR)
1490 Euclidean Jordan algebra of dimension 3 over Real Field with
1491 53 bits of precision
1495 The dimension of this algebra is `(n^2 + n) / 2`::
1497 sage: set_random_seed()
1498 sage: n_max = RealSymmetricEJA._max_random_instance_size()
1499 sage: n = ZZ.random_element(1, n_max)
1500 sage: J = RealSymmetricEJA(n)
1501 sage: J.dimension() == (n^2 + n)/2
1504 The Jordan multiplication is what we think it is::
1506 sage: set_random_seed()
1507 sage: J = RealSymmetricEJA.random_instance()
1508 sage: x,y = J.random_elements(2)
1509 sage: actual = (x*y).to_matrix()
1510 sage: X = x.to_matrix()
1511 sage: Y = y.to_matrix()
1512 sage: expected = (X*Y + Y*X)/2
1513 sage: actual == expected
1515 sage: J(expected) == x*y
1518 We can change the generator prefix::
1520 sage: RealSymmetricEJA(3, prefix='q').gens()
1521 (q0, q1, q2, q3, q4, q5)
1523 We can construct the (trivial) algebra of rank zero::
1525 sage: RealSymmetricEJA(0)
1526 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1530 def _denormalized_basis(cls
, n
):
1532 Return a basis for the space of real symmetric n-by-n matrices.
1536 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1540 sage: set_random_seed()
1541 sage: n = ZZ.random_element(1,5)
1542 sage: B = RealSymmetricEJA._denormalized_basis(n)
1543 sage: all( M.is_symmetric() for M in B)
1547 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1551 for j
in range(i
+1):
1552 Eij
= matrix(ZZ
, n
, lambda k
,l
: k
==i
and l
==j
)
1556 Sij
= Eij
+ Eij
.transpose()
1562 def _max_random_instance_size():
1563 return 4 # Dimension 10
1566 def random_instance(cls
, **kwargs
):
1568 Return a random instance of this type of algebra.
1570 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1571 return cls(n
, **kwargs
)
1573 def __init__(self
, n
, **kwargs
):
1574 # We know this is a valid EJA, but will double-check
1575 # if the user passes check_axioms=True.
1576 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
1578 super(RealSymmetricEJA
, self
).__init
__(self
._denormalized
_basis
(n
),
1579 self
.jordan_product
,
1580 self
.trace_inner_product
,
1583 # TODO: this could be factored out somehow, but is left here
1584 # because the MatrixEJA is not presently a subclass of the
1585 # FDEJA class that defines rank() and one().
1586 self
.rank
.set_cache(n
)
1587 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
1588 self
.one
.set_cache(self(idV
))
1592 class ComplexMatrixEJA(MatrixEJA
):
1594 def dimension_over_reals():
1598 def real_embed(cls
,M
):
1600 Embed the n-by-n complex matrix ``M`` into the space of real
1601 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1602 bi` to the block matrix ``[[a,b],[-b,a]]``.
1606 sage: from mjo.eja.eja_algebra import ComplexMatrixEJA
1610 sage: F = QuadraticField(-1, 'I')
1611 sage: x1 = F(4 - 2*i)
1612 sage: x2 = F(1 + 2*i)
1615 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1616 sage: ComplexMatrixEJA.real_embed(M)
1625 Embedding is a homomorphism (isomorphism, in fact)::
1627 sage: set_random_seed()
1628 sage: n = ZZ.random_element(3)
1629 sage: F = QuadraticField(-1, 'I')
1630 sage: X = random_matrix(F, n)
1631 sage: Y = random_matrix(F, n)
1632 sage: Xe = ComplexMatrixEJA.real_embed(X)
1633 sage: Ye = ComplexMatrixEJA.real_embed(Y)
1634 sage: XYe = ComplexMatrixEJA.real_embed(X*Y)
1639 super(ComplexMatrixEJA
,cls
).real_embed(M
)
1642 # We don't need any adjoined elements...
1643 field
= M
.base_ring().base_ring()
1647 a
= z
.list()[0] # real part, I guess
1648 b
= z
.list()[1] # imag part, I guess
1649 blocks
.append(matrix(field
, 2, [[a
,b
],[-b
,a
]]))
1651 return matrix
.block(field
, n
, blocks
)
1655 def real_unembed(cls
,M
):
1657 The inverse of _embed_complex_matrix().
1661 sage: from mjo.eja.eja_algebra import ComplexMatrixEJA
1665 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1666 ....: [-2, 1, -4, 3],
1667 ....: [ 9, 10, 11, 12],
1668 ....: [-10, 9, -12, 11] ])
1669 sage: ComplexMatrixEJA.real_unembed(A)
1671 [ 10*I + 9 12*I + 11]
1675 Unembedding is the inverse of embedding::
1677 sage: set_random_seed()
1678 sage: F = QuadraticField(-1, 'I')
1679 sage: M = random_matrix(F, 3)
1680 sage: Me = ComplexMatrixEJA.real_embed(M)
1681 sage: ComplexMatrixEJA.real_unembed(Me) == M
1685 super(ComplexMatrixEJA
,cls
).real_unembed(M
)
1687 d
= cls
.dimension_over_reals()
1689 # If "M" was normalized, its base ring might have roots
1690 # adjoined and they can stick around after unembedding.
1691 field
= M
.base_ring()
1692 R
= PolynomialRing(field
, 'z')
1695 # Sage doesn't know how to adjoin the complex "i" (the root of
1696 # x^2 + 1) to a field in a general way. Here, we just enumerate
1697 # all of the cases that I have cared to support so far.
1699 # Sage doesn't know how to embed AA into QQbar, i.e. how
1700 # to adjoin sqrt(-1) to AA.
1702 elif not field
.is_exact():
1704 F
= field
.complex_field()
1706 # Works for QQ and... maybe some other fields.
1707 F
= field
.extension(z
**2 + 1, 'I', embedding
=CLF(-1).sqrt())
1710 # Go top-left to bottom-right (reading order), converting every
1711 # 2-by-2 block we see to a single complex element.
1713 for k
in range(n
/d
):
1714 for j
in range(n
/d
):
1715 submat
= M
[d
*k
:d
*k
+d
,d
*j
:d
*j
+d
]
1716 if submat
[0,0] != submat
[1,1]:
1717 raise ValueError('bad on-diagonal submatrix')
1718 if submat
[0,1] != -submat
[1,0]:
1719 raise ValueError('bad off-diagonal submatrix')
1720 z
= submat
[0,0] + submat
[0,1]*i
1723 return matrix(F
, n
/d
, elements
)
1726 class ComplexHermitianEJA(ConcreteEJA
, ComplexMatrixEJA
):
1728 The rank-n simple EJA consisting of complex Hermitian n-by-n
1729 matrices over the real numbers, the usual symmetric Jordan product,
1730 and the real-part-of-trace inner product. It has dimension `n^2` over
1735 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1739 In theory, our "field" can be any subfield of the reals::
1741 sage: ComplexHermitianEJA(2, field=RDF)
1742 Euclidean Jordan algebra of dimension 4 over Real Double Field
1743 sage: ComplexHermitianEJA(2, field=RR)
1744 Euclidean Jordan algebra of dimension 4 over Real Field with
1745 53 bits of precision
1749 The dimension of this algebra is `n^2`::
1751 sage: set_random_seed()
1752 sage: n_max = ComplexHermitianEJA._max_random_instance_size()
1753 sage: n = ZZ.random_element(1, n_max)
1754 sage: J = ComplexHermitianEJA(n)
1755 sage: J.dimension() == n^2
1758 The Jordan multiplication is what we think it is::
1760 sage: set_random_seed()
1761 sage: J = ComplexHermitianEJA.random_instance()
1762 sage: x,y = J.random_elements(2)
1763 sage: actual = (x*y).to_matrix()
1764 sage: X = x.to_matrix()
1765 sage: Y = y.to_matrix()
1766 sage: expected = (X*Y + Y*X)/2
1767 sage: actual == expected
1769 sage: J(expected) == x*y
1772 We can change the generator prefix::
1774 sage: ComplexHermitianEJA(2, prefix='z').gens()
1777 We can construct the (trivial) algebra of rank zero::
1779 sage: ComplexHermitianEJA(0)
1780 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1785 def _denormalized_basis(cls
, n
):
1787 Returns a basis for the space of complex Hermitian n-by-n matrices.
1789 Why do we embed these? Basically, because all of numerical linear
1790 algebra assumes that you're working with vectors consisting of `n`
1791 entries from a field and scalars from the same field. There's no way
1792 to tell SageMath that (for example) the vectors contain complex
1793 numbers, while the scalar field is real.
1797 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1801 sage: set_random_seed()
1802 sage: n = ZZ.random_element(1,5)
1803 sage: field = QuadraticField(2, 'sqrt2')
1804 sage: B = ComplexHermitianEJA._denormalized_basis(n)
1805 sage: all( M.is_symmetric() for M in B)
1810 R
= PolynomialRing(field
, 'z')
1812 F
= field
.extension(z
**2 + 1, 'I')
1815 # This is like the symmetric case, but we need to be careful:
1817 # * We want conjugate-symmetry, not just symmetry.
1818 # * The diagonal will (as a result) be real.
1822 for j
in range(i
+1):
1823 Eij
= matrix(F
, n
, lambda k
,l
: k
==i
and l
==j
)
1825 Sij
= cls
.real_embed(Eij
)
1828 # The second one has a minus because it's conjugated.
1829 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
1831 Sij_imag
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
1834 # Since we embedded these, we can drop back to the "field" that we
1835 # started with instead of the complex extension "F".
1836 return tuple( s
.change_ring(field
) for s
in S
)
1839 def __init__(self
, n
, **kwargs
):
1840 # We know this is a valid EJA, but will double-check
1841 # if the user passes check_axioms=True.
1842 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
1844 super(ComplexHermitianEJA
, self
).__init
__(self
._denormalized
_basis
(n
),
1845 self
.jordan_product
,
1846 self
.trace_inner_product
,
1848 # TODO: this could be factored out somehow, but is left here
1849 # because the MatrixEJA is not presently a subclass of the
1850 # FDEJA class that defines rank() and one().
1851 self
.rank
.set_cache(n
)
1852 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
1853 self
.one
.set_cache(self(idV
))
1856 def _max_random_instance_size():
1857 return 3 # Dimension 9
1860 def random_instance(cls
, **kwargs
):
1862 Return a random instance of this type of algebra.
1864 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
1865 return cls(n
, **kwargs
)
1867 class QuaternionMatrixEJA(MatrixEJA
):
1869 def dimension_over_reals():
1873 def real_embed(cls
,M
):
1875 Embed the n-by-n quaternion matrix ``M`` into the space of real
1876 matrices of size 4n-by-4n by first sending each quaternion entry `z
1877 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
1878 c+di],[-c + di, a-bi]]`, and then embedding those into a real
1883 sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA
1887 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1888 sage: i,j,k = Q.gens()
1889 sage: x = 1 + 2*i + 3*j + 4*k
1890 sage: M = matrix(Q, 1, [[x]])
1891 sage: QuaternionMatrixEJA.real_embed(M)
1897 Embedding is a homomorphism (isomorphism, in fact)::
1899 sage: set_random_seed()
1900 sage: n = ZZ.random_element(2)
1901 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1902 sage: X = random_matrix(Q, n)
1903 sage: Y = random_matrix(Q, n)
1904 sage: Xe = QuaternionMatrixEJA.real_embed(X)
1905 sage: Ye = QuaternionMatrixEJA.real_embed(Y)
1906 sage: XYe = QuaternionMatrixEJA.real_embed(X*Y)
1911 super(QuaternionMatrixEJA
,cls
).real_embed(M
)
1912 quaternions
= M
.base_ring()
1915 F
= QuadraticField(-1, 'I')
1920 t
= z
.coefficient_tuple()
1925 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
1926 [-c
+ d
*i
, a
- b
*i
]])
1927 realM
= ComplexMatrixEJA
.real_embed(cplxM
)
1928 blocks
.append(realM
)
1930 # We should have real entries by now, so use the realest field
1931 # we've got for the return value.
1932 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
1937 def real_unembed(cls
,M
):
1939 The inverse of _embed_quaternion_matrix().
1943 sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA
1947 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
1948 ....: [-2, 1, -4, 3],
1949 ....: [-3, 4, 1, -2],
1950 ....: [-4, -3, 2, 1]])
1951 sage: QuaternionMatrixEJA.real_unembed(M)
1952 [1 + 2*i + 3*j + 4*k]
1956 Unembedding is the inverse of embedding::
1958 sage: set_random_seed()
1959 sage: Q = QuaternionAlgebra(QQ, -1, -1)
1960 sage: M = random_matrix(Q, 3)
1961 sage: Me = QuaternionMatrixEJA.real_embed(M)
1962 sage: QuaternionMatrixEJA.real_unembed(Me) == M
1966 super(QuaternionMatrixEJA
,cls
).real_unembed(M
)
1968 d
= cls
.dimension_over_reals()
1970 # Use the base ring of the matrix to ensure that its entries can be
1971 # multiplied by elements of the quaternion algebra.
1972 field
= M
.base_ring()
1973 Q
= QuaternionAlgebra(field
,-1,-1)
1976 # Go top-left to bottom-right (reading order), converting every
1977 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
1980 for l
in range(n
/d
):
1981 for m
in range(n
/d
):
1982 submat
= ComplexMatrixEJA
.real_unembed(
1983 M
[d
*l
:d
*l
+d
,d
*m
:d
*m
+d
] )
1984 if submat
[0,0] != submat
[1,1].conjugate():
1985 raise ValueError('bad on-diagonal submatrix')
1986 if submat
[0,1] != -submat
[1,0].conjugate():
1987 raise ValueError('bad off-diagonal submatrix')
1988 z
= submat
[0,0].real()
1989 z
+= submat
[0,0].imag()*i
1990 z
+= submat
[0,1].real()*j
1991 z
+= submat
[0,1].imag()*k
1994 return matrix(Q
, n
/d
, elements
)
1997 class QuaternionHermitianEJA(ConcreteEJA
, QuaternionMatrixEJA
):
1999 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
2000 matrices, the usual symmetric Jordan product, and the
2001 real-part-of-trace inner product. It has dimension `2n^2 - n` over
2006 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2010 In theory, our "field" can be any subfield of the reals::
2012 sage: QuaternionHermitianEJA(2, field=RDF)
2013 Euclidean Jordan algebra of dimension 6 over Real Double Field
2014 sage: QuaternionHermitianEJA(2, field=RR)
2015 Euclidean Jordan algebra of dimension 6 over Real Field with
2016 53 bits of precision
2020 The dimension of this algebra is `2*n^2 - n`::
2022 sage: set_random_seed()
2023 sage: n_max = QuaternionHermitianEJA._max_random_instance_size()
2024 sage: n = ZZ.random_element(1, n_max)
2025 sage: J = QuaternionHermitianEJA(n)
2026 sage: J.dimension() == 2*(n^2) - n
2029 The Jordan multiplication is what we think it is::
2031 sage: set_random_seed()
2032 sage: J = QuaternionHermitianEJA.random_instance()
2033 sage: x,y = J.random_elements(2)
2034 sage: actual = (x*y).to_matrix()
2035 sage: X = x.to_matrix()
2036 sage: Y = y.to_matrix()
2037 sage: expected = (X*Y + Y*X)/2
2038 sage: actual == expected
2040 sage: J(expected) == x*y
2043 We can change the generator prefix::
2045 sage: QuaternionHermitianEJA(2, prefix='a').gens()
2046 (a0, a1, a2, a3, a4, a5)
2048 We can construct the (trivial) algebra of rank zero::
2050 sage: QuaternionHermitianEJA(0)
2051 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2055 def _denormalized_basis(cls
, n
):
2057 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
2059 Why do we embed these? Basically, because all of numerical
2060 linear algebra assumes that you're working with vectors consisting
2061 of `n` entries from a field and scalars from the same field. There's
2062 no way to tell SageMath that (for example) the vectors contain
2063 complex numbers, while the scalar field is real.
2067 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2071 sage: set_random_seed()
2072 sage: n = ZZ.random_element(1,5)
2073 sage: B = QuaternionHermitianEJA._denormalized_basis(n)
2074 sage: all( M.is_symmetric() for M in B )
2079 Q
= QuaternionAlgebra(QQ
,-1,-1)
2082 # This is like the symmetric case, but we need to be careful:
2084 # * We want conjugate-symmetry, not just symmetry.
2085 # * The diagonal will (as a result) be real.
2089 for j
in range(i
+1):
2090 Eij
= matrix(Q
, n
, lambda k
,l
: k
==i
and l
==j
)
2092 Sij
= cls
.real_embed(Eij
)
2095 # The second, third, and fourth ones have a minus
2096 # because they're conjugated.
2097 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
2099 Sij_I
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
2101 Sij_J
= cls
.real_embed(J
*Eij
- J
*Eij
.transpose())
2103 Sij_K
= cls
.real_embed(K
*Eij
- K
*Eij
.transpose())
2106 # Since we embedded these, we can drop back to the "field" that we
2107 # started with instead of the quaternion algebra "Q".
2108 return tuple( s
.change_ring(field
) for s
in S
)
2111 def __init__(self
, n
, **kwargs
):
2112 # We know this is a valid EJA, but will double-check
2113 # if the user passes check_axioms=True.
2114 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2116 super(QuaternionHermitianEJA
, self
).__init
__(self
._denormalized
_basis
(n
),
2117 self
.jordan_product
,
2118 self
.trace_inner_product
,
2120 # TODO: this could be factored out somehow, but is left here
2121 # because the MatrixEJA is not presently a subclass of the
2122 # FDEJA class that defines rank() and one().
2123 self
.rank
.set_cache(n
)
2124 idV
= matrix
.identity(ZZ
, self
.dimension_over_reals()*n
)
2125 self
.one
.set_cache(self(idV
))
2129 def _max_random_instance_size():
2131 The maximum rank of a random QuaternionHermitianEJA.
2133 return 2 # Dimension 6
2136 def random_instance(cls
, **kwargs
):
2138 Return a random instance of this type of algebra.
2140 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2141 return cls(n
, **kwargs
)
2144 class HadamardEJA(ConcreteEJA
):
2146 Return the Euclidean Jordan Algebra corresponding to the set
2147 `R^n` under the Hadamard product.
2149 Note: this is nothing more than the Cartesian product of ``n``
2150 copies of the spin algebra. Once Cartesian product algebras
2151 are implemented, this can go.
2155 sage: from mjo.eja.eja_algebra import HadamardEJA
2159 This multiplication table can be verified by hand::
2161 sage: J = HadamardEJA(3)
2162 sage: e0,e1,e2 = J.gens()
2178 We can change the generator prefix::
2180 sage: HadamardEJA(3, prefix='r').gens()
2184 def __init__(self
, n
, **kwargs
):
2186 jordan_product
= lambda x
,y
: x
2187 inner_product
= lambda x
,y
: x
2189 def jordan_product(x
,y
):
2191 return P( xi
*yi
for (xi
,yi
) in zip(x
,y
) )
2193 def inner_product(x
,y
):
2196 # New defaults for keyword arguments. Don't orthonormalize
2197 # because our basis is already orthonormal with respect to our
2198 # inner-product. Don't check the axioms, because we know this
2199 # is a valid EJA... but do double-check if the user passes
2200 # check_axioms=True. Note: we DON'T override the "check_field"
2201 # default here, because the user can pass in a field!
2202 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2203 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2205 column_basis
= tuple( b
.column() for b
in FreeModule(ZZ
, n
).basis() )
2206 super().__init
__(column_basis
, jordan_product
, inner_product
, **kwargs
)
2207 self
.rank
.set_cache(n
)
2210 self
.one
.set_cache( self
.zero() )
2212 self
.one
.set_cache( sum(self
.gens()) )
2215 def _max_random_instance_size():
2217 The maximum dimension of a random HadamardEJA.
2222 def random_instance(cls
, **kwargs
):
2224 Return a random instance of this type of algebra.
2226 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2227 return cls(n
, **kwargs
)
2230 class BilinearFormEJA(ConcreteEJA
):
2232 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2233 with the half-trace inner product and jordan product ``x*y =
2234 (<Bx,y>,y_bar>, x0*y_bar + y0*x_bar)`` where `B = 1 \times B22` is
2235 a symmetric positive-definite "bilinear form" matrix. Its
2236 dimension is the size of `B`, and it has rank two in dimensions
2237 larger than two. It reduces to the ``JordanSpinEJA`` when `B` is
2238 the identity matrix of order ``n``.
2240 We insist that the one-by-one upper-left identity block of `B` be
2241 passed in as well so that we can be passed a matrix of size zero
2242 to construct a trivial algebra.
2246 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
2247 ....: JordanSpinEJA)
2251 When no bilinear form is specified, the identity matrix is used,
2252 and the resulting algebra is the Jordan spin algebra::
2254 sage: B = matrix.identity(AA,3)
2255 sage: J0 = BilinearFormEJA(B)
2256 sage: J1 = JordanSpinEJA(3)
2257 sage: J0.multiplication_table() == J0.multiplication_table()
2260 An error is raised if the matrix `B` does not correspond to a
2261 positive-definite bilinear form::
2263 sage: B = matrix.random(QQ,2,3)
2264 sage: J = BilinearFormEJA(B)
2265 Traceback (most recent call last):
2267 ValueError: bilinear form is not positive-definite
2268 sage: B = matrix.zero(QQ,3)
2269 sage: J = BilinearFormEJA(B)
2270 Traceback (most recent call last):
2272 ValueError: bilinear form is not positive-definite
2276 We can create a zero-dimensional algebra::
2278 sage: B = matrix.identity(AA,0)
2279 sage: J = BilinearFormEJA(B)
2283 We can check the multiplication condition given in the Jordan, von
2284 Neumann, and Wigner paper (and also discussed on my "On the
2285 symmetry..." paper). Note that this relies heavily on the standard
2286 choice of basis, as does anything utilizing the bilinear form
2287 matrix. We opt not to orthonormalize the basis, because if we
2288 did, we would have to normalize the `s_{i}` in a similar manner::
2290 sage: set_random_seed()
2291 sage: n = ZZ.random_element(5)
2292 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2293 sage: B11 = matrix.identity(QQ,1)
2294 sage: B22 = M.transpose()*M
2295 sage: B = block_matrix(2,2,[ [B11,0 ],
2297 sage: J = BilinearFormEJA(B, orthonormalize=False)
2298 sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
2299 sage: V = J.vector_space()
2300 sage: sis = [ J( V([0] + (M.inverse()*ei).list()).column() )
2301 ....: for ei in eis ]
2302 sage: actual = [ sis[i]*sis[j]
2303 ....: for i in range(n-1)
2304 ....: for j in range(n-1) ]
2305 sage: expected = [ J.one() if i == j else J.zero()
2306 ....: for i in range(n-1)
2307 ....: for j in range(n-1) ]
2308 sage: actual == expected
2312 def __init__(self
, B
, **kwargs
):
2313 # The matrix "B" is supplied by the user in most cases,
2314 # so it makes sense to check whether or not its positive-
2315 # definite unless we are specifically asked not to...
2316 if ("check_axioms" not in kwargs
) or kwargs
["check_axioms"]:
2317 if not B
.is_positive_definite():
2318 raise ValueError("bilinear form is not positive-definite")
2320 # However, all of the other data for this EJA is computed
2321 # by us in manner that guarantees the axioms are
2322 # satisfied. So, again, unless we are specifically asked to
2323 # verify things, we'll skip the rest of the checks.
2324 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2326 def inner_product(x
,y
):
2327 return (y
.T
*B
*x
)[0,0]
2329 def jordan_product(x
,y
):
2335 z0
= inner_product(y
,x
)
2336 zbar
= y0
*xbar
+ x0
*ybar
2337 return P([z0
] + zbar
.list())
2340 column_basis
= tuple( b
.column() for b
in FreeModule(ZZ
, n
).basis() )
2341 super(BilinearFormEJA
, self
).__init
__(column_basis
,
2346 # The rank of this algebra is two, unless we're in a
2347 # one-dimensional ambient space (because the rank is bounded
2348 # by the ambient dimension).
2349 self
.rank
.set_cache(min(n
,2))
2352 self
.one
.set_cache( self
.zero() )
2354 self
.one
.set_cache( self
.monomial(0) )
2357 def _max_random_instance_size():
2359 The maximum dimension of a random BilinearFormEJA.
2364 def random_instance(cls
, **kwargs
):
2366 Return a random instance of this algebra.
2368 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2370 B
= matrix
.identity(ZZ
, n
)
2371 return cls(B
, **kwargs
)
2373 B11
= matrix
.identity(ZZ
, 1)
2374 M
= matrix
.random(ZZ
, n
-1)
2375 I
= matrix
.identity(ZZ
, n
-1)
2377 while alpha
.is_zero():
2378 alpha
= ZZ
.random_element().abs()
2379 B22
= M
.transpose()*M
+ alpha
*I
2381 from sage
.matrix
.special
import block_matrix
2382 B
= block_matrix(2,2, [ [B11
, ZZ(0) ],
2385 return cls(B
, **kwargs
)
2388 class JordanSpinEJA(BilinearFormEJA
):
2390 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2391 with the usual inner product and jordan product ``x*y =
2392 (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2397 sage: from mjo.eja.eja_algebra import JordanSpinEJA
2401 This multiplication table can be verified by hand::
2403 sage: J = JordanSpinEJA(4)
2404 sage: e0,e1,e2,e3 = J.gens()
2420 We can change the generator prefix::
2422 sage: JordanSpinEJA(2, prefix='B').gens()
2427 Ensure that we have the usual inner product on `R^n`::
2429 sage: set_random_seed()
2430 sage: J = JordanSpinEJA.random_instance()
2431 sage: x,y = J.random_elements(2)
2432 sage: actual = x.inner_product(y)
2433 sage: expected = x.to_vector().inner_product(y.to_vector())
2434 sage: actual == expected
2438 def __init__(self
, n
, **kwargs
):
2439 # This is a special case of the BilinearFormEJA with the
2440 # identity matrix as its bilinear form.
2441 B
= matrix
.identity(ZZ
, n
)
2443 # Don't orthonormalize because our basis is already
2444 # orthonormal with respect to our inner-product.
2445 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2447 # But also don't pass check_field=False here, because the user
2448 # can pass in a field!
2449 super(JordanSpinEJA
, self
).__init
__(B
, **kwargs
)
2452 def _max_random_instance_size():
2454 The maximum dimension of a random JordanSpinEJA.
2459 def random_instance(cls
, **kwargs
):
2461 Return a random instance of this type of algebra.
2463 Needed here to override the implementation for ``BilinearFormEJA``.
2465 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2466 return cls(n
, **kwargs
)
2469 class TrivialEJA(ConcreteEJA
):
2471 The trivial Euclidean Jordan algebra consisting of only a zero element.
2475 sage: from mjo.eja.eja_algebra import TrivialEJA
2479 sage: J = TrivialEJA()
2486 sage: 7*J.one()*12*J.one()
2488 sage: J.one().inner_product(J.one())
2490 sage: J.one().norm()
2492 sage: J.one().subalgebra_generated_by()
2493 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2498 def __init__(self
, **kwargs
):
2499 jordan_product
= lambda x
,y
: x
2500 inner_product
= lambda x
,y
: 0
2503 # New defaults for keyword arguments
2504 if "orthonormalize" not in kwargs
: kwargs
["orthonormalize"] = False
2505 if "check_axioms" not in kwargs
: kwargs
["check_axioms"] = False
2507 super(TrivialEJA
, self
).__init
__(basis
,
2511 # The rank is zero using my definition, namely the dimension of the
2512 # largest subalgebra generated by any element.
2513 self
.rank
.set_cache(0)
2514 self
.one
.set_cache( self
.zero() )
2517 def random_instance(cls
, **kwargs
):
2518 # We don't take a "size" argument so the superclass method is
2519 # inappropriate for us.
2520 return cls(**kwargs
)
2522 # class DirectSumEJA(ConcreteEJA):
2524 # The external (orthogonal) direct sum of two other Euclidean Jordan
2525 # algebras. Essentially the Cartesian product of its two factors.
2526 # Every Euclidean Jordan algebra decomposes into an orthogonal
2527 # direct sum of simple Euclidean Jordan algebras, so no generality
2528 # is lost by providing only this construction.
2532 # sage: from mjo.eja.eja_algebra import (random_eja,
2533 # ....: HadamardEJA,
2534 # ....: RealSymmetricEJA,
2535 # ....: DirectSumEJA)
2539 # sage: J1 = HadamardEJA(2)
2540 # sage: J2 = RealSymmetricEJA(3)
2541 # sage: J = DirectSumEJA(J1,J2)
2542 # sage: J.dimension()
2549 # The external direct sum construction is only valid when the two factors
2550 # have the same base ring; an error is raised otherwise::
2552 # sage: set_random_seed()
2553 # sage: J1 = random_eja(field=AA)
2554 # sage: J2 = random_eja(field=QQ,orthonormalize=False)
2555 # sage: J = DirectSumEJA(J1,J2)
2556 # Traceback (most recent call last):
2558 # ValueError: algebras must share the same base field
2561 # def __init__(self, J1, J2, **kwargs):
2562 # if J1.base_ring() != J2.base_ring():
2563 # raise ValueError("algebras must share the same base field")
2564 # field = J1.base_ring()
2566 # self._factors = (J1, J2)
2567 # n1 = J1.dimension()
2568 # n2 = J2.dimension()
2570 # V = VectorSpace(field, n)
2571 # mult_table = [ [ V.zero() for j in range(i+1) ]
2572 # for i in range(n) ]
2573 # for i in range(n1):
2574 # for j in range(i+1):
2575 # p = (J1.monomial(i)*J1.monomial(j)).to_vector()
2576 # mult_table[i][j] = V(p.list() + [field.zero()]*n2)
2578 # for i in range(n2):
2579 # for j in range(i+1):
2580 # p = (J2.monomial(i)*J2.monomial(j)).to_vector()
2581 # mult_table[n1+i][n1+j] = V([field.zero()]*n1 + p.list())
2583 # # TODO: build the IP table here from the two constituent IP
2584 # # matrices (it'll be block diagonal, I think).
2585 # ip_table = [ [ field.zero() for j in range(i+1) ]
2586 # for i in range(n) ]
2587 # super(DirectSumEJA, self).__init__(field,
2590 # check_axioms=False,
2592 # self.rank.set_cache(J1.rank() + J2.rank())
2595 # def factors(self):
2597 # Return the pair of this algebra's factors.
2601 # sage: from mjo.eja.eja_algebra import (HadamardEJA,
2602 # ....: JordanSpinEJA,
2603 # ....: DirectSumEJA)
2607 # sage: J1 = HadamardEJA(2, field=QQ)
2608 # sage: J2 = JordanSpinEJA(3, field=QQ)
2609 # sage: J = DirectSumEJA(J1,J2)
2611 # (Euclidean Jordan algebra of dimension 2 over Rational Field,
2612 # Euclidean Jordan algebra of dimension 3 over Rational Field)
2615 # return self._factors
2617 # def projections(self):
2619 # Return a pair of projections onto this algebra's factors.
2623 # sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
2624 # ....: ComplexHermitianEJA,
2625 # ....: DirectSumEJA)
2629 # sage: J1 = JordanSpinEJA(2)
2630 # sage: J2 = ComplexHermitianEJA(2)
2631 # sage: J = DirectSumEJA(J1,J2)
2632 # sage: (pi_left, pi_right) = J.projections()
2633 # sage: J.one().to_vector()
2634 # (1, 0, 1, 0, 0, 1)
2635 # sage: pi_left(J.one()).to_vector()
2637 # sage: pi_right(J.one()).to_vector()
2641 # (J1,J2) = self.factors()
2642 # m = J1.dimension()
2643 # n = J2.dimension()
2644 # V_basis = self.vector_space().basis()
2645 # # Need to specify the dimensions explicitly so that we don't
2646 # # wind up with a zero-by-zero matrix when we want e.g. a
2647 # # zero-by-two matrix (important for composing things).
2648 # P1 = matrix(self.base_ring(), m, m+n, V_basis[:m])
2649 # P2 = matrix(self.base_ring(), n, m+n, V_basis[m:])
2650 # pi_left = FiniteDimensionalEJAOperator(self,J1,P1)
2651 # pi_right = FiniteDimensionalEJAOperator(self,J2,P2)
2652 # return (pi_left, pi_right)
2654 # def inclusions(self):
2656 # Return the pair of inclusion maps from our factors into us.
2660 # sage: from mjo.eja.eja_algebra import (random_eja,
2661 # ....: JordanSpinEJA,
2662 # ....: RealSymmetricEJA,
2663 # ....: DirectSumEJA)
2667 # sage: J1 = JordanSpinEJA(3)
2668 # sage: J2 = RealSymmetricEJA(2)
2669 # sage: J = DirectSumEJA(J1,J2)
2670 # sage: (iota_left, iota_right) = J.inclusions()
2671 # sage: iota_left(J1.zero()) == J.zero()
2673 # sage: iota_right(J2.zero()) == J.zero()
2675 # sage: J1.one().to_vector()
2677 # sage: iota_left(J1.one()).to_vector()
2678 # (1, 0, 0, 0, 0, 0)
2679 # sage: J2.one().to_vector()
2681 # sage: iota_right(J2.one()).to_vector()
2682 # (0, 0, 0, 1, 0, 1)
2683 # sage: J.one().to_vector()
2684 # (1, 0, 0, 1, 0, 1)
2688 # Composing a projection with the corresponding inclusion should
2689 # produce the identity map, and mismatching them should produce
2692 # sage: set_random_seed()
2693 # sage: J1 = random_eja()
2694 # sage: J2 = random_eja()
2695 # sage: J = DirectSumEJA(J1,J2)
2696 # sage: (iota_left, iota_right) = J.inclusions()
2697 # sage: (pi_left, pi_right) = J.projections()
2698 # sage: pi_left*iota_left == J1.one().operator()
2700 # sage: pi_right*iota_right == J2.one().operator()
2702 # sage: (pi_left*iota_right).is_zero()
2704 # sage: (pi_right*iota_left).is_zero()
2708 # (J1,J2) = self.factors()
2709 # m = J1.dimension()
2710 # n = J2.dimension()
2711 # V_basis = self.vector_space().basis()
2712 # # Need to specify the dimensions explicitly so that we don't
2713 # # wind up with a zero-by-zero matrix when we want e.g. a
2714 # # two-by-zero matrix (important for composing things).
2715 # I1 = matrix.column(self.base_ring(), m, m+n, V_basis[:m])
2716 # I2 = matrix.column(self.base_ring(), n, m+n, V_basis[m:])
2717 # iota_left = FiniteDimensionalEJAOperator(J1,self,I1)
2718 # iota_right = FiniteDimensionalEJAOperator(J2,self,I2)
2719 # return (iota_left, iota_right)
2721 # def inner_product(self, x, y):
2723 # The standard Cartesian inner-product.
2725 # We project ``x`` and ``y`` onto our factors, and add up the
2726 # inner-products from the subalgebras.
2731 # sage: from mjo.eja.eja_algebra import (HadamardEJA,
2732 # ....: QuaternionHermitianEJA,
2733 # ....: DirectSumEJA)
2737 # sage: J1 = HadamardEJA(3,field=QQ)
2738 # sage: J2 = QuaternionHermitianEJA(2,field=QQ,orthonormalize=False)
2739 # sage: J = DirectSumEJA(J1,J2)
2740 # sage: x1 = J1.one()
2742 # sage: y1 = J2.one()
2744 # sage: x1.inner_product(x2)
2746 # sage: y1.inner_product(y2)
2748 # sage: J.one().inner_product(J.one())
2752 # (pi_left, pi_right) = self.projections()
2758 # return (x1.inner_product(y1) + x2.inner_product(y2))
2762 random_eja
= ConcreteEJA
.random_instance