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 izip
, 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
.prandom
import choice
17 from sage
.misc
.table
import table
18 from sage
.modules
.free_module
import FreeModule
, VectorSpace
19 from sage
.rings
.integer_ring
import ZZ
20 from sage
.rings
.number_field
.number_field
import NumberField
, QuadraticField
21 from sage
.rings
.polynomial
.polynomial_ring_constructor
import PolynomialRing
22 from sage
.rings
.rational_field
import QQ
23 from sage
.rings
.real_lazy
import CLF
, RLF
25 from mjo
.eja
.eja_element
import FiniteDimensionalEuclideanJordanAlgebraElement
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
47 sage: from mjo.eja.eja_algebra import random_eja
51 By definition, Jordan multiplication commutes::
53 sage: set_random_seed()
54 sage: J = random_eja()
55 sage: x,y = J.random_elements(2)
61 self
._natural
_basis
= natural_basis
64 category
= MagmaticAlgebras(field
).FiniteDimensional()
65 category
= category
.WithBasis().Unital()
67 fda
= super(FiniteDimensionalEuclideanJordanAlgebra
, self
)
69 range(len(mult_table
)),
72 self
.print_options(bracket
='')
74 # The multiplication table we're given is necessarily in terms
75 # of vectors, because we don't have an algebra yet for
76 # anything to be an element of. However, it's faster in the
77 # long run to have the multiplication table be in terms of
78 # algebra elements. We do this after calling the superclass
79 # constructor so that from_vector() knows what to do.
80 self
._multiplication
_table
= [ map(lambda x
: self
.from_vector(x
), ls
)
81 for ls
in mult_table
]
84 def _element_constructor_(self
, elt
):
86 Construct an element of this algebra from its natural
89 This gets called only after the parent element _call_ method
90 fails to find a coercion for the argument.
94 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
95 ....: RealCartesianProductEJA,
96 ....: RealSymmetricEJA)
100 The identity in `S^n` is converted to the identity in the EJA::
102 sage: J = RealSymmetricEJA(3)
103 sage: I = matrix.identity(QQ,3)
104 sage: J(I) == J.one()
107 This skew-symmetric matrix can't be represented in the EJA::
109 sage: J = RealSymmetricEJA(3)
110 sage: A = matrix(QQ,3, lambda i,j: i-j)
112 Traceback (most recent call last):
114 ArithmeticError: vector is not in free module
118 Ensure that we can convert any element of the two non-matrix
119 simple algebras (whose natural representations are their usual
120 vector representations) back and forth faithfully::
122 sage: set_random_seed()
123 sage: J = RealCartesianProductEJA.random_instance()
124 sage: x = J.random_element()
125 sage: J(x.to_vector().column()) == x
127 sage: J = JordanSpinEJA.random_instance()
128 sage: x = J.random_element()
129 sage: J(x.to_vector().column()) == x
134 # The superclass implementation of random_element()
135 # needs to be able to coerce "0" into the algebra.
138 natural_basis
= self
.natural_basis()
139 basis_space
= natural_basis
[0].matrix_space()
140 if elt
not in basis_space
:
141 raise ValueError("not a naturally-represented algebra element")
143 # Thanks for nothing! Matrix spaces aren't vector spaces in
144 # Sage, so we have to figure out its natural-basis coordinates
145 # ourselves. We use the basis space's ring instead of the
146 # element's ring because the basis space might be an algebraic
147 # closure whereas the base ring of the 3-by-3 identity matrix
148 # could be QQ instead of QQbar.
149 V
= VectorSpace(basis_space
.base_ring(), elt
.nrows()*elt
.ncols())
150 W
= V
.span_of_basis( _mat2vec(s
) for s
in natural_basis
)
151 coords
= W
.coordinate_vector(_mat2vec(elt
))
152 return self
.from_vector(coords
)
157 Return a string representation of ``self``.
161 sage: from mjo.eja.eja_algebra import JordanSpinEJA
165 Ensure that it says what we think it says::
167 sage: JordanSpinEJA(2, field=QQ)
168 Euclidean Jordan algebra of dimension 2 over Rational Field
169 sage: JordanSpinEJA(3, field=RDF)
170 Euclidean Jordan algebra of dimension 3 over Real Double Field
173 fmt
= "Euclidean Jordan algebra of dimension {} over {}"
174 return fmt
.format(self
.dimension(), self
.base_ring())
176 def product_on_basis(self
, i
, j
):
177 return self
._multiplication
_table
[i
][j
]
179 def _a_regular_element(self
):
181 Guess a regular element. Needed to compute the basis for our
182 characteristic polynomial coefficients.
186 sage: from mjo.eja.eja_algebra import random_eja
190 Ensure that this hacky method succeeds for every algebra that we
191 know how to construct::
193 sage: set_random_seed()
194 sage: J = random_eja()
195 sage: J._a_regular_element().is_regular()
200 z
= self
.sum( (i
+1)*gs
[i
] for i
in range(len(gs
)) )
201 if not z
.is_regular():
202 raise ValueError("don't know a regular element")
207 def _charpoly_basis_space(self
):
209 Return the vector space spanned by the basis used in our
210 characteristic polynomial coefficients. This is used not only to
211 compute those coefficients, but also any time we need to
212 evaluate the coefficients (like when we compute the trace or
215 z
= self
._a
_regular
_element
()
216 # Don't use the parent vector space directly here in case this
217 # happens to be a subalgebra. In that case, we would be e.g.
218 # two-dimensional but span_of_basis() would expect three
220 V
= VectorSpace(self
.base_ring(), self
.vector_space().dimension())
221 basis
= [ (z
**k
).to_vector() for k
in range(self
.rank()) ]
222 V1
= V
.span_of_basis( basis
)
223 b
= (V1
.basis() + V1
.complement().basis())
224 return V
.span_of_basis(b
)
229 def _charpoly_coeff(self
, i
):
231 Return the coefficient polynomial "a_{i}" of this algebra's
232 general characteristic polynomial.
234 Having this be a separate cached method lets us compute and
235 store the trace/determinant (a_{r-1} and a_{0} respectively)
236 separate from the entire characteristic polynomial.
238 (A_of_x
, x
, xr
, detA
) = self
._charpoly
_matrix
_system
()
239 R
= A_of_x
.base_ring()
241 # Guaranteed by theory
244 # Danger: the in-place modification is done for performance
245 # reasons (reconstructing a matrix with huge polynomial
246 # entries is slow), but I don't know how cached_method works,
247 # so it's highly possible that we're modifying some global
248 # list variable by reference, here. In other words, you
249 # probably shouldn't call this method twice on the same
250 # algebra, at the same time, in two threads
251 Ai_orig
= A_of_x
.column(i
)
252 A_of_x
.set_column(i
,xr
)
253 numerator
= A_of_x
.det()
254 A_of_x
.set_column(i
,Ai_orig
)
256 # We're relying on the theory here to ensure that each a_i is
257 # indeed back in R, and the added negative signs are to make
258 # the whole charpoly expression sum to zero.
259 return R(-numerator
/detA
)
263 def _charpoly_matrix_system(self
):
265 Compute the matrix whose entries A_ij are polynomials in
266 X1,...,XN, the vector ``x`` of variables X1,...,XN, the vector
267 corresponding to `x^r` and the determinent of the matrix A =
268 [A_ij]. In other words, all of the fixed (cachable) data needed
269 to compute the coefficients of the characteristic polynomial.
274 # Turn my vector space into a module so that "vectors" can
275 # have multivatiate polynomial entries.
276 names
= tuple('X' + str(i
) for i
in range(1,n
+1))
277 R
= PolynomialRing(self
.base_ring(), names
)
279 # Using change_ring() on the parent's vector space doesn't work
280 # here because, in a subalgebra, that vector space has a basis
281 # and change_ring() tries to bring the basis along with it. And
282 # that doesn't work unless the new ring is a PID, which it usually
286 # Now let x = (X1,X2,...,Xn) be the vector whose entries are
290 # And figure out the "left multiplication by x" matrix in
293 monomial_matrices
= [ self
.monomial(i
).operator().matrix()
294 for i
in range(n
) ] # don't recompute these!
296 ek
= self
.monomial(k
).to_vector()
298 sum( x
[i
]*(monomial_matrices
[i
]*ek
)
299 for i
in range(n
) ) )
300 Lx
= matrix
.column(R
, lmbx_cols
)
302 # Now we can compute powers of x "symbolically"
303 x_powers
= [self
.one().to_vector(), x
]
304 for d
in range(2, r
+1):
305 x_powers
.append( Lx
*(x_powers
[-1]) )
307 idmat
= matrix
.identity(R
, n
)
309 W
= self
._charpoly
_basis
_space
()
310 W
= W
.change_ring(R
.fraction_field())
312 # Starting with the standard coordinates x = (X1,X2,...,Xn)
313 # and then converting the entries to W-coordinates allows us
314 # to pass in the standard coordinates to the charpoly and get
315 # back the right answer. Specifically, with x = (X1,X2,...,Xn),
318 # W.coordinates(x^2) eval'd at (standard z-coords)
322 # W-coords of (standard coords of x^2 eval'd at std-coords of z)
324 # We want the middle equivalent thing in our matrix, but use
325 # the first equivalent thing instead so that we can pass in
326 # standard coordinates.
327 x_powers
= [ W
.coordinate_vector(xp
) for xp
in x_powers
]
328 l2
= [idmat
.column(k
-1) for k
in range(r
+1, n
+1)]
329 A_of_x
= matrix
.column(R
, n
, (x_powers
[:r
] + l2
))
330 return (A_of_x
, x
, x_powers
[r
], A_of_x
.det())
334 def characteristic_polynomial(self
):
336 Return a characteristic polynomial that works for all elements
339 The resulting polynomial has `n+1` variables, where `n` is the
340 dimension of this algebra. The first `n` variables correspond to
341 the coordinates of an algebra element: when evaluated at the
342 coordinates of an algebra element with respect to a certain
343 basis, the result is a univariate polynomial (in the one
344 remaining variable ``t``), namely the characteristic polynomial
349 sage: from mjo.eja.eja_algebra import JordanSpinEJA
353 The characteristic polynomial in the spin algebra is given in
354 Alizadeh, Example 11.11::
356 sage: J = JordanSpinEJA(3)
357 sage: p = J.characteristic_polynomial(); p
358 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
359 sage: xvec = J.one().to_vector()
367 # The list of coefficient polynomials a_1, a_2, ..., a_n.
368 a
= [ self
._charpoly
_coeff
(i
) for i
in range(n
) ]
370 # We go to a bit of trouble here to reorder the
371 # indeterminates, so that it's easier to evaluate the
372 # characteristic polynomial at x's coordinates and get back
373 # something in terms of t, which is what we want.
375 S
= PolynomialRing(self
.base_ring(),'t')
377 S
= PolynomialRing(S
, R
.variable_names())
380 # Note: all entries past the rth should be zero. The
381 # coefficient of the highest power (x^r) is 1, but it doesn't
382 # appear in the solution vector which contains coefficients
383 # for the other powers (to make them sum to x^r).
385 a
[r
] = 1 # corresponds to x^r
387 # When the rank is equal to the dimension, trying to
388 # assign a[r] goes out-of-bounds.
389 a
.append(1) # corresponds to x^r
391 return sum( a
[k
]*(t
**k
) for k
in xrange(len(a
)) )
394 def inner_product(self
, x
, y
):
396 The inner product associated with this Euclidean Jordan algebra.
398 Defaults to the trace inner product, but can be overridden by
399 subclasses if they are sure that the necessary properties are
404 sage: from mjo.eja.eja_algebra import random_eja
408 Our inner product is "associative," which means the following for
409 a symmetric bilinear form::
411 sage: set_random_seed()
412 sage: J = random_eja()
413 sage: x,y,z = J.random_elements(3)
414 sage: (x*y).inner_product(z) == y.inner_product(x*z)
418 X
= x
.natural_representation()
419 Y
= y
.natural_representation()
420 return self
.natural_inner_product(X
,Y
)
423 def is_trivial(self
):
425 Return whether or not this algebra is trivial.
427 A trivial algebra contains only the zero element.
431 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
435 sage: J = ComplexHermitianEJA(3)
438 sage: A = J.zero().subalgebra_generated_by()
443 return self
.dimension() == 0
446 def multiplication_table(self
):
448 Return a visual representation of this algebra's multiplication
449 table (on basis elements).
453 sage: from mjo.eja.eja_algebra import JordanSpinEJA
457 sage: J = JordanSpinEJA(4)
458 sage: J.multiplication_table()
459 +----++----+----+----+----+
460 | * || e0 | e1 | e2 | e3 |
461 +====++====+====+====+====+
462 | e0 || e0 | e1 | e2 | e3 |
463 +----++----+----+----+----+
464 | e1 || e1 | e0 | 0 | 0 |
465 +----++----+----+----+----+
466 | e2 || e2 | 0 | e0 | 0 |
467 +----++----+----+----+----+
468 | e3 || e3 | 0 | 0 | e0 |
469 +----++----+----+----+----+
472 M
= list(self
._multiplication
_table
) # copy
473 for i
in xrange(len(M
)):
474 # M had better be "square"
475 M
[i
] = [self
.monomial(i
)] + M
[i
]
476 M
= [["*"] + list(self
.gens())] + M
477 return table(M
, header_row
=True, header_column
=True, frame
=True)
480 def natural_basis(self
):
482 Return a more-natural representation of this algebra's basis.
484 Every finite-dimensional Euclidean Jordan Algebra is a direct
485 sum of five simple algebras, four of which comprise Hermitian
486 matrices. This method returns the original "natural" basis
487 for our underlying vector space. (Typically, the natural basis
488 is used to construct the multiplication table in the first place.)
490 Note that this will always return a matrix. The standard basis
491 in `R^n` will be returned as `n`-by-`1` column matrices.
495 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
496 ....: RealSymmetricEJA)
500 sage: J = RealSymmetricEJA(2)
502 Finite family {0: e0, 1: e1, 2: e2}
503 sage: J.natural_basis()
505 [1 0] [ 0 1/2*sqrt2] [0 0]
506 [0 0], [1/2*sqrt2 0], [0 1]
511 sage: J = JordanSpinEJA(2)
513 Finite family {0: e0, 1: e1}
514 sage: J.natural_basis()
521 if self
._natural
_basis
is None:
522 M
= self
.natural_basis_space()
523 return tuple( M(b
.to_vector()) for b
in self
.basis() )
525 return self
._natural
_basis
528 def natural_basis_space(self
):
530 Return the matrix space in which this algebra's natural basis
533 if self
._natural
_basis
is None or len(self
._natural
_basis
) == 0:
534 return MatrixSpace(self
.base_ring(), self
.dimension(), 1)
536 return self
._natural
_basis
[0].matrix_space()
540 def natural_inner_product(X
,Y
):
542 Compute the inner product of two naturally-represented elements.
544 For example in the real symmetric matrix EJA, this will compute
545 the trace inner-product of two n-by-n symmetric matrices. The
546 default should work for the real cartesian product EJA, the
547 Jordan spin EJA, and the real symmetric matrices. The others
548 will have to be overridden.
550 return (X
.conjugate_transpose()*Y
).trace()
556 Return the unit element of this algebra.
560 sage: from mjo.eja.eja_algebra import (RealCartesianProductEJA,
565 sage: J = RealCartesianProductEJA(5)
567 e0 + e1 + e2 + e3 + e4
571 The identity element acts like the identity::
573 sage: set_random_seed()
574 sage: J = random_eja()
575 sage: x = J.random_element()
576 sage: J.one()*x == x and x*J.one() == x
579 The matrix of the unit element's operator is the identity::
581 sage: set_random_seed()
582 sage: J = random_eja()
583 sage: actual = J.one().operator().matrix()
584 sage: expected = matrix.identity(J.base_ring(), J.dimension())
585 sage: actual == expected
589 # We can brute-force compute the matrices of the operators
590 # that correspond to the basis elements of this algebra.
591 # If some linear combination of those basis elements is the
592 # algebra identity, then the same linear combination of
593 # their matrices has to be the identity matrix.
595 # Of course, matrices aren't vectors in sage, so we have to
596 # appeal to the "long vectors" isometry.
597 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
599 # Now we use basis linear algebra to find the coefficients,
600 # of the matrices-as-vectors-linear-combination, which should
601 # work for the original algebra basis too.
602 A
= matrix
.column(self
.base_ring(), oper_vecs
)
604 # We used the isometry on the left-hand side already, but we
605 # still need to do it for the right-hand side. Recall that we
606 # wanted something that summed to the identity matrix.
607 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
609 # Now if there's an identity element in the algebra, this should work.
610 coeffs
= A
.solve_right(b
)
611 return self
.linear_combination(zip(self
.gens(), coeffs
))
614 def random_element(self
):
615 # Temporary workaround for https://trac.sagemath.org/ticket/28327
616 if self
.is_trivial():
619 s
= super(FiniteDimensionalEuclideanJordanAlgebra
, self
)
620 return s
.random_element()
622 def random_elements(self
, count
):
624 Return ``count`` random elements as a tuple.
628 sage: from mjo.eja.eja_algebra import JordanSpinEJA
632 sage: J = JordanSpinEJA(3)
633 sage: x,y,z = J.random_elements(3)
634 sage: all( [ x in J, y in J, z in J ])
636 sage: len( J.random_elements(10) ) == 10
640 return tuple( self
.random_element() for idx
in xrange(count
) )
645 Return the rank of this EJA.
649 The author knows of no algorithm to compute the rank of an EJA
650 where only the multiplication table is known. In lieu of one, we
651 require the rank to be specified when the algebra is created,
652 and simply pass along that number here.
656 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
657 ....: RealSymmetricEJA,
658 ....: ComplexHermitianEJA,
659 ....: QuaternionHermitianEJA,
664 The rank of the Jordan spin algebra is always two::
666 sage: JordanSpinEJA(2).rank()
668 sage: JordanSpinEJA(3).rank()
670 sage: JordanSpinEJA(4).rank()
673 The rank of the `n`-by-`n` Hermitian real, complex, or
674 quaternion matrices is `n`::
676 sage: RealSymmetricEJA(4).rank()
678 sage: ComplexHermitianEJA(3).rank()
680 sage: QuaternionHermitianEJA(2).rank()
685 Ensure that every EJA that we know how to construct has a
686 positive integer rank::
688 sage: set_random_seed()
689 sage: r = random_eja().rank()
690 sage: r in ZZ and r > 0
697 def vector_space(self
):
699 Return the vector space that underlies this algebra.
703 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
707 sage: J = RealSymmetricEJA(2)
708 sage: J.vector_space()
709 Vector space of dimension 3 over...
712 return self
.zero().to_vector().parent().ambient_vector_space()
715 Element
= FiniteDimensionalEuclideanJordanAlgebraElement
718 class KnownRankEJA(object):
720 A class for algebras that we actually know we can construct. The
721 main issue is that, for most of our methods to make sense, we need
722 to know the rank of our algebra. Thus we can't simply generate a
723 "random" algebra, or even check that a given basis and product
724 satisfy the axioms; because even if everything looks OK, we wouldn't
725 know the rank we need to actuallty build the thing.
727 Not really a subclass of FDEJA because doing that causes method
728 resolution errors, e.g.
730 TypeError: Error when calling the metaclass bases
731 Cannot create a consistent method resolution
732 order (MRO) for bases FiniteDimensionalEuclideanJordanAlgebra,
737 def _max_test_case_size():
739 Return an integer "size" that is an upper bound on the size of
740 this algebra when it is used in a random test
741 case. Unfortunately, the term "size" is quite vague -- when
742 dealing with `R^n` under either the Hadamard or Jordan spin
743 product, the "size" refers to the dimension `n`. When dealing
744 with a matrix algebra (real symmetric or complex/quaternion
745 Hermitian), it refers to the size of the matrix, which is
746 far less than the dimension of the underlying vector space.
748 We default to five in this class, which is safe in `R^n`. The
749 matrix algebra subclasses (or any class where the "size" is
750 interpreted to be far less than the dimension) should override
751 with a smaller number.
756 def random_instance(cls
, field
=QQ
, **kwargs
):
758 Return a random instance of this type of algebra.
760 Beware, this will crash for "most instances" because the
761 constructor below looks wrong.
763 n
= ZZ
.random_element(cls
._max
_test
_case
_size
()) + 1
764 return cls(n
, field
, **kwargs
)
767 class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra
,
770 Return the Euclidean Jordan Algebra corresponding to the set
771 `R^n` under the Hadamard product.
773 Note: this is nothing more than the Cartesian product of ``n``
774 copies of the spin algebra. Once Cartesian product algebras
775 are implemented, this can go.
779 sage: from mjo.eja.eja_algebra import RealCartesianProductEJA
783 This multiplication table can be verified by hand::
785 sage: J = RealCartesianProductEJA(3)
786 sage: e0,e1,e2 = J.gens()
802 We can change the generator prefix::
804 sage: RealCartesianProductEJA(3, prefix='r').gens()
808 def __init__(self
, n
, field
=QQ
, **kwargs
):
809 V
= VectorSpace(field
, n
)
810 mult_table
= [ [ V
.gen(i
)*(i
== j
) for j
in xrange(n
) ]
813 fdeja
= super(RealCartesianProductEJA
, self
)
814 return fdeja
.__init
__(field
, mult_table
, rank
=n
, **kwargs
)
816 def inner_product(self
, x
, y
):
818 Faster to reimplement than to use natural representations.
822 sage: from mjo.eja.eja_algebra import RealCartesianProductEJA
826 Ensure that this is the usual inner product for the algebras
829 sage: set_random_seed()
830 sage: J = RealCartesianProductEJA.random_instance()
831 sage: x,y = J.random_elements(2)
832 sage: X = x.natural_representation()
833 sage: Y = y.natural_representation()
834 sage: x.inner_product(y) == J.natural_inner_product(X,Y)
838 return x
.to_vector().inner_product(y
.to_vector())
843 Return a "random" finite-dimensional Euclidean Jordan Algebra.
847 For now, we choose a random natural number ``n`` (greater than zero)
848 and then give you back one of the following:
850 * The cartesian product of the rational numbers ``n`` times; this is
851 ``QQ^n`` with the Hadamard product.
853 * The Jordan spin algebra on ``QQ^n``.
855 * The ``n``-by-``n`` rational symmetric matrices with the symmetric
858 * The ``n``-by-``n`` complex-rational Hermitian matrices embedded
859 in the space of ``2n``-by-``2n`` real symmetric matrices.
861 * The ``n``-by-``n`` quaternion-rational Hermitian matrices embedded
862 in the space of ``4n``-by-``4n`` real symmetric matrices.
864 Later this might be extended to return Cartesian products of the
869 sage: from mjo.eja.eja_algebra import random_eja
874 Euclidean Jordan algebra of dimension...
877 classname
= choice(KnownRankEJA
.__subclasses
__())
878 return classname
.random_instance()
885 class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra
):
887 def _max_test_case_size():
888 # Play it safe, since this will be squared and the underlying
889 # field can have dimension 4 (quaternions) too.
892 def __init__(self
, field
, basis
, rank
, normalize_basis
=True, **kwargs
):
894 Compared to the superclass constructor, we take a basis instead of
895 a multiplication table because the latter can be computed in terms
896 of the former when the product is known (like it is here).
898 # Used in this class's fast _charpoly_coeff() override.
899 self
._basis
_normalizers
= None
901 # We're going to loop through this a few times, so now's a good
902 # time to ensure that it isn't a generator expression.
905 if rank
> 1 and normalize_basis
:
906 # We'll need sqrt(2) to normalize the basis, and this
907 # winds up in the multiplication table, so the whole
908 # algebra needs to be over the field extension.
909 R
= PolynomialRing(field
, 'z')
912 if p
.is_irreducible():
913 field
= NumberField(p
, 'sqrt2', embedding
=RLF(2).sqrt())
914 basis
= tuple( s
.change_ring(field
) for s
in basis
)
915 self
._basis
_normalizers
= tuple(
916 ~
(self
.natural_inner_product(s
,s
).sqrt()) for s
in basis
)
917 basis
= tuple(s
*c
for (s
,c
) in izip(basis
,self
._basis
_normalizers
))
919 Qs
= self
.multiplication_table_from_matrix_basis(basis
)
921 fdeja
= super(MatrixEuclideanJordanAlgebra
, self
)
922 return fdeja
.__init
__(field
,
930 def _charpoly_coeff(self
, i
):
932 Override the parent method with something that tries to compute
933 over a faster (non-extension) field.
935 if self
._basis
_normalizers
is None:
936 # We didn't normalize, so assume that the basis we started
937 # with had entries in a nice field.
938 return super(MatrixEuclideanJordanAlgebra
, self
)._charpoly
_coeff
(i
)
940 basis
= ( (b
/n
) for (b
,n
) in izip(self
.natural_basis(),
941 self
._basis
_normalizers
) )
943 # Do this over the rationals and convert back at the end.
944 J
= MatrixEuclideanJordanAlgebra(QQ
,
947 normalize_basis
=False)
948 (_
,x
,_
,_
) = J
._charpoly
_matrix
_system
()
949 p
= J
._charpoly
_coeff
(i
)
950 # p might be missing some vars, have to substitute "optionally"
951 pairs
= izip(x
.base_ring().gens(), self
._basis
_normalizers
)
952 substitutions
= { v: v*c for (v,c) in pairs }
953 result
= p
.subs(substitutions
)
955 # The result of "subs" can be either a coefficient-ring
956 # element or a polynomial. Gotta handle both cases.
958 return self
.base_ring()(result
)
960 return result
.change_ring(self
.base_ring())
964 def multiplication_table_from_matrix_basis(basis
):
966 At least three of the five simple Euclidean Jordan algebras have the
967 symmetric multiplication (A,B) |-> (AB + BA)/2, where the
968 multiplication on the right is matrix multiplication. Given a basis
969 for the underlying matrix space, this function returns a
970 multiplication table (obtained by looping through the basis
971 elements) for an algebra of those matrices.
973 # In S^2, for example, we nominally have four coordinates even
974 # though the space is of dimension three only. The vector space V
975 # is supposed to hold the entire long vector, and the subspace W
976 # of V will be spanned by the vectors that arise from symmetric
977 # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
978 field
= basis
[0].base_ring()
979 dimension
= basis
[0].nrows()
981 V
= VectorSpace(field
, dimension
**2)
982 W
= V
.span_of_basis( _mat2vec(s
) for s
in basis
)
984 mult_table
= [[W
.zero() for j
in xrange(n
)] for i
in xrange(n
)]
987 mat_entry
= (basis
[i
]*basis
[j
] + basis
[j
]*basis
[i
])/2
988 mult_table
[i
][j
] = W
.coordinate_vector(_mat2vec(mat_entry
))
996 Embed the matrix ``M`` into a space of real matrices.
998 The matrix ``M`` can have entries in any field at the moment:
999 the real numbers, complex numbers, or quaternions. And although
1000 they are not a field, we can probably support octonions at some
1001 point, too. This function returns a real matrix that "acts like"
1002 the original with respect to matrix multiplication; i.e.
1004 real_embed(M*N) = real_embed(M)*real_embed(N)
1007 raise NotImplementedError
1011 def real_unembed(M
):
1013 The inverse of :meth:`real_embed`.
1015 raise NotImplementedError
1019 def natural_inner_product(cls
,X
,Y
):
1020 Xu
= cls
.real_unembed(X
)
1021 Yu
= cls
.real_unembed(Y
)
1022 tr
= (Xu
*Yu
).trace()
1024 # It's real already.
1027 # Otherwise, try the thing that works for complex numbers; and
1028 # if that doesn't work, the thing that works for quaternions.
1030 return tr
.vector()[0] # real part, imag part is index 1
1031 except AttributeError:
1032 # A quaternions doesn't have a vector() method, but does
1033 # have coefficient_tuple() method that returns the
1034 # coefficients of 1, i, j, and k -- in that order.
1035 return tr
.coefficient_tuple()[0]
1038 class RealMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1042 The identity function, for embedding real matrices into real
1048 def real_unembed(M
):
1050 The identity function, for unembedding real matrices from real
1056 class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra
, KnownRankEJA
):
1058 The rank-n simple EJA consisting of real symmetric n-by-n
1059 matrices, the usual symmetric Jordan product, and the trace inner
1060 product. It has dimension `(n^2 + n)/2` over the reals.
1064 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1068 sage: J = RealSymmetricEJA(2)
1069 sage: e0, e1, e2 = J.gens()
1079 The dimension of this algebra is `(n^2 + n) / 2`::
1081 sage: set_random_seed()
1082 sage: n_max = RealSymmetricEJA._max_test_case_size()
1083 sage: n = ZZ.random_element(1, n_max)
1084 sage: J = RealSymmetricEJA(n)
1085 sage: J.dimension() == (n^2 + n)/2
1088 The Jordan multiplication is what we think it is::
1090 sage: set_random_seed()
1091 sage: J = RealSymmetricEJA.random_instance()
1092 sage: x,y = J.random_elements(2)
1093 sage: actual = (x*y).natural_representation()
1094 sage: X = x.natural_representation()
1095 sage: Y = y.natural_representation()
1096 sage: expected = (X*Y + Y*X)/2
1097 sage: actual == expected
1099 sage: J(expected) == x*y
1102 We can change the generator prefix::
1104 sage: RealSymmetricEJA(3, prefix='q').gens()
1105 (q0, q1, q2, q3, q4, q5)
1107 Our natural basis is normalized with respect to the natural inner
1108 product unless we specify otherwise::
1110 sage: set_random_seed()
1111 sage: J = RealSymmetricEJA.random_instance()
1112 sage: all( b.norm() == 1 for b in J.gens() )
1115 Since our natural basis is normalized with respect to the natural
1116 inner product, and since we know that this algebra is an EJA, any
1117 left-multiplication operator's matrix will be symmetric because
1118 natural->EJA basis representation is an isometry and within the EJA
1119 the operator is self-adjoint by the Jordan axiom::
1121 sage: set_random_seed()
1122 sage: x = RealSymmetricEJA.random_instance().random_element()
1123 sage: x.operator().matrix().is_symmetric()
1128 def _denormalized_basis(cls
, n
, field
):
1130 Return a basis for the space of real symmetric n-by-n matrices.
1134 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1138 sage: set_random_seed()
1139 sage: n = ZZ.random_element(1,5)
1140 sage: B = RealSymmetricEJA._denormalized_basis(n,QQ)
1141 sage: all( M.is_symmetric() for M in B)
1145 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1149 for j
in xrange(i
+1):
1150 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
1154 Sij
= Eij
+ Eij
.transpose()
1160 def _max_test_case_size():
1161 return 4 # Dimension 10
1164 def __init__(self
, n
, field
=QQ
, **kwargs
):
1165 basis
= self
._denormalized
_basis
(n
, field
)
1166 super(RealSymmetricEJA
, self
).__init
__(field
, basis
, n
, **kwargs
)
1169 class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1173 Embed the n-by-n complex matrix ``M`` into the space of real
1174 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1175 bi` to the block matrix ``[[a,b],[-b,a]]``.
1179 sage: from mjo.eja.eja_algebra import \
1180 ....: ComplexMatrixEuclideanJordanAlgebra
1184 sage: F = QuadraticField(-1, 'i')
1185 sage: x1 = F(4 - 2*i)
1186 sage: x2 = F(1 + 2*i)
1189 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1190 sage: ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1199 Embedding is a homomorphism (isomorphism, in fact)::
1201 sage: set_random_seed()
1202 sage: n_max = ComplexMatrixEuclideanJordanAlgebra._max_test_case_size()
1203 sage: n = ZZ.random_element(n_max)
1204 sage: F = QuadraticField(-1, 'i')
1205 sage: X = random_matrix(F, n)
1206 sage: Y = random_matrix(F, n)
1207 sage: Xe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X)
1208 sage: Ye = ComplexMatrixEuclideanJordanAlgebra.real_embed(Y)
1209 sage: XYe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1216 raise ValueError("the matrix 'M' must be square")
1217 field
= M
.base_ring()
1220 a
= z
.vector()[0] # real part, I guess
1221 b
= z
.vector()[1] # imag part, I guess
1222 blocks
.append(matrix(field
, 2, [[a
,b
],[-b
,a
]]))
1224 # We can drop the imaginaries here.
1225 return matrix
.block(field
.base_ring(), n
, blocks
)
1229 def real_unembed(M
):
1231 The inverse of _embed_complex_matrix().
1235 sage: from mjo.eja.eja_algebra import \
1236 ....: ComplexMatrixEuclideanJordanAlgebra
1240 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1241 ....: [-2, 1, -4, 3],
1242 ....: [ 9, 10, 11, 12],
1243 ....: [-10, 9, -12, 11] ])
1244 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(A)
1246 [ 10*i + 9 12*i + 11]
1250 Unembedding is the inverse of embedding::
1252 sage: set_random_seed()
1253 sage: F = QuadraticField(-1, 'i')
1254 sage: M = random_matrix(F, 3)
1255 sage: Me = ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1256 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1262 raise ValueError("the matrix 'M' must be square")
1263 if not n
.mod(2).is_zero():
1264 raise ValueError("the matrix 'M' must be a complex embedding")
1266 field
= M
.base_ring() # This should already have sqrt2
1267 R
= PolynomialRing(field
, 'z')
1269 F
= NumberField(z
**2 + 1,'i', embedding
=CLF(-1).sqrt())
1272 # Go top-left to bottom-right (reading order), converting every
1273 # 2-by-2 block we see to a single complex element.
1275 for k
in xrange(n
/2):
1276 for j
in xrange(n
/2):
1277 submat
= M
[2*k
:2*k
+2,2*j
:2*j
+2]
1278 if submat
[0,0] != submat
[1,1]:
1279 raise ValueError('bad on-diagonal submatrix')
1280 if submat
[0,1] != -submat
[1,0]:
1281 raise ValueError('bad off-diagonal submatrix')
1282 z
= submat
[0,0] + submat
[0,1]*i
1285 return matrix(F
, n
/2, elements
)
1289 def natural_inner_product(cls
,X
,Y
):
1291 Compute a natural inner product in this algebra directly from
1296 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1300 This gives the same answer as the slow, default method implemented
1301 in :class:`MatrixEuclideanJordanAlgebra`::
1303 sage: set_random_seed()
1304 sage: J = ComplexHermitianEJA.random_instance()
1305 sage: x,y = J.random_elements(2)
1306 sage: Xe = x.natural_representation()
1307 sage: Ye = y.natural_representation()
1308 sage: X = ComplexHermitianEJA.real_unembed(Xe)
1309 sage: Y = ComplexHermitianEJA.real_unembed(Ye)
1310 sage: expected = (X*Y).trace().vector()[0]
1311 sage: actual = ComplexHermitianEJA.natural_inner_product(Xe,Ye)
1312 sage: actual == expected
1316 return RealMatrixEuclideanJordanAlgebra
.natural_inner_product(X
,Y
)/2
1319 class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra
, KnownRankEJA
):
1321 The rank-n simple EJA consisting of complex Hermitian n-by-n
1322 matrices over the real numbers, the usual symmetric Jordan product,
1323 and the real-part-of-trace inner product. It has dimension `n^2` over
1328 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1332 The dimension of this algebra is `n^2`::
1334 sage: set_random_seed()
1335 sage: n_max = ComplexHermitianEJA._max_test_case_size()
1336 sage: n = ZZ.random_element(1, n_max)
1337 sage: J = ComplexHermitianEJA(n)
1338 sage: J.dimension() == n^2
1341 The Jordan multiplication is what we think it is::
1343 sage: set_random_seed()
1344 sage: J = ComplexHermitianEJA.random_instance()
1345 sage: x,y = J.random_elements(2)
1346 sage: actual = (x*y).natural_representation()
1347 sage: X = x.natural_representation()
1348 sage: Y = y.natural_representation()
1349 sage: expected = (X*Y + Y*X)/2
1350 sage: actual == expected
1352 sage: J(expected) == x*y
1355 We can change the generator prefix::
1357 sage: ComplexHermitianEJA(2, prefix='z').gens()
1360 Our natural basis is normalized with respect to the natural inner
1361 product unless we specify otherwise::
1363 sage: set_random_seed()
1364 sage: J = ComplexHermitianEJA.random_instance()
1365 sage: all( b.norm() == 1 for b in J.gens() )
1368 Since our natural basis is normalized with respect to the natural
1369 inner product, and since we know that this algebra is an EJA, any
1370 left-multiplication operator's matrix will be symmetric because
1371 natural->EJA basis representation is an isometry and within the EJA
1372 the operator is self-adjoint by the Jordan axiom::
1374 sage: set_random_seed()
1375 sage: x = ComplexHermitianEJA.random_instance().random_element()
1376 sage: x.operator().matrix().is_symmetric()
1382 def _denormalized_basis(cls
, n
, field
):
1384 Returns a basis for the space of complex Hermitian n-by-n matrices.
1386 Why do we embed these? Basically, because all of numerical linear
1387 algebra assumes that you're working with vectors consisting of `n`
1388 entries from a field and scalars from the same field. There's no way
1389 to tell SageMath that (for example) the vectors contain complex
1390 numbers, while the scalar field is real.
1394 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1398 sage: set_random_seed()
1399 sage: n = ZZ.random_element(1,5)
1400 sage: field = QuadraticField(2, 'sqrt2')
1401 sage: B = ComplexHermitianEJA._denormalized_basis(n, field)
1402 sage: all( M.is_symmetric() for M in B)
1406 R
= PolynomialRing(field
, 'z')
1408 F
= NumberField(z
**2 + 1, 'I', embedding
=CLF(-1).sqrt())
1411 # This is like the symmetric case, but we need to be careful:
1413 # * We want conjugate-symmetry, not just symmetry.
1414 # * The diagonal will (as a result) be real.
1418 for j
in xrange(i
+1):
1419 Eij
= matrix(F
, n
, lambda k
,l
: k
==i
and l
==j
)
1421 Sij
= cls
.real_embed(Eij
)
1424 # The second one has a minus because it's conjugated.
1425 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
1427 Sij_imag
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
1430 # Since we embedded these, we can drop back to the "field" that we
1431 # started with instead of the complex extension "F".
1432 return ( s
.change_ring(field
) for s
in S
)
1435 def __init__(self
, n
, field
=QQ
, **kwargs
):
1436 basis
= self
._denormalized
_basis
(n
,field
)
1437 super(ComplexHermitianEJA
,self
).__init
__(field
, basis
, n
, **kwargs
)
1440 class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1444 Embed the n-by-n quaternion matrix ``M`` into the space of real
1445 matrices of size 4n-by-4n by first sending each quaternion entry `z
1446 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
1447 c+di],[-c + di, a-bi]]`, and then embedding those into a real
1452 sage: from mjo.eja.eja_algebra import \
1453 ....: QuaternionMatrixEuclideanJordanAlgebra
1457 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1458 sage: i,j,k = Q.gens()
1459 sage: x = 1 + 2*i + 3*j + 4*k
1460 sage: M = matrix(Q, 1, [[x]])
1461 sage: QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1467 Embedding is a homomorphism (isomorphism, in fact)::
1469 sage: set_random_seed()
1470 sage: n_max = QuaternionMatrixEuclideanJordanAlgebra._max_test_case_size()
1471 sage: n = ZZ.random_element(n_max)
1472 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1473 sage: X = random_matrix(Q, n)
1474 sage: Y = random_matrix(Q, n)
1475 sage: Xe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X)
1476 sage: Ye = QuaternionMatrixEuclideanJordanAlgebra.real_embed(Y)
1477 sage: XYe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1482 quaternions
= M
.base_ring()
1485 raise ValueError("the matrix 'M' must be square")
1487 F
= QuadraticField(-1, 'i')
1492 t
= z
.coefficient_tuple()
1497 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
1498 [-c
+ d
*i
, a
- b
*i
]])
1499 realM
= ComplexMatrixEuclideanJordanAlgebra
.real_embed(cplxM
)
1500 blocks
.append(realM
)
1502 # We should have real entries by now, so use the realest field
1503 # we've got for the return value.
1504 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
1509 def real_unembed(M
):
1511 The inverse of _embed_quaternion_matrix().
1515 sage: from mjo.eja.eja_algebra import \
1516 ....: QuaternionMatrixEuclideanJordanAlgebra
1520 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
1521 ....: [-2, 1, -4, 3],
1522 ....: [-3, 4, 1, -2],
1523 ....: [-4, -3, 2, 1]])
1524 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(M)
1525 [1 + 2*i + 3*j + 4*k]
1529 Unembedding is the inverse of embedding::
1531 sage: set_random_seed()
1532 sage: Q = QuaternionAlgebra(QQ, -1, -1)
1533 sage: M = random_matrix(Q, 3)
1534 sage: Me = QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1535 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1541 raise ValueError("the matrix 'M' must be square")
1542 if not n
.mod(4).is_zero():
1543 raise ValueError("the matrix 'M' must be a quaternion embedding")
1545 # Use the base ring of the matrix to ensure that its entries can be
1546 # multiplied by elements of the quaternion algebra.
1547 field
= M
.base_ring()
1548 Q
= QuaternionAlgebra(field
,-1,-1)
1551 # Go top-left to bottom-right (reading order), converting every
1552 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
1555 for l
in xrange(n
/4):
1556 for m
in xrange(n
/4):
1557 submat
= ComplexMatrixEuclideanJordanAlgebra
.real_unembed(
1558 M
[4*l
:4*l
+4,4*m
:4*m
+4] )
1559 if submat
[0,0] != submat
[1,1].conjugate():
1560 raise ValueError('bad on-diagonal submatrix')
1561 if submat
[0,1] != -submat
[1,0].conjugate():
1562 raise ValueError('bad off-diagonal submatrix')
1563 z
= submat
[0,0].vector()[0] # real part
1564 z
+= submat
[0,0].vector()[1]*i
# imag part
1565 z
+= submat
[0,1].vector()[0]*j
# real part
1566 z
+= submat
[0,1].vector()[1]*k
# imag part
1569 return matrix(Q
, n
/4, elements
)
1573 def natural_inner_product(cls
,X
,Y
):
1575 Compute a natural inner product in this algebra directly from
1580 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1584 This gives the same answer as the slow, default method implemented
1585 in :class:`MatrixEuclideanJordanAlgebra`::
1587 sage: set_random_seed()
1588 sage: J = QuaternionHermitianEJA.random_instance()
1589 sage: x,y = J.random_elements(2)
1590 sage: Xe = x.natural_representation()
1591 sage: Ye = y.natural_representation()
1592 sage: X = QuaternionHermitianEJA.real_unembed(Xe)
1593 sage: Y = QuaternionHermitianEJA.real_unembed(Ye)
1594 sage: expected = (X*Y).trace().coefficient_tuple()[0]
1595 sage: actual = QuaternionHermitianEJA.natural_inner_product(Xe,Ye)
1596 sage: actual == expected
1600 return RealMatrixEuclideanJordanAlgebra
.natural_inner_product(X
,Y
)/4
1603 class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra
,
1606 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
1607 matrices, the usual symmetric Jordan product, and the
1608 real-part-of-trace inner product. It has dimension `2n^2 - n` over
1613 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1617 The dimension of this algebra is `2*n^2 - n`::
1619 sage: set_random_seed()
1620 sage: n_max = QuaternionHermitianEJA._max_test_case_size()
1621 sage: n = ZZ.random_element(1, n_max)
1622 sage: J = QuaternionHermitianEJA(n)
1623 sage: J.dimension() == 2*(n^2) - n
1626 The Jordan multiplication is what we think it is::
1628 sage: set_random_seed()
1629 sage: J = QuaternionHermitianEJA.random_instance()
1630 sage: x,y = J.random_elements(2)
1631 sage: actual = (x*y).natural_representation()
1632 sage: X = x.natural_representation()
1633 sage: Y = y.natural_representation()
1634 sage: expected = (X*Y + Y*X)/2
1635 sage: actual == expected
1637 sage: J(expected) == x*y
1640 We can change the generator prefix::
1642 sage: QuaternionHermitianEJA(2, prefix='a').gens()
1643 (a0, a1, a2, a3, a4, a5)
1645 Our natural basis is normalized with respect to the natural inner
1646 product unless we specify otherwise::
1648 sage: set_random_seed()
1649 sage: J = QuaternionHermitianEJA.random_instance()
1650 sage: all( b.norm() == 1 for b in J.gens() )
1653 Since our natural basis is normalized with respect to the natural
1654 inner product, and since we know that this algebra is an EJA, any
1655 left-multiplication operator's matrix will be symmetric because
1656 natural->EJA basis representation is an isometry and within the EJA
1657 the operator is self-adjoint by the Jordan axiom::
1659 sage: set_random_seed()
1660 sage: x = QuaternionHermitianEJA.random_instance().random_element()
1661 sage: x.operator().matrix().is_symmetric()
1666 def _denormalized_basis(cls
, n
, field
):
1668 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
1670 Why do we embed these? Basically, because all of numerical
1671 linear algebra assumes that you're working with vectors consisting
1672 of `n` entries from a field and scalars from the same field. There's
1673 no way to tell SageMath that (for example) the vectors contain
1674 complex numbers, while the scalar field is real.
1678 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1682 sage: set_random_seed()
1683 sage: n = ZZ.random_element(1,5)
1684 sage: B = QuaternionHermitianEJA._denormalized_basis(n,QQ)
1685 sage: all( M.is_symmetric() for M in B )
1689 Q
= QuaternionAlgebra(QQ
,-1,-1)
1692 # This is like the symmetric case, but we need to be careful:
1694 # * We want conjugate-symmetry, not just symmetry.
1695 # * The diagonal will (as a result) be real.
1699 for j
in xrange(i
+1):
1700 Eij
= matrix(Q
, n
, lambda k
,l
: k
==i
and l
==j
)
1702 Sij
= cls
.real_embed(Eij
)
1705 # The second, third, and fourth ones have a minus
1706 # because they're conjugated.
1707 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
1709 Sij_I
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
1711 Sij_J
= cls
.real_embed(J
*Eij
- J
*Eij
.transpose())
1713 Sij_K
= cls
.real_embed(K
*Eij
- K
*Eij
.transpose())
1716 # Since we embedded these, we can drop back to the "field" that we
1717 # started with instead of the quaternion algebra "Q".
1718 return ( s
.change_ring(field
) for s
in S
)
1721 def __init__(self
, n
, field
=QQ
, **kwargs
):
1722 basis
= self
._denormalized
_basis
(n
,field
)
1723 super(QuaternionHermitianEJA
,self
).__init
__(field
, basis
, n
, **kwargs
)
1726 class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra
, KnownRankEJA
):
1728 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
1729 with the usual inner product and jordan product ``x*y =
1730 (<x_bar,y_bar>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
1735 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1739 This multiplication table can be verified by hand::
1741 sage: J = JordanSpinEJA(4)
1742 sage: e0,e1,e2,e3 = J.gens()
1758 We can change the generator prefix::
1760 sage: JordanSpinEJA(2, prefix='B').gens()
1764 def __init__(self
, n
, field
=QQ
, **kwargs
):
1765 V
= VectorSpace(field
, n
)
1766 mult_table
= [[V
.zero() for j
in xrange(n
)] for i
in xrange(n
)]
1776 z0
= x
.inner_product(y
)
1777 zbar
= y0
*xbar
+ x0
*ybar
1778 z
= V([z0
] + zbar
.list())
1779 mult_table
[i
][j
] = z
1781 # The rank of the spin algebra is two, unless we're in a
1782 # one-dimensional ambient space (because the rank is bounded by
1783 # the ambient dimension).
1784 fdeja
= super(JordanSpinEJA
, self
)
1785 return fdeja
.__init
__(field
, mult_table
, rank
=min(n
,2), **kwargs
)
1787 def inner_product(self
, x
, y
):
1789 Faster to reimplement than to use natural representations.
1793 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1797 Ensure that this is the usual inner product for the algebras
1800 sage: set_random_seed()
1801 sage: J = JordanSpinEJA.random_instance()
1802 sage: x,y = J.random_elements(2)
1803 sage: X = x.natural_representation()
1804 sage: Y = y.natural_representation()
1805 sage: x.inner_product(y) == J.natural_inner_product(X,Y)
1809 return x
.to_vector().inner_product(y
.to_vector())