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.
8 from itertools
import repeat
10 from sage
.algebras
.quatalg
.quaternion_algebra
import QuaternionAlgebra
11 from sage
.categories
.magmatic_algebras
import MagmaticAlgebras
12 from sage
.combinat
.free_module
import CombinatorialFreeModule
13 from sage
.matrix
.constructor
import matrix
14 from sage
.matrix
.matrix_space
import MatrixSpace
15 from sage
.misc
.cachefunc
import cached_method
16 from sage
.misc
.lazy_import
import lazy_import
17 from sage
.misc
.prandom
import choice
18 from sage
.misc
.table
import table
19 from sage
.modules
.free_module
import FreeModule
, VectorSpace
20 from sage
.rings
.all
import (ZZ
, QQ
, RR
, RLF
, CLF
,
23 from mjo
.eja
.eja_element
import FiniteDimensionalEuclideanJordanAlgebraElement
24 lazy_import('mjo.eja.eja_subalgebra',
25 'FiniteDimensionalEuclideanJordanSubalgebra')
26 from mjo
.eja
.eja_utils
import _mat2vec
28 class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule
):
29 # This is an ugly hack needed to prevent the category framework
30 # from implementing a coercion from our base ring (e.g. the
31 # rationals) into the algebra. First of all -- such a coercion is
32 # nonsense to begin with. But more importantly, it tries to do so
33 # in the category of rings, and since our algebras aren't
34 # associative they generally won't be rings.
35 _no_generic_basering_coercion
= True
48 sage: from mjo.eja.eja_algebra import (JordanSpinEJA, random_eja)
52 By definition, Jordan multiplication commutes::
54 sage: set_random_seed()
55 sage: J = random_eja()
56 sage: x,y = J.random_elements(2)
62 The ``field`` we're given must be real::
64 sage: JordanSpinEJA(2,QQbar)
65 Traceback (most recent call last):
67 ValueError: field is not real
71 if not field
.is_subring(RR
):
72 # Note: this does return true for the real algebraic
73 # field, and any quadratic field where we've specified
75 raise ValueError('field is not real')
78 self
._natural
_basis
= natural_basis
81 category
= MagmaticAlgebras(field
).FiniteDimensional()
82 category
= category
.WithBasis().Unital()
84 fda
= super(FiniteDimensionalEuclideanJordanAlgebra
, self
)
86 range(len(mult_table
)),
89 self
.print_options(bracket
='')
91 # The multiplication table we're given is necessarily in terms
92 # of vectors, because we don't have an algebra yet for
93 # anything to be an element of. However, it's faster in the
94 # long run to have the multiplication table be in terms of
95 # algebra elements. We do this after calling the superclass
96 # constructor so that from_vector() knows what to do.
97 self
._multiplication
_table
= [
98 list(map(lambda x
: self
.from_vector(x
), ls
))
103 def _element_constructor_(self
, elt
):
105 Construct an element of this algebra from its natural
108 This gets called only after the parent element _call_ method
109 fails to find a coercion for the argument.
113 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
114 ....: RealCartesianProductEJA,
115 ....: RealSymmetricEJA)
119 The identity in `S^n` is converted to the identity in the EJA::
121 sage: J = RealSymmetricEJA(3)
122 sage: I = matrix.identity(QQ,3)
123 sage: J(I) == J.one()
126 This skew-symmetric matrix can't be represented in the EJA::
128 sage: J = RealSymmetricEJA(3)
129 sage: A = matrix(QQ,3, lambda i,j: i-j)
131 Traceback (most recent call last):
133 ArithmeticError: vector is not in free module
137 Ensure that we can convert any element of the two non-matrix
138 simple algebras (whose natural representations are their usual
139 vector representations) back and forth faithfully::
141 sage: set_random_seed()
142 sage: J = RealCartesianProductEJA.random_instance()
143 sage: x = J.random_element()
144 sage: J(x.to_vector().column()) == x
146 sage: J = JordanSpinEJA.random_instance()
147 sage: x = J.random_element()
148 sage: J(x.to_vector().column()) == x
153 # The superclass implementation of random_element()
154 # needs to be able to coerce "0" into the algebra.
157 natural_basis
= self
.natural_basis()
158 basis_space
= natural_basis
[0].matrix_space()
159 if elt
not in basis_space
:
160 raise ValueError("not a naturally-represented algebra element")
162 # Thanks for nothing! Matrix spaces aren't vector spaces in
163 # Sage, so we have to figure out its natural-basis coordinates
164 # ourselves. We use the basis space's ring instead of the
165 # element's ring because the basis space might be an algebraic
166 # closure whereas the base ring of the 3-by-3 identity matrix
167 # could be QQ instead of QQbar.
168 V
= VectorSpace(basis_space
.base_ring(), elt
.nrows()*elt
.ncols())
169 W
= V
.span_of_basis( _mat2vec(s
) for s
in natural_basis
)
170 coords
= W
.coordinate_vector(_mat2vec(elt
))
171 return self
.from_vector(coords
)
176 Return a string representation of ``self``.
180 sage: from mjo.eja.eja_algebra import JordanSpinEJA
184 Ensure that it says what we think it says::
186 sage: JordanSpinEJA(2, field=QQ)
187 Euclidean Jordan algebra of dimension 2 over Rational Field
188 sage: JordanSpinEJA(3, field=RDF)
189 Euclidean Jordan algebra of dimension 3 over Real Double Field
192 fmt
= "Euclidean Jordan algebra of dimension {} over {}"
193 return fmt
.format(self
.dimension(), self
.base_ring())
195 def product_on_basis(self
, i
, j
):
196 return self
._multiplication
_table
[i
][j
]
198 def _a_regular_element(self
):
200 Guess a regular element. Needed to compute the basis for our
201 characteristic polynomial coefficients.
205 sage: from mjo.eja.eja_algebra import random_eja
209 Ensure that this hacky method succeeds for every algebra that we
210 know how to construct::
212 sage: set_random_seed()
213 sage: J = random_eja()
214 sage: J._a_regular_element().is_regular()
219 z
= self
.sum( (i
+1)*gs
[i
] for i
in range(len(gs
)) )
220 if not z
.is_regular():
221 raise ValueError("don't know a regular element")
226 def _charpoly_basis_space(self
):
228 Return the vector space spanned by the basis used in our
229 characteristic polynomial coefficients. This is used not only to
230 compute those coefficients, but also any time we need to
231 evaluate the coefficients (like when we compute the trace or
234 z
= self
._a
_regular
_element
()
235 # Don't use the parent vector space directly here in case this
236 # happens to be a subalgebra. In that case, we would be e.g.
237 # two-dimensional but span_of_basis() would expect three
239 V
= VectorSpace(self
.base_ring(), self
.vector_space().dimension())
240 basis
= [ (z
**k
).to_vector() for k
in range(self
.rank()) ]
241 V1
= V
.span_of_basis( basis
)
242 b
= (V1
.basis() + V1
.complement().basis())
243 return V
.span_of_basis(b
)
248 def _charpoly_coeff(self
, i
):
250 Return the coefficient polynomial "a_{i}" of this algebra's
251 general characteristic polynomial.
253 Having this be a separate cached method lets us compute and
254 store the trace/determinant (a_{r-1} and a_{0} respectively)
255 separate from the entire characteristic polynomial.
257 (A_of_x
, x
, xr
, detA
) = self
._charpoly
_matrix
_system
()
258 R
= A_of_x
.base_ring()
263 # Guaranteed by theory
266 # Danger: the in-place modification is done for performance
267 # reasons (reconstructing a matrix with huge polynomial
268 # entries is slow), but I don't know how cached_method works,
269 # so it's highly possible that we're modifying some global
270 # list variable by reference, here. In other words, you
271 # probably shouldn't call this method twice on the same
272 # algebra, at the same time, in two threads
273 Ai_orig
= A_of_x
.column(i
)
274 A_of_x
.set_column(i
,xr
)
275 numerator
= A_of_x
.det()
276 A_of_x
.set_column(i
,Ai_orig
)
278 # We're relying on the theory here to ensure that each a_i is
279 # indeed back in R, and the added negative signs are to make
280 # the whole charpoly expression sum to zero.
281 return R(-numerator
/detA
)
285 def _charpoly_matrix_system(self
):
287 Compute the matrix whose entries A_ij are polynomials in
288 X1,...,XN, the vector ``x`` of variables X1,...,XN, the vector
289 corresponding to `x^r` and the determinent of the matrix A =
290 [A_ij]. In other words, all of the fixed (cachable) data needed
291 to compute the coefficients of the characteristic polynomial.
296 # Turn my vector space into a module so that "vectors" can
297 # have multivatiate polynomial entries.
298 names
= tuple('X' + str(i
) for i
in range(1,n
+1))
299 R
= PolynomialRing(self
.base_ring(), names
)
301 # Using change_ring() on the parent's vector space doesn't work
302 # here because, in a subalgebra, that vector space has a basis
303 # and change_ring() tries to bring the basis along with it. And
304 # that doesn't work unless the new ring is a PID, which it usually
308 # Now let x = (X1,X2,...,Xn) be the vector whose entries are
312 # And figure out the "left multiplication by x" matrix in
315 monomial_matrices
= [ self
.monomial(i
).operator().matrix()
316 for i
in range(n
) ] # don't recompute these!
318 ek
= self
.monomial(k
).to_vector()
320 sum( x
[i
]*(monomial_matrices
[i
]*ek
)
321 for i
in range(n
) ) )
322 Lx
= matrix
.column(R
, lmbx_cols
)
324 # Now we can compute powers of x "symbolically"
325 x_powers
= [self
.one().to_vector(), x
]
326 for d
in range(2, r
+1):
327 x_powers
.append( Lx
*(x_powers
[-1]) )
329 idmat
= matrix
.identity(R
, n
)
331 W
= self
._charpoly
_basis
_space
()
332 W
= W
.change_ring(R
.fraction_field())
334 # Starting with the standard coordinates x = (X1,X2,...,Xn)
335 # and then converting the entries to W-coordinates allows us
336 # to pass in the standard coordinates to the charpoly and get
337 # back the right answer. Specifically, with x = (X1,X2,...,Xn),
340 # W.coordinates(x^2) eval'd at (standard z-coords)
344 # W-coords of (standard coords of x^2 eval'd at std-coords of z)
346 # We want the middle equivalent thing in our matrix, but use
347 # the first equivalent thing instead so that we can pass in
348 # standard coordinates.
349 x_powers
= [ W
.coordinate_vector(xp
) for xp
in x_powers
]
350 l2
= [idmat
.column(k
-1) for k
in range(r
+1, n
+1)]
351 A_of_x
= matrix
.column(R
, n
, (x_powers
[:r
] + l2
))
352 return (A_of_x
, x
, x_powers
[r
], A_of_x
.det())
356 def characteristic_polynomial(self
):
358 Return a characteristic polynomial that works for all elements
361 The resulting polynomial has `n+1` variables, where `n` is the
362 dimension of this algebra. The first `n` variables correspond to
363 the coordinates of an algebra element: when evaluated at the
364 coordinates of an algebra element with respect to a certain
365 basis, the result is a univariate polynomial (in the one
366 remaining variable ``t``), namely the characteristic polynomial
371 sage: from mjo.eja.eja_algebra import JordanSpinEJA, TrivialEJA
375 The characteristic polynomial in the spin algebra is given in
376 Alizadeh, Example 11.11::
378 sage: J = JordanSpinEJA(3)
379 sage: p = J.characteristic_polynomial(); p
380 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
381 sage: xvec = J.one().to_vector()
385 By definition, the characteristic polynomial is a monic
386 degree-zero polynomial in a rank-zero algebra. Note that
387 Cayley-Hamilton is indeed satisfied since the polynomial
388 ``1`` evaluates to the identity element of the algebra on
391 sage: J = TrivialEJA()
392 sage: J.characteristic_polynomial()
399 # The list of coefficient polynomials a_0, a_1, a_2, ..., a_n.
400 a
= [ self
._charpoly
_coeff
(i
) for i
in range(r
+1) ]
402 # We go to a bit of trouble here to reorder the
403 # indeterminates, so that it's easier to evaluate the
404 # characteristic polynomial at x's coordinates and get back
405 # something in terms of t, which is what we want.
407 S
= PolynomialRing(self
.base_ring(),'t')
409 S
= PolynomialRing(S
, R
.variable_names())
412 return sum( a
[k
]*(t
**k
) for k
in range(len(a
)) )
415 def inner_product(self
, x
, y
):
417 The inner product associated with this Euclidean Jordan algebra.
419 Defaults to the trace inner product, but can be overridden by
420 subclasses if they are sure that the necessary properties are
425 sage: from mjo.eja.eja_algebra import random_eja
429 Our inner product is "associative," which means the following for
430 a symmetric bilinear form::
432 sage: set_random_seed()
433 sage: J = random_eja()
434 sage: x,y,z = J.random_elements(3)
435 sage: (x*y).inner_product(z) == y.inner_product(x*z)
439 X
= x
.natural_representation()
440 Y
= y
.natural_representation()
441 return self
.natural_inner_product(X
,Y
)
444 def is_trivial(self
):
446 Return whether or not this algebra is trivial.
448 A trivial algebra contains only the zero element.
452 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
457 sage: J = ComplexHermitianEJA(3)
463 sage: J = TrivialEJA()
468 return self
.dimension() == 0
471 def multiplication_table(self
):
473 Return a visual representation of this algebra's multiplication
474 table (on basis elements).
478 sage: from mjo.eja.eja_algebra import JordanSpinEJA
482 sage: J = JordanSpinEJA(4)
483 sage: J.multiplication_table()
484 +----++----+----+----+----+
485 | * || e0 | e1 | e2 | e3 |
486 +====++====+====+====+====+
487 | e0 || e0 | e1 | e2 | e3 |
488 +----++----+----+----+----+
489 | e1 || e1 | e0 | 0 | 0 |
490 +----++----+----+----+----+
491 | e2 || e2 | 0 | e0 | 0 |
492 +----++----+----+----+----+
493 | e3 || e3 | 0 | 0 | e0 |
494 +----++----+----+----+----+
497 M
= list(self
._multiplication
_table
) # copy
498 for i
in range(len(M
)):
499 # M had better be "square"
500 M
[i
] = [self
.monomial(i
)] + M
[i
]
501 M
= [["*"] + list(self
.gens())] + M
502 return table(M
, header_row
=True, header_column
=True, frame
=True)
505 def natural_basis(self
):
507 Return a more-natural representation of this algebra's basis.
509 Every finite-dimensional Euclidean Jordan Algebra is a direct
510 sum of five simple algebras, four of which comprise Hermitian
511 matrices. This method returns the original "natural" basis
512 for our underlying vector space. (Typically, the natural basis
513 is used to construct the multiplication table in the first place.)
515 Note that this will always return a matrix. The standard basis
516 in `R^n` will be returned as `n`-by-`1` column matrices.
520 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
521 ....: RealSymmetricEJA)
525 sage: J = RealSymmetricEJA(2)
527 Finite family {0: e0, 1: e1, 2: e2}
528 sage: J.natural_basis()
530 [1 0] [ 0 1/2*sqrt2] [0 0]
531 [0 0], [1/2*sqrt2 0], [0 1]
536 sage: J = JordanSpinEJA(2)
538 Finite family {0: e0, 1: e1}
539 sage: J.natural_basis()
546 if self
._natural
_basis
is None:
547 M
= self
.natural_basis_space()
548 return tuple( M(b
.to_vector()) for b
in self
.basis() )
550 return self
._natural
_basis
553 def natural_basis_space(self
):
555 Return the matrix space in which this algebra's natural basis
558 if self
._natural
_basis
is None or len(self
._natural
_basis
) == 0:
559 return MatrixSpace(self
.base_ring(), self
.dimension(), 1)
561 return self
._natural
_basis
[0].matrix_space()
565 def natural_inner_product(X
,Y
):
567 Compute the inner product of two naturally-represented elements.
569 For example in the real symmetric matrix EJA, this will compute
570 the trace inner-product of two n-by-n symmetric matrices. The
571 default should work for the real cartesian product EJA, the
572 Jordan spin EJA, and the real symmetric matrices. The others
573 will have to be overridden.
575 return (X
.conjugate_transpose()*Y
).trace()
581 Return the unit element of this algebra.
585 sage: from mjo.eja.eja_algebra import (RealCartesianProductEJA,
590 sage: J = RealCartesianProductEJA(5)
592 e0 + e1 + e2 + e3 + e4
596 The identity element acts like the identity::
598 sage: set_random_seed()
599 sage: J = random_eja()
600 sage: x = J.random_element()
601 sage: J.one()*x == x and x*J.one() == x
604 The matrix of the unit element's operator is the identity::
606 sage: set_random_seed()
607 sage: J = random_eja()
608 sage: actual = J.one().operator().matrix()
609 sage: expected = matrix.identity(J.base_ring(), J.dimension())
610 sage: actual == expected
614 # We can brute-force compute the matrices of the operators
615 # that correspond to the basis elements of this algebra.
616 # If some linear combination of those basis elements is the
617 # algebra identity, then the same linear combination of
618 # their matrices has to be the identity matrix.
620 # Of course, matrices aren't vectors in sage, so we have to
621 # appeal to the "long vectors" isometry.
622 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
624 # Now we use basis linear algebra to find the coefficients,
625 # of the matrices-as-vectors-linear-combination, which should
626 # work for the original algebra basis too.
627 A
= matrix
.column(self
.base_ring(), oper_vecs
)
629 # We used the isometry on the left-hand side already, but we
630 # still need to do it for the right-hand side. Recall that we
631 # wanted something that summed to the identity matrix.
632 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
634 # Now if there's an identity element in the algebra, this should work.
635 coeffs
= A
.solve_right(b
)
636 return self
.linear_combination(zip(self
.gens(), coeffs
))
639 def peirce_decomposition(self
, c
):
641 The Peirce decomposition of this algebra relative to the
644 In the future, this can be extended to a complete system of
645 orthogonal idempotents.
649 - ``c`` -- an idempotent of this algebra.
653 A triple (J0, J5, J1) containing two subalgebras and one subspace
656 - ``J0`` -- the algebra on the eigenspace of ``c.operator()``
657 corresponding to the eigenvalue zero.
659 - ``J5`` -- the eigenspace (NOT a subalgebra) of ``c.operator()``
660 corresponding to the eigenvalue one-half.
662 - ``J1`` -- the algebra on the eigenspace of ``c.operator()``
663 corresponding to the eigenvalue one.
665 These are the only possible eigenspaces for that operator, and this
666 algebra is a direct sum of them. The spaces ``J0`` and ``J1`` are
667 orthogonal, and are subalgebras of this algebra with the appropriate
672 sage: from mjo.eja.eja_algebra import random_eja, RealSymmetricEJA
676 The canonical example comes from the symmetric matrices, which
677 decompose into diagonal and off-diagonal parts::
679 sage: J = RealSymmetricEJA(3)
680 sage: C = matrix(QQ, [ [1,0,0],
684 sage: J0,J5,J1 = J.peirce_decomposition(c)
686 Euclidean Jordan algebra of dimension 1...
688 Vector space of degree 6 and dimension 2...
690 Euclidean Jordan algebra of dimension 3...
694 Every algebra decomposes trivially with respect to its identity
697 sage: set_random_seed()
698 sage: J = random_eja()
699 sage: J0,J5,J1 = J.peirce_decomposition(J.one())
700 sage: J0.dimension() == 0 and J5.dimension() == 0
702 sage: J1.superalgebra() == J and J1.dimension() == J.dimension()
705 The identity elements in the two subalgebras are the
706 projections onto their respective subspaces of the
707 superalgebra's identity element::
709 sage: set_random_seed()
710 sage: J = random_eja()
711 sage: x = J.random_element()
712 sage: if not J.is_trivial():
713 ....: while x.is_nilpotent():
714 ....: x = J.random_element()
715 sage: c = x.subalgebra_idempotent()
716 sage: J0,J5,J1 = J.peirce_decomposition(c)
717 sage: J1(c) == J1.one()
719 sage: J0(J.one() - c) == J0.one()
723 if not c
.is_idempotent():
724 raise ValueError("element is not idempotent: %s" % c
)
726 # Default these to what they should be if they turn out to be
727 # trivial, because eigenspaces_left() won't return eigenvalues
728 # corresponding to trivial spaces (e.g. it returns only the
729 # eigenspace corresponding to lambda=1 if you take the
730 # decomposition relative to the identity element).
731 trivial
= FiniteDimensionalEuclideanJordanSubalgebra(self
, ())
732 J0
= trivial
# eigenvalue zero
733 J5
= VectorSpace(self
.base_ring(), 0) # eigenvalue one-half
734 J1
= trivial
# eigenvalue one
736 for (eigval
, eigspace
) in c
.operator().matrix().left_eigenspaces():
737 if eigval
== ~
(self
.base_ring()(2)):
740 gens
= tuple( self
.from_vector(b
) for b
in eigspace
.basis() )
741 subalg
= FiniteDimensionalEuclideanJordanSubalgebra(self
, gens
)
747 raise ValueError("unexpected eigenvalue: %s" % eigval
)
752 def random_elements(self
, count
):
754 Return ``count`` random elements as a tuple.
758 sage: from mjo.eja.eja_algebra import JordanSpinEJA
762 sage: J = JordanSpinEJA(3)
763 sage: x,y,z = J.random_elements(3)
764 sage: all( [ x in J, y in J, z in J ])
766 sage: len( J.random_elements(10) ) == 10
770 return tuple( self
.random_element() for idx
in range(count
) )
775 Return the rank of this EJA.
779 The author knows of no algorithm to compute the rank of an EJA
780 where only the multiplication table is known. In lieu of one, we
781 require the rank to be specified when the algebra is created,
782 and simply pass along that number here.
786 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
787 ....: RealSymmetricEJA,
788 ....: ComplexHermitianEJA,
789 ....: QuaternionHermitianEJA,
794 The rank of the Jordan spin algebra is always two::
796 sage: JordanSpinEJA(2).rank()
798 sage: JordanSpinEJA(3).rank()
800 sage: JordanSpinEJA(4).rank()
803 The rank of the `n`-by-`n` Hermitian real, complex, or
804 quaternion matrices is `n`::
806 sage: RealSymmetricEJA(4).rank()
808 sage: ComplexHermitianEJA(3).rank()
810 sage: QuaternionHermitianEJA(2).rank()
815 Ensure that every EJA that we know how to construct has a
816 positive integer rank, unless the algebra is trivial in
817 which case its rank will be zero::
819 sage: set_random_seed()
820 sage: J = random_eja()
824 sage: r > 0 or (r == 0 and J.is_trivial())
831 def vector_space(self
):
833 Return the vector space that underlies this algebra.
837 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
841 sage: J = RealSymmetricEJA(2)
842 sage: J.vector_space()
843 Vector space of dimension 3 over...
846 return self
.zero().to_vector().parent().ambient_vector_space()
849 Element
= FiniteDimensionalEuclideanJordanAlgebraElement
852 class KnownRankEJA(object):
854 A class for algebras that we actually know we can construct. The
855 main issue is that, for most of our methods to make sense, we need
856 to know the rank of our algebra. Thus we can't simply generate a
857 "random" algebra, or even check that a given basis and product
858 satisfy the axioms; because even if everything looks OK, we wouldn't
859 know the rank we need to actuallty build the thing.
861 Not really a subclass of FDEJA because doing that causes method
862 resolution errors, e.g.
864 TypeError: Error when calling the metaclass bases
865 Cannot create a consistent method resolution
866 order (MRO) for bases FiniteDimensionalEuclideanJordanAlgebra,
871 def _max_test_case_size():
873 Return an integer "size" that is an upper bound on the size of
874 this algebra when it is used in a random test
875 case. Unfortunately, the term "size" is quite vague -- when
876 dealing with `R^n` under either the Hadamard or Jordan spin
877 product, the "size" refers to the dimension `n`. When dealing
878 with a matrix algebra (real symmetric or complex/quaternion
879 Hermitian), it refers to the size of the matrix, which is
880 far less than the dimension of the underlying vector space.
882 We default to five in this class, which is safe in `R^n`. The
883 matrix algebra subclasses (or any class where the "size" is
884 interpreted to be far less than the dimension) should override
885 with a smaller number.
890 def random_instance(cls
, field
=QQ
, **kwargs
):
892 Return a random instance of this type of algebra.
894 Beware, this will crash for "most instances" because the
895 constructor below looks wrong.
897 if cls
is TrivialEJA
:
898 # The TrivialEJA class doesn't take an "n" argument because
902 n
= ZZ
.random_element(cls
._max
_test
_case
_size
()) + 1
903 return cls(n
, field
, **kwargs
)
906 class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra
,
909 Return the Euclidean Jordan Algebra corresponding to the set
910 `R^n` under the Hadamard product.
912 Note: this is nothing more than the Cartesian product of ``n``
913 copies of the spin algebra. Once Cartesian product algebras
914 are implemented, this can go.
918 sage: from mjo.eja.eja_algebra import RealCartesianProductEJA
922 This multiplication table can be verified by hand::
924 sage: J = RealCartesianProductEJA(3)
925 sage: e0,e1,e2 = J.gens()
941 We can change the generator prefix::
943 sage: RealCartesianProductEJA(3, prefix='r').gens()
947 def __init__(self
, n
, field
=QQ
, **kwargs
):
948 V
= VectorSpace(field
, n
)
949 mult_table
= [ [ V
.gen(i
)*(i
== j
) for j
in range(n
) ]
952 fdeja
= super(RealCartesianProductEJA
, self
)
953 return fdeja
.__init
__(field
, mult_table
, rank
=n
, **kwargs
)
955 def inner_product(self
, x
, y
):
957 Faster to reimplement than to use natural representations.
961 sage: from mjo.eja.eja_algebra import RealCartesianProductEJA
965 Ensure that this is the usual inner product for the algebras
968 sage: set_random_seed()
969 sage: J = RealCartesianProductEJA.random_instance()
970 sage: x,y = J.random_elements(2)
971 sage: X = x.natural_representation()
972 sage: Y = y.natural_representation()
973 sage: x.inner_product(y) == J.natural_inner_product(X,Y)
977 return x
.to_vector().inner_product(y
.to_vector())
980 def random_eja(field
=QQ
, nontrivial
=False):
982 Return a "random" finite-dimensional Euclidean Jordan Algebra.
986 sage: from mjo.eja.eja_algebra import random_eja
991 Euclidean Jordan algebra of dimension...
994 eja_classes
= KnownRankEJA
.__subclasses
__()
996 eja_classes
.remove(TrivialEJA
)
997 classname
= choice(eja_classes
)
998 return classname
.random_instance(field
=field
)
1005 class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra
):
1007 def _max_test_case_size():
1008 # Play it safe, since this will be squared and the underlying
1009 # field can have dimension 4 (quaternions) too.
1012 def __init__(self
, field
, basis
, rank
, normalize_basis
=True, **kwargs
):
1014 Compared to the superclass constructor, we take a basis instead of
1015 a multiplication table because the latter can be computed in terms
1016 of the former when the product is known (like it is here).
1018 # Used in this class's fast _charpoly_coeff() override.
1019 self
._basis
_normalizers
= None
1021 # We're going to loop through this a few times, so now's a good
1022 # time to ensure that it isn't a generator expression.
1023 basis
= tuple(basis
)
1025 if rank
> 1 and normalize_basis
:
1026 # We'll need sqrt(2) to normalize the basis, and this
1027 # winds up in the multiplication table, so the whole
1028 # algebra needs to be over the field extension.
1029 R
= PolynomialRing(field
, 'z')
1032 if p
.is_irreducible():
1033 field
= field
.extension(p
, 'sqrt2', embedding
=RLF(2).sqrt())
1034 basis
= tuple( s
.change_ring(field
) for s
in basis
)
1035 self
._basis
_normalizers
= tuple(
1036 ~
(self
.natural_inner_product(s
,s
).sqrt()) for s
in basis
)
1037 basis
= tuple(s
*c
for (s
,c
) in zip(basis
,self
._basis
_normalizers
))
1039 Qs
= self
.multiplication_table_from_matrix_basis(basis
)
1041 fdeja
= super(MatrixEuclideanJordanAlgebra
, self
)
1042 return fdeja
.__init
__(field
,
1045 natural_basis
=basis
,
1050 def _charpoly_coeff(self
, i
):
1052 Override the parent method with something that tries to compute
1053 over a faster (non-extension) field.
1055 if self
._basis
_normalizers
is None:
1056 # We didn't normalize, so assume that the basis we started
1057 # with had entries in a nice field.
1058 return super(MatrixEuclideanJordanAlgebra
, self
)._charpoly
_coeff
(i
)
1060 basis
= ( (b
/n
) for (b
,n
) in zip(self
.natural_basis(),
1061 self
._basis
_normalizers
) )
1063 # Do this over the rationals and convert back at the end.
1064 J
= MatrixEuclideanJordanAlgebra(QQ
,
1067 normalize_basis
=False)
1068 (_
,x
,_
,_
) = J
._charpoly
_matrix
_system
()
1069 p
= J
._charpoly
_coeff
(i
)
1070 # p might be missing some vars, have to substitute "optionally"
1071 pairs
= zip(x
.base_ring().gens(), self
._basis
_normalizers
)
1072 substitutions
= { v: v*c for (v,c) in pairs }
1073 result
= p
.subs(substitutions
)
1075 # The result of "subs" can be either a coefficient-ring
1076 # element or a polynomial. Gotta handle both cases.
1078 return self
.base_ring()(result
)
1080 return result
.change_ring(self
.base_ring())
1084 def multiplication_table_from_matrix_basis(basis
):
1086 At least three of the five simple Euclidean Jordan algebras have the
1087 symmetric multiplication (A,B) |-> (AB + BA)/2, where the
1088 multiplication on the right is matrix multiplication. Given a basis
1089 for the underlying matrix space, this function returns a
1090 multiplication table (obtained by looping through the basis
1091 elements) for an algebra of those matrices.
1093 # In S^2, for example, we nominally have four coordinates even
1094 # though the space is of dimension three only. The vector space V
1095 # is supposed to hold the entire long vector, and the subspace W
1096 # of V will be spanned by the vectors that arise from symmetric
1097 # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
1098 field
= basis
[0].base_ring()
1099 dimension
= basis
[0].nrows()
1101 V
= VectorSpace(field
, dimension
**2)
1102 W
= V
.span_of_basis( _mat2vec(s
) for s
in basis
)
1104 mult_table
= [[W
.zero() for j
in range(n
)] for i
in range(n
)]
1107 mat_entry
= (basis
[i
]*basis
[j
] + basis
[j
]*basis
[i
])/2
1108 mult_table
[i
][j
] = W
.coordinate_vector(_mat2vec(mat_entry
))
1116 Embed the matrix ``M`` into a space of real matrices.
1118 The matrix ``M`` can have entries in any field at the moment:
1119 the real numbers, complex numbers, or quaternions. And although
1120 they are not a field, we can probably support octonions at some
1121 point, too. This function returns a real matrix that "acts like"
1122 the original with respect to matrix multiplication; i.e.
1124 real_embed(M*N) = real_embed(M)*real_embed(N)
1127 raise NotImplementedError
1131 def real_unembed(M
):
1133 The inverse of :meth:`real_embed`.
1135 raise NotImplementedError
1139 def natural_inner_product(cls
,X
,Y
):
1140 Xu
= cls
.real_unembed(X
)
1141 Yu
= cls
.real_unembed(Y
)
1142 tr
= (Xu
*Yu
).trace()
1145 # It's real already.
1148 # Otherwise, try the thing that works for complex numbers; and
1149 # if that doesn't work, the thing that works for quaternions.
1151 return tr
.vector()[0] # real part, imag part is index 1
1152 except AttributeError:
1153 # A quaternions doesn't have a vector() method, but does
1154 # have coefficient_tuple() method that returns the
1155 # coefficients of 1, i, j, and k -- in that order.
1156 return tr
.coefficient_tuple()[0]
1159 class RealMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1163 The identity function, for embedding real matrices into real
1169 def real_unembed(M
):
1171 The identity function, for unembedding real matrices from real
1177 class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra
, KnownRankEJA
):
1179 The rank-n simple EJA consisting of real symmetric n-by-n
1180 matrices, the usual symmetric Jordan product, and the trace inner
1181 product. It has dimension `(n^2 + n)/2` over the reals.
1185 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1189 sage: J = RealSymmetricEJA(2)
1190 sage: e0, e1, e2 = J.gens()
1198 In theory, our "field" can be any subfield of the reals::
1200 sage: RealSymmetricEJA(2, AA)
1201 Euclidean Jordan algebra of dimension 3 over Algebraic Real Field
1202 sage: RealSymmetricEJA(2, RR)
1203 Euclidean Jordan algebra of dimension 3 over Real Field with
1204 53 bits of precision
1208 The dimension of this algebra is `(n^2 + n) / 2`::
1210 sage: set_random_seed()
1211 sage: n_max = RealSymmetricEJA._max_test_case_size()
1212 sage: n = ZZ.random_element(1, n_max)
1213 sage: J = RealSymmetricEJA(n)
1214 sage: J.dimension() == (n^2 + n)/2
1217 The Jordan multiplication is what we think it is::
1219 sage: set_random_seed()
1220 sage: J = RealSymmetricEJA.random_instance()
1221 sage: x,y = J.random_elements(2)
1222 sage: actual = (x*y).natural_representation()
1223 sage: X = x.natural_representation()
1224 sage: Y = y.natural_representation()
1225 sage: expected = (X*Y + Y*X)/2
1226 sage: actual == expected
1228 sage: J(expected) == x*y
1231 We can change the generator prefix::
1233 sage: RealSymmetricEJA(3, prefix='q').gens()
1234 (q0, q1, q2, q3, q4, q5)
1236 Our natural basis is normalized with respect to the natural inner
1237 product unless we specify otherwise::
1239 sage: set_random_seed()
1240 sage: J = RealSymmetricEJA.random_instance()
1241 sage: all( b.norm() == 1 for b in J.gens() )
1244 Since our natural basis is normalized with respect to the natural
1245 inner product, and since we know that this algebra is an EJA, any
1246 left-multiplication operator's matrix will be symmetric because
1247 natural->EJA basis representation is an isometry and within the EJA
1248 the operator is self-adjoint by the Jordan axiom::
1250 sage: set_random_seed()
1251 sage: x = RealSymmetricEJA.random_instance().random_element()
1252 sage: x.operator().matrix().is_symmetric()
1257 def _denormalized_basis(cls
, n
, field
):
1259 Return a basis for the space of real symmetric n-by-n matrices.
1263 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1267 sage: set_random_seed()
1268 sage: n = ZZ.random_element(1,5)
1269 sage: B = RealSymmetricEJA._denormalized_basis(n,QQ)
1270 sage: all( M.is_symmetric() for M in B)
1274 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1278 for j
in range(i
+1):
1279 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
1283 Sij
= Eij
+ Eij
.transpose()
1289 def _max_test_case_size():
1290 return 4 # Dimension 10
1293 def __init__(self
, n
, field
=QQ
, **kwargs
):
1294 basis
= self
._denormalized
_basis
(n
, field
)
1295 super(RealSymmetricEJA
, self
).__init
__(field
, basis
, n
, **kwargs
)
1298 class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1302 Embed the n-by-n complex matrix ``M`` into the space of real
1303 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1304 bi` to the block matrix ``[[a,b],[-b,a]]``.
1308 sage: from mjo.eja.eja_algebra import \
1309 ....: ComplexMatrixEuclideanJordanAlgebra
1313 sage: F = QuadraticField(-1, 'i')
1314 sage: x1 = F(4 - 2*i)
1315 sage: x2 = F(1 + 2*i)
1318 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1319 sage: ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1328 Embedding is a homomorphism (isomorphism, in fact)::
1330 sage: set_random_seed()
1331 sage: n_max = ComplexMatrixEuclideanJordanAlgebra._max_test_case_size()
1332 sage: n = ZZ.random_element(n_max)
1333 sage: F = QuadraticField(-1, 'i')
1334 sage: X = random_matrix(F, n)
1335 sage: Y = random_matrix(F, n)
1336 sage: Xe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X)
1337 sage: Ye = ComplexMatrixEuclideanJordanAlgebra.real_embed(Y)
1338 sage: XYe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1345 raise ValueError("the matrix 'M' must be square")
1347 # We don't need any adjoined elements...
1348 field
= M
.base_ring().base_ring()
1352 a
= z
.list()[0] # real part, I guess
1353 b
= z
.list()[1] # imag part, I guess
1354 blocks
.append(matrix(field
, 2, [[a
,b
],[-b
,a
]]))
1356 return matrix
.block(field
, n
, blocks
)
1360 def real_unembed(M
):
1362 The inverse of _embed_complex_matrix().
1366 sage: from mjo.eja.eja_algebra import \
1367 ....: ComplexMatrixEuclideanJordanAlgebra
1371 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1372 ....: [-2, 1, -4, 3],
1373 ....: [ 9, 10, 11, 12],
1374 ....: [-10, 9, -12, 11] ])
1375 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(A)
1377 [ 10*i + 9 12*i + 11]
1381 Unembedding is the inverse of embedding::
1383 sage: set_random_seed()
1384 sage: F = QuadraticField(-1, 'i')
1385 sage: M = random_matrix(F, 3)
1386 sage: Me = ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1387 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1393 raise ValueError("the matrix 'M' must be square")
1394 if not n
.mod(2).is_zero():
1395 raise ValueError("the matrix 'M' must be a complex embedding")
1397 # If "M" was normalized, its base ring might have roots
1398 # adjoined and they can stick around after unembedding.
1399 field
= M
.base_ring()
1400 R
= PolynomialRing(field
, 'z')
1402 F
= field
.extension(z
**2 + 1, 'i', embedding
=CLF(-1).sqrt())
1405 # Go top-left to bottom-right (reading order), converting every
1406 # 2-by-2 block we see to a single complex element.
1408 for k
in range(n
/2):
1409 for j
in range(n
/2):
1410 submat
= M
[2*k
:2*k
+2,2*j
:2*j
+2]
1411 if submat
[0,0] != submat
[1,1]:
1412 raise ValueError('bad on-diagonal submatrix')
1413 if submat
[0,1] != -submat
[1,0]:
1414 raise ValueError('bad off-diagonal submatrix')
1415 z
= submat
[0,0] + submat
[0,1]*i
1418 return matrix(F
, n
/2, elements
)
1422 def natural_inner_product(cls
,X
,Y
):
1424 Compute a natural inner product in this algebra directly from
1429 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1433 This gives the same answer as the slow, default method implemented
1434 in :class:`MatrixEuclideanJordanAlgebra`::
1436 sage: set_random_seed()
1437 sage: J = ComplexHermitianEJA.random_instance()
1438 sage: x,y = J.random_elements(2)
1439 sage: Xe = x.natural_representation()
1440 sage: Ye = y.natural_representation()
1441 sage: X = ComplexHermitianEJA.real_unembed(Xe)
1442 sage: Y = ComplexHermitianEJA.real_unembed(Ye)
1443 sage: expected = (X*Y).trace().vector()[0]
1444 sage: actual = ComplexHermitianEJA.natural_inner_product(Xe,Ye)
1445 sage: actual == expected
1449 return RealMatrixEuclideanJordanAlgebra
.natural_inner_product(X
,Y
)/2
1452 class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra
, KnownRankEJA
):
1454 The rank-n simple EJA consisting of complex Hermitian n-by-n
1455 matrices over the real numbers, the usual symmetric Jordan product,
1456 and the real-part-of-trace inner product. It has dimension `n^2` over
1461 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1465 In theory, our "field" can be any subfield of the reals::
1467 sage: ComplexHermitianEJA(2, AA)
1468 Euclidean Jordan algebra of dimension 4 over Algebraic Real Field
1469 sage: ComplexHermitianEJA(2, RR)
1470 Euclidean Jordan algebra of dimension 4 over Real Field with
1471 53 bits of precision
1475 The dimension of this algebra is `n^2`::
1477 sage: set_random_seed()
1478 sage: n_max = ComplexHermitianEJA._max_test_case_size()
1479 sage: n = ZZ.random_element(1, n_max)
1480 sage: J = ComplexHermitianEJA(n)
1481 sage: J.dimension() == n^2
1484 The Jordan multiplication is what we think it is::
1486 sage: set_random_seed()
1487 sage: J = ComplexHermitianEJA.random_instance()
1488 sage: x,y = J.random_elements(2)
1489 sage: actual = (x*y).natural_representation()
1490 sage: X = x.natural_representation()
1491 sage: Y = y.natural_representation()
1492 sage: expected = (X*Y + Y*X)/2
1493 sage: actual == expected
1495 sage: J(expected) == x*y
1498 We can change the generator prefix::
1500 sage: ComplexHermitianEJA(2, prefix='z').gens()
1503 Our natural basis is normalized with respect to the natural inner
1504 product unless we specify otherwise::
1506 sage: set_random_seed()
1507 sage: J = ComplexHermitianEJA.random_instance()
1508 sage: all( b.norm() == 1 for b in J.gens() )
1511 Since our natural basis is normalized with respect to the natural
1512 inner product, and since we know that this algebra is an EJA, any
1513 left-multiplication operator's matrix will be symmetric because
1514 natural->EJA basis representation is an isometry and within the EJA
1515 the operator is self-adjoint by the Jordan axiom::
1517 sage: set_random_seed()
1518 sage: x = ComplexHermitianEJA.random_instance().random_element()
1519 sage: x.operator().matrix().is_symmetric()
1525 def _denormalized_basis(cls
, n
, field
):
1527 Returns a basis for the space of complex Hermitian n-by-n matrices.
1529 Why do we embed these? Basically, because all of numerical linear
1530 algebra assumes that you're working with vectors consisting of `n`
1531 entries from a field and scalars from the same field. There's no way
1532 to tell SageMath that (for example) the vectors contain complex
1533 numbers, while the scalar field is real.
1537 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1541 sage: set_random_seed()
1542 sage: n = ZZ.random_element(1,5)
1543 sage: field = QuadraticField(2, 'sqrt2')
1544 sage: B = ComplexHermitianEJA._denormalized_basis(n, field)
1545 sage: all( M.is_symmetric() for M in B)
1549 R
= PolynomialRing(field
, 'z')
1551 F
= field
.extension(z
**2 + 1, 'I')
1554 # This is like the symmetric case, but we need to be careful:
1556 # * We want conjugate-symmetry, not just symmetry.
1557 # * The diagonal will (as a result) be real.
1561 for j
in range(i
+1):
1562 Eij
= matrix(F
, n
, lambda k
,l
: k
==i
and l
==j
)
1564 Sij
= cls
.real_embed(Eij
)
1567 # The second one has a minus because it's conjugated.
1568 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
1570 Sij_imag
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
1573 # Since we embedded these, we can drop back to the "field" that we
1574 # started with instead of the complex extension "F".
1575 return ( s
.change_ring(field
) for s
in S
)
1578 def __init__(self
, n
, field
=QQ
, **kwargs
):
1579 basis
= self
._denormalized
_basis
(n
,field
)
1580 super(ComplexHermitianEJA
,self
).__init
__(field
, basis
, n
, **kwargs
)
1583 class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1587 Embed the n-by-n quaternion matrix ``M`` into the space of real
1588 matrices of size 4n-by-4n by first sending each quaternion entry `z
1589 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
1590 c+di],[-c + di, a-bi]]`, and then embedding those into a real
1595 sage: from mjo.eja.eja_algebra import \
1596 ....: QuaternionMatrixEuclideanJordanAlgebra
1600 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1601 sage: i,j,k = Q.gens()
1602 sage: x = 1 + 2*i + 3*j + 4*k
1603 sage: M = matrix(Q, 1, [[x]])
1604 sage: QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1610 Embedding is a homomorphism (isomorphism, in fact)::
1612 sage: set_random_seed()
1613 sage: n_max = QuaternionMatrixEuclideanJordanAlgebra._max_test_case_size()
1614 sage: n = ZZ.random_element(n_max)
1615 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1616 sage: X = random_matrix(Q, n)
1617 sage: Y = random_matrix(Q, n)
1618 sage: Xe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X)
1619 sage: Ye = QuaternionMatrixEuclideanJordanAlgebra.real_embed(Y)
1620 sage: XYe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1625 quaternions
= M
.base_ring()
1628 raise ValueError("the matrix 'M' must be square")
1630 F
= QuadraticField(-1, 'i')
1635 t
= z
.coefficient_tuple()
1640 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
1641 [-c
+ d
*i
, a
- b
*i
]])
1642 realM
= ComplexMatrixEuclideanJordanAlgebra
.real_embed(cplxM
)
1643 blocks
.append(realM
)
1645 # We should have real entries by now, so use the realest field
1646 # we've got for the return value.
1647 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
1652 def real_unembed(M
):
1654 The inverse of _embed_quaternion_matrix().
1658 sage: from mjo.eja.eja_algebra import \
1659 ....: QuaternionMatrixEuclideanJordanAlgebra
1663 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
1664 ....: [-2, 1, -4, 3],
1665 ....: [-3, 4, 1, -2],
1666 ....: [-4, -3, 2, 1]])
1667 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(M)
1668 [1 + 2*i + 3*j + 4*k]
1672 Unembedding is the inverse of embedding::
1674 sage: set_random_seed()
1675 sage: Q = QuaternionAlgebra(QQ, -1, -1)
1676 sage: M = random_matrix(Q, 3)
1677 sage: Me = QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1678 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1684 raise ValueError("the matrix 'M' must be square")
1685 if not n
.mod(4).is_zero():
1686 raise ValueError("the matrix 'M' must be a quaternion embedding")
1688 # Use the base ring of the matrix to ensure that its entries can be
1689 # multiplied by elements of the quaternion algebra.
1690 field
= M
.base_ring()
1691 Q
= QuaternionAlgebra(field
,-1,-1)
1694 # Go top-left to bottom-right (reading order), converting every
1695 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
1698 for l
in range(n
/4):
1699 for m
in range(n
/4):
1700 submat
= ComplexMatrixEuclideanJordanAlgebra
.real_unembed(
1701 M
[4*l
:4*l
+4,4*m
:4*m
+4] )
1702 if submat
[0,0] != submat
[1,1].conjugate():
1703 raise ValueError('bad on-diagonal submatrix')
1704 if submat
[0,1] != -submat
[1,0].conjugate():
1705 raise ValueError('bad off-diagonal submatrix')
1706 z
= submat
[0,0].vector()[0] # real part
1707 z
+= submat
[0,0].vector()[1]*i
# imag part
1708 z
+= submat
[0,1].vector()[0]*j
# real part
1709 z
+= submat
[0,1].vector()[1]*k
# imag part
1712 return matrix(Q
, n
/4, elements
)
1716 def natural_inner_product(cls
,X
,Y
):
1718 Compute a natural inner product in this algebra directly from
1723 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1727 This gives the same answer as the slow, default method implemented
1728 in :class:`MatrixEuclideanJordanAlgebra`::
1730 sage: set_random_seed()
1731 sage: J = QuaternionHermitianEJA.random_instance()
1732 sage: x,y = J.random_elements(2)
1733 sage: Xe = x.natural_representation()
1734 sage: Ye = y.natural_representation()
1735 sage: X = QuaternionHermitianEJA.real_unembed(Xe)
1736 sage: Y = QuaternionHermitianEJA.real_unembed(Ye)
1737 sage: expected = (X*Y).trace().coefficient_tuple()[0]
1738 sage: actual = QuaternionHermitianEJA.natural_inner_product(Xe,Ye)
1739 sage: actual == expected
1743 return RealMatrixEuclideanJordanAlgebra
.natural_inner_product(X
,Y
)/4
1746 class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra
,
1749 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
1750 matrices, the usual symmetric Jordan product, and the
1751 real-part-of-trace inner product. It has dimension `2n^2 - n` over
1756 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1760 In theory, our "field" can be any subfield of the reals::
1762 sage: QuaternionHermitianEJA(2, AA)
1763 Euclidean Jordan algebra of dimension 6 over Algebraic Real Field
1764 sage: QuaternionHermitianEJA(2, RR)
1765 Euclidean Jordan algebra of dimension 6 over Real Field with
1766 53 bits of precision
1770 The dimension of this algebra is `2*n^2 - n`::
1772 sage: set_random_seed()
1773 sage: n_max = QuaternionHermitianEJA._max_test_case_size()
1774 sage: n = ZZ.random_element(1, n_max)
1775 sage: J = QuaternionHermitianEJA(n)
1776 sage: J.dimension() == 2*(n^2) - n
1779 The Jordan multiplication is what we think it is::
1781 sage: set_random_seed()
1782 sage: J = QuaternionHermitianEJA.random_instance()
1783 sage: x,y = J.random_elements(2)
1784 sage: actual = (x*y).natural_representation()
1785 sage: X = x.natural_representation()
1786 sage: Y = y.natural_representation()
1787 sage: expected = (X*Y + Y*X)/2
1788 sage: actual == expected
1790 sage: J(expected) == x*y
1793 We can change the generator prefix::
1795 sage: QuaternionHermitianEJA(2, prefix='a').gens()
1796 (a0, a1, a2, a3, a4, a5)
1798 Our natural basis is normalized with respect to the natural inner
1799 product unless we specify otherwise::
1801 sage: set_random_seed()
1802 sage: J = QuaternionHermitianEJA.random_instance()
1803 sage: all( b.norm() == 1 for b in J.gens() )
1806 Since our natural basis is normalized with respect to the natural
1807 inner product, and since we know that this algebra is an EJA, any
1808 left-multiplication operator's matrix will be symmetric because
1809 natural->EJA basis representation is an isometry and within the EJA
1810 the operator is self-adjoint by the Jordan axiom::
1812 sage: set_random_seed()
1813 sage: x = QuaternionHermitianEJA.random_instance().random_element()
1814 sage: x.operator().matrix().is_symmetric()
1819 def _denormalized_basis(cls
, n
, field
):
1821 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
1823 Why do we embed these? Basically, because all of numerical
1824 linear algebra assumes that you're working with vectors consisting
1825 of `n` entries from a field and scalars from the same field. There's
1826 no way to tell SageMath that (for example) the vectors contain
1827 complex numbers, while the scalar field is real.
1831 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1835 sage: set_random_seed()
1836 sage: n = ZZ.random_element(1,5)
1837 sage: B = QuaternionHermitianEJA._denormalized_basis(n,QQ)
1838 sage: all( M.is_symmetric() for M in B )
1842 Q
= QuaternionAlgebra(QQ
,-1,-1)
1845 # This is like the symmetric case, but we need to be careful:
1847 # * We want conjugate-symmetry, not just symmetry.
1848 # * The diagonal will (as a result) be real.
1852 for j
in range(i
+1):
1853 Eij
= matrix(Q
, n
, lambda k
,l
: k
==i
and l
==j
)
1855 Sij
= cls
.real_embed(Eij
)
1858 # The second, third, and fourth ones have a minus
1859 # because they're conjugated.
1860 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
1862 Sij_I
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
1864 Sij_J
= cls
.real_embed(J
*Eij
- J
*Eij
.transpose())
1866 Sij_K
= cls
.real_embed(K
*Eij
- K
*Eij
.transpose())
1869 # Since we embedded these, we can drop back to the "field" that we
1870 # started with instead of the quaternion algebra "Q".
1871 return ( s
.change_ring(field
) for s
in S
)
1874 def __init__(self
, n
, field
=QQ
, **kwargs
):
1875 basis
= self
._denormalized
_basis
(n
,field
)
1876 super(QuaternionHermitianEJA
,self
).__init
__(field
, basis
, n
, **kwargs
)
1879 class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra
, KnownRankEJA
):
1881 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
1882 with the usual inner product and jordan product ``x*y =
1883 (<x_bar,y_bar>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
1888 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1892 This multiplication table can be verified by hand::
1894 sage: J = JordanSpinEJA(4)
1895 sage: e0,e1,e2,e3 = J.gens()
1911 We can change the generator prefix::
1913 sage: JordanSpinEJA(2, prefix='B').gens()
1917 def __init__(self
, n
, field
=QQ
, **kwargs
):
1918 V
= VectorSpace(field
, n
)
1919 mult_table
= [[V
.zero() for j
in range(n
)] for i
in range(n
)]
1929 z0
= x
.inner_product(y
)
1930 zbar
= y0
*xbar
+ x0
*ybar
1931 z
= V([z0
] + zbar
.list())
1932 mult_table
[i
][j
] = z
1934 # The rank of the spin algebra is two, unless we're in a
1935 # one-dimensional ambient space (because the rank is bounded by
1936 # the ambient dimension).
1937 fdeja
= super(JordanSpinEJA
, self
)
1938 return fdeja
.__init
__(field
, mult_table
, rank
=min(n
,2), **kwargs
)
1940 def inner_product(self
, x
, y
):
1942 Faster to reimplement than to use natural representations.
1946 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1950 Ensure that this is the usual inner product for the algebras
1953 sage: set_random_seed()
1954 sage: J = JordanSpinEJA.random_instance()
1955 sage: x,y = J.random_elements(2)
1956 sage: X = x.natural_representation()
1957 sage: Y = y.natural_representation()
1958 sage: x.inner_product(y) == J.natural_inner_product(X,Y)
1962 return x
.to_vector().inner_product(y
.to_vector())
1965 class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra
, KnownRankEJA
):
1967 The trivial Euclidean Jordan algebra consisting of only a zero element.
1971 sage: from mjo.eja.eja_algebra import TrivialEJA
1975 sage: J = TrivialEJA()
1982 sage: 7*J.one()*12*J.one()
1984 sage: J.one().inner_product(J.one())
1986 sage: J.one().norm()
1988 sage: J.one().subalgebra_generated_by()
1989 Euclidean Jordan algebra of dimension 0 over Rational Field
1994 def __init__(self
, field
=QQ
, **kwargs
):
1996 fdeja
= super(TrivialEJA
, self
)
1997 # The rank is zero using my definition, namely the dimension of the
1998 # largest subalgebra generated by any element.
1999 return fdeja
.__init
__(field
, mult_table
, rank
=0, **kwargs
)