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
.all
import (ZZ
, QQ
, RR
, RLF
, CLF
,
22 from mjo
.eja
.eja_element
import FiniteDimensionalEuclideanJordanAlgebraElement
23 from mjo
.eja
.eja_utils
import _mat2vec
25 class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule
):
26 # This is an ugly hack needed to prevent the category framework
27 # from implementing a coercion from our base ring (e.g. the
28 # rationals) into the algebra. First of all -- such a coercion is
29 # nonsense to begin with. But more importantly, it tries to do so
30 # in the category of rings, and since our algebras aren't
31 # associative they generally won't be rings.
32 _no_generic_basering_coercion
= True
45 sage: from mjo.eja.eja_algebra import (JordanSpinEJA, random_eja)
49 By definition, Jordan multiplication commutes::
51 sage: set_random_seed()
52 sage: J = random_eja()
53 sage: x,y = J.random_elements(2)
59 The ``field`` we're given must be real::
61 sage: JordanSpinEJA(2,QQbar)
62 Traceback (most recent call last):
64 ValueError: field is not real
68 if not field
.is_subring(RR
):
69 # Note: this does return true for the real algebraic
70 # field, and any quadratic field where we've specified
72 raise ValueError('field is not real')
75 self
._natural
_basis
= natural_basis
78 category
= MagmaticAlgebras(field
).FiniteDimensional()
79 category
= category
.WithBasis().Unital()
81 fda
= super(FiniteDimensionalEuclideanJordanAlgebra
, self
)
83 range(len(mult_table
)),
86 self
.print_options(bracket
='')
88 # The multiplication table we're given is necessarily in terms
89 # of vectors, because we don't have an algebra yet for
90 # anything to be an element of. However, it's faster in the
91 # long run to have the multiplication table be in terms of
92 # algebra elements. We do this after calling the superclass
93 # constructor so that from_vector() knows what to do.
94 self
._multiplication
_table
= [ map(lambda x
: self
.from_vector(x
), ls
)
95 for ls
in mult_table
]
98 def _element_constructor_(self
, elt
):
100 Construct an element of this algebra from its natural
103 This gets called only after the parent element _call_ method
104 fails to find a coercion for the argument.
108 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
109 ....: RealCartesianProductEJA,
110 ....: RealSymmetricEJA)
114 The identity in `S^n` is converted to the identity in the EJA::
116 sage: J = RealSymmetricEJA(3)
117 sage: I = matrix.identity(QQ,3)
118 sage: J(I) == J.one()
121 This skew-symmetric matrix can't be represented in the EJA::
123 sage: J = RealSymmetricEJA(3)
124 sage: A = matrix(QQ,3, lambda i,j: i-j)
126 Traceback (most recent call last):
128 ArithmeticError: vector is not in free module
132 Ensure that we can convert any element of the two non-matrix
133 simple algebras (whose natural representations are their usual
134 vector representations) back and forth faithfully::
136 sage: set_random_seed()
137 sage: J = RealCartesianProductEJA.random_instance()
138 sage: x = J.random_element()
139 sage: J(x.to_vector().column()) == x
141 sage: J = JordanSpinEJA.random_instance()
142 sage: x = J.random_element()
143 sage: J(x.to_vector().column()) == x
148 # The superclass implementation of random_element()
149 # needs to be able to coerce "0" into the algebra.
152 natural_basis
= self
.natural_basis()
153 basis_space
= natural_basis
[0].matrix_space()
154 if elt
not in basis_space
:
155 raise ValueError("not a naturally-represented algebra element")
157 # Thanks for nothing! Matrix spaces aren't vector spaces in
158 # Sage, so we have to figure out its natural-basis coordinates
159 # ourselves. We use the basis space's ring instead of the
160 # element's ring because the basis space might be an algebraic
161 # closure whereas the base ring of the 3-by-3 identity matrix
162 # could be QQ instead of QQbar.
163 V
= VectorSpace(basis_space
.base_ring(), elt
.nrows()*elt
.ncols())
164 W
= V
.span_of_basis( _mat2vec(s
) for s
in natural_basis
)
165 coords
= W
.coordinate_vector(_mat2vec(elt
))
166 return self
.from_vector(coords
)
171 Return a string representation of ``self``.
175 sage: from mjo.eja.eja_algebra import JordanSpinEJA
179 Ensure that it says what we think it says::
181 sage: JordanSpinEJA(2, field=QQ)
182 Euclidean Jordan algebra of dimension 2 over Rational Field
183 sage: JordanSpinEJA(3, field=RDF)
184 Euclidean Jordan algebra of dimension 3 over Real Double Field
187 fmt
= "Euclidean Jordan algebra of dimension {} over {}"
188 return fmt
.format(self
.dimension(), self
.base_ring())
190 def product_on_basis(self
, i
, j
):
191 return self
._multiplication
_table
[i
][j
]
193 def _a_regular_element(self
):
195 Guess a regular element. Needed to compute the basis for our
196 characteristic polynomial coefficients.
200 sage: from mjo.eja.eja_algebra import random_eja
204 Ensure that this hacky method succeeds for every algebra that we
205 know how to construct::
207 sage: set_random_seed()
208 sage: J = random_eja()
209 sage: J._a_regular_element().is_regular()
214 z
= self
.sum( (i
+1)*gs
[i
] for i
in range(len(gs
)) )
215 if not z
.is_regular():
216 raise ValueError("don't know a regular element")
221 def _charpoly_basis_space(self
):
223 Return the vector space spanned by the basis used in our
224 characteristic polynomial coefficients. This is used not only to
225 compute those coefficients, but also any time we need to
226 evaluate the coefficients (like when we compute the trace or
229 z
= self
._a
_regular
_element
()
230 # Don't use the parent vector space directly here in case this
231 # happens to be a subalgebra. In that case, we would be e.g.
232 # two-dimensional but span_of_basis() would expect three
234 V
= VectorSpace(self
.base_ring(), self
.vector_space().dimension())
235 basis
= [ (z
**k
).to_vector() for k
in range(self
.rank()) ]
236 V1
= V
.span_of_basis( basis
)
237 b
= (V1
.basis() + V1
.complement().basis())
238 return V
.span_of_basis(b
)
243 def _charpoly_coeff(self
, i
):
245 Return the coefficient polynomial "a_{i}" of this algebra's
246 general characteristic polynomial.
248 Having this be a separate cached method lets us compute and
249 store the trace/determinant (a_{r-1} and a_{0} respectively)
250 separate from the entire characteristic polynomial.
252 (A_of_x
, x
, xr
, detA
) = self
._charpoly
_matrix
_system
()
253 R
= A_of_x
.base_ring()
255 # Guaranteed by theory
258 # Danger: the in-place modification is done for performance
259 # reasons (reconstructing a matrix with huge polynomial
260 # entries is slow), but I don't know how cached_method works,
261 # so it's highly possible that we're modifying some global
262 # list variable by reference, here. In other words, you
263 # probably shouldn't call this method twice on the same
264 # algebra, at the same time, in two threads
265 Ai_orig
= A_of_x
.column(i
)
266 A_of_x
.set_column(i
,xr
)
267 numerator
= A_of_x
.det()
268 A_of_x
.set_column(i
,Ai_orig
)
270 # We're relying on the theory here to ensure that each a_i is
271 # indeed back in R, and the added negative signs are to make
272 # the whole charpoly expression sum to zero.
273 return R(-numerator
/detA
)
277 def _charpoly_matrix_system(self
):
279 Compute the matrix whose entries A_ij are polynomials in
280 X1,...,XN, the vector ``x`` of variables X1,...,XN, the vector
281 corresponding to `x^r` and the determinent of the matrix A =
282 [A_ij]. In other words, all of the fixed (cachable) data needed
283 to compute the coefficients of the characteristic polynomial.
288 # Turn my vector space into a module so that "vectors" can
289 # have multivatiate polynomial entries.
290 names
= tuple('X' + str(i
) for i
in range(1,n
+1))
291 R
= PolynomialRing(self
.base_ring(), names
)
293 # Using change_ring() on the parent's vector space doesn't work
294 # here because, in a subalgebra, that vector space has a basis
295 # and change_ring() tries to bring the basis along with it. And
296 # that doesn't work unless the new ring is a PID, which it usually
300 # Now let x = (X1,X2,...,Xn) be the vector whose entries are
304 # And figure out the "left multiplication by x" matrix in
307 monomial_matrices
= [ self
.monomial(i
).operator().matrix()
308 for i
in range(n
) ] # don't recompute these!
310 ek
= self
.monomial(k
).to_vector()
312 sum( x
[i
]*(monomial_matrices
[i
]*ek
)
313 for i
in range(n
) ) )
314 Lx
= matrix
.column(R
, lmbx_cols
)
316 # Now we can compute powers of x "symbolically"
317 x_powers
= [self
.one().to_vector(), x
]
318 for d
in range(2, r
+1):
319 x_powers
.append( Lx
*(x_powers
[-1]) )
321 idmat
= matrix
.identity(R
, n
)
323 W
= self
._charpoly
_basis
_space
()
324 W
= W
.change_ring(R
.fraction_field())
326 # Starting with the standard coordinates x = (X1,X2,...,Xn)
327 # and then converting the entries to W-coordinates allows us
328 # to pass in the standard coordinates to the charpoly and get
329 # back the right answer. Specifically, with x = (X1,X2,...,Xn),
332 # W.coordinates(x^2) eval'd at (standard z-coords)
336 # W-coords of (standard coords of x^2 eval'd at std-coords of z)
338 # We want the middle equivalent thing in our matrix, but use
339 # the first equivalent thing instead so that we can pass in
340 # standard coordinates.
341 x_powers
= [ W
.coordinate_vector(xp
) for xp
in x_powers
]
342 l2
= [idmat
.column(k
-1) for k
in range(r
+1, n
+1)]
343 A_of_x
= matrix
.column(R
, n
, (x_powers
[:r
] + l2
))
344 return (A_of_x
, x
, x_powers
[r
], A_of_x
.det())
348 def characteristic_polynomial(self
):
350 Return a characteristic polynomial that works for all elements
353 The resulting polynomial has `n+1` variables, where `n` is the
354 dimension of this algebra. The first `n` variables correspond to
355 the coordinates of an algebra element: when evaluated at the
356 coordinates of an algebra element with respect to a certain
357 basis, the result is a univariate polynomial (in the one
358 remaining variable ``t``), namely the characteristic polynomial
363 sage: from mjo.eja.eja_algebra import JordanSpinEJA
367 The characteristic polynomial in the spin algebra is given in
368 Alizadeh, Example 11.11::
370 sage: J = JordanSpinEJA(3)
371 sage: p = J.characteristic_polynomial(); p
372 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
373 sage: xvec = J.one().to_vector()
381 # The list of coefficient polynomials a_1, a_2, ..., a_n.
382 a
= [ self
._charpoly
_coeff
(i
) for i
in range(n
) ]
384 # We go to a bit of trouble here to reorder the
385 # indeterminates, so that it's easier to evaluate the
386 # characteristic polynomial at x's coordinates and get back
387 # something in terms of t, which is what we want.
389 S
= PolynomialRing(self
.base_ring(),'t')
391 S
= PolynomialRing(S
, R
.variable_names())
394 # Note: all entries past the rth should be zero. The
395 # coefficient of the highest power (x^r) is 1, but it doesn't
396 # appear in the solution vector which contains coefficients
397 # for the other powers (to make them sum to x^r).
399 a
[r
] = 1 # corresponds to x^r
401 # When the rank is equal to the dimension, trying to
402 # assign a[r] goes out-of-bounds.
403 a
.append(1) # corresponds to x^r
405 return sum( a
[k
]*(t
**k
) for k
in xrange(len(a
)) )
408 def inner_product(self
, x
, y
):
410 The inner product associated with this Euclidean Jordan algebra.
412 Defaults to the trace inner product, but can be overridden by
413 subclasses if they are sure that the necessary properties are
418 sage: from mjo.eja.eja_algebra import random_eja
422 Our inner product is "associative," which means the following for
423 a symmetric bilinear form::
425 sage: set_random_seed()
426 sage: J = random_eja()
427 sage: x,y,z = J.random_elements(3)
428 sage: (x*y).inner_product(z) == y.inner_product(x*z)
432 X
= x
.natural_representation()
433 Y
= y
.natural_representation()
434 return self
.natural_inner_product(X
,Y
)
437 def is_trivial(self
):
439 Return whether or not this algebra is trivial.
441 A trivial algebra contains only the zero element.
445 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
449 sage: J = ComplexHermitianEJA(3)
452 sage: A = J.zero().subalgebra_generated_by()
457 return self
.dimension() == 0
460 def multiplication_table(self
):
462 Return a visual representation of this algebra's multiplication
463 table (on basis elements).
467 sage: from mjo.eja.eja_algebra import JordanSpinEJA
471 sage: J = JordanSpinEJA(4)
472 sage: J.multiplication_table()
473 +----++----+----+----+----+
474 | * || e0 | e1 | e2 | e3 |
475 +====++====+====+====+====+
476 | e0 || e0 | e1 | e2 | e3 |
477 +----++----+----+----+----+
478 | e1 || e1 | e0 | 0 | 0 |
479 +----++----+----+----+----+
480 | e2 || e2 | 0 | e0 | 0 |
481 +----++----+----+----+----+
482 | e3 || e3 | 0 | 0 | e0 |
483 +----++----+----+----+----+
486 M
= list(self
._multiplication
_table
) # copy
487 for i
in xrange(len(M
)):
488 # M had better be "square"
489 M
[i
] = [self
.monomial(i
)] + M
[i
]
490 M
= [["*"] + list(self
.gens())] + M
491 return table(M
, header_row
=True, header_column
=True, frame
=True)
494 def natural_basis(self
):
496 Return a more-natural representation of this algebra's basis.
498 Every finite-dimensional Euclidean Jordan Algebra is a direct
499 sum of five simple algebras, four of which comprise Hermitian
500 matrices. This method returns the original "natural" basis
501 for our underlying vector space. (Typically, the natural basis
502 is used to construct the multiplication table in the first place.)
504 Note that this will always return a matrix. The standard basis
505 in `R^n` will be returned as `n`-by-`1` column matrices.
509 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
510 ....: RealSymmetricEJA)
514 sage: J = RealSymmetricEJA(2)
516 Finite family {0: e0, 1: e1, 2: e2}
517 sage: J.natural_basis()
519 [1 0] [ 0 1/2*sqrt2] [0 0]
520 [0 0], [1/2*sqrt2 0], [0 1]
525 sage: J = JordanSpinEJA(2)
527 Finite family {0: e0, 1: e1}
528 sage: J.natural_basis()
535 if self
._natural
_basis
is None:
536 M
= self
.natural_basis_space()
537 return tuple( M(b
.to_vector()) for b
in self
.basis() )
539 return self
._natural
_basis
542 def natural_basis_space(self
):
544 Return the matrix space in which this algebra's natural basis
547 if self
._natural
_basis
is None or len(self
._natural
_basis
) == 0:
548 return MatrixSpace(self
.base_ring(), self
.dimension(), 1)
550 return self
._natural
_basis
[0].matrix_space()
554 def natural_inner_product(X
,Y
):
556 Compute the inner product of two naturally-represented elements.
558 For example in the real symmetric matrix EJA, this will compute
559 the trace inner-product of two n-by-n symmetric matrices. The
560 default should work for the real cartesian product EJA, the
561 Jordan spin EJA, and the real symmetric matrices. The others
562 will have to be overridden.
564 return (X
.conjugate_transpose()*Y
).trace()
570 Return the unit element of this algebra.
574 sage: from mjo.eja.eja_algebra import (RealCartesianProductEJA,
579 sage: J = RealCartesianProductEJA(5)
581 e0 + e1 + e2 + e3 + e4
585 The identity element acts like the identity::
587 sage: set_random_seed()
588 sage: J = random_eja()
589 sage: x = J.random_element()
590 sage: J.one()*x == x and x*J.one() == x
593 The matrix of the unit element's operator is the identity::
595 sage: set_random_seed()
596 sage: J = random_eja()
597 sage: actual = J.one().operator().matrix()
598 sage: expected = matrix.identity(J.base_ring(), J.dimension())
599 sage: actual == expected
603 # We can brute-force compute the matrices of the operators
604 # that correspond to the basis elements of this algebra.
605 # If some linear combination of those basis elements is the
606 # algebra identity, then the same linear combination of
607 # their matrices has to be the identity matrix.
609 # Of course, matrices aren't vectors in sage, so we have to
610 # appeal to the "long vectors" isometry.
611 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
613 # Now we use basis linear algebra to find the coefficients,
614 # of the matrices-as-vectors-linear-combination, which should
615 # work for the original algebra basis too.
616 A
= matrix
.column(self
.base_ring(), oper_vecs
)
618 # We used the isometry on the left-hand side already, but we
619 # still need to do it for the right-hand side. Recall that we
620 # wanted something that summed to the identity matrix.
621 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
623 # Now if there's an identity element in the algebra, this should work.
624 coeffs
= A
.solve_right(b
)
625 return self
.linear_combination(zip(self
.gens(), coeffs
))
628 def random_elements(self
, count
):
630 Return ``count`` random elements as a tuple.
634 sage: from mjo.eja.eja_algebra import JordanSpinEJA
638 sage: J = JordanSpinEJA(3)
639 sage: x,y,z = J.random_elements(3)
640 sage: all( [ x in J, y in J, z in J ])
642 sage: len( J.random_elements(10) ) == 10
646 return tuple( self
.random_element() for idx
in xrange(count
) )
651 Return the rank of this EJA.
655 The author knows of no algorithm to compute the rank of an EJA
656 where only the multiplication table is known. In lieu of one, we
657 require the rank to be specified when the algebra is created,
658 and simply pass along that number here.
662 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
663 ....: RealSymmetricEJA,
664 ....: ComplexHermitianEJA,
665 ....: QuaternionHermitianEJA,
670 The rank of the Jordan spin algebra is always two::
672 sage: JordanSpinEJA(2).rank()
674 sage: JordanSpinEJA(3).rank()
676 sage: JordanSpinEJA(4).rank()
679 The rank of the `n`-by-`n` Hermitian real, complex, or
680 quaternion matrices is `n`::
682 sage: RealSymmetricEJA(4).rank()
684 sage: ComplexHermitianEJA(3).rank()
686 sage: QuaternionHermitianEJA(2).rank()
691 Ensure that every EJA that we know how to construct has a
692 positive integer rank::
694 sage: set_random_seed()
695 sage: r = random_eja().rank()
696 sage: r in ZZ and r > 0
703 def vector_space(self
):
705 Return the vector space that underlies this algebra.
709 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
713 sage: J = RealSymmetricEJA(2)
714 sage: J.vector_space()
715 Vector space of dimension 3 over...
718 return self
.zero().to_vector().parent().ambient_vector_space()
721 Element
= FiniteDimensionalEuclideanJordanAlgebraElement
724 class KnownRankEJA(object):
726 A class for algebras that we actually know we can construct. The
727 main issue is that, for most of our methods to make sense, we need
728 to know the rank of our algebra. Thus we can't simply generate a
729 "random" algebra, or even check that a given basis and product
730 satisfy the axioms; because even if everything looks OK, we wouldn't
731 know the rank we need to actuallty build the thing.
733 Not really a subclass of FDEJA because doing that causes method
734 resolution errors, e.g.
736 TypeError: Error when calling the metaclass bases
737 Cannot create a consistent method resolution
738 order (MRO) for bases FiniteDimensionalEuclideanJordanAlgebra,
743 def _max_test_case_size():
745 Return an integer "size" that is an upper bound on the size of
746 this algebra when it is used in a random test
747 case. Unfortunately, the term "size" is quite vague -- when
748 dealing with `R^n` under either the Hadamard or Jordan spin
749 product, the "size" refers to the dimension `n`. When dealing
750 with a matrix algebra (real symmetric or complex/quaternion
751 Hermitian), it refers to the size of the matrix, which is
752 far less than the dimension of the underlying vector space.
754 We default to five in this class, which is safe in `R^n`. The
755 matrix algebra subclasses (or any class where the "size" is
756 interpreted to be far less than the dimension) should override
757 with a smaller number.
762 def random_instance(cls
, field
=QQ
, **kwargs
):
764 Return a random instance of this type of algebra.
766 Beware, this will crash for "most instances" because the
767 constructor below looks wrong.
769 n
= ZZ
.random_element(cls
._max
_test
_case
_size
()) + 1
770 return cls(n
, field
, **kwargs
)
773 class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra
,
776 Return the Euclidean Jordan Algebra corresponding to the set
777 `R^n` under the Hadamard product.
779 Note: this is nothing more than the Cartesian product of ``n``
780 copies of the spin algebra. Once Cartesian product algebras
781 are implemented, this can go.
785 sage: from mjo.eja.eja_algebra import RealCartesianProductEJA
789 This multiplication table can be verified by hand::
791 sage: J = RealCartesianProductEJA(3)
792 sage: e0,e1,e2 = J.gens()
808 We can change the generator prefix::
810 sage: RealCartesianProductEJA(3, prefix='r').gens()
814 def __init__(self
, n
, field
=QQ
, **kwargs
):
815 V
= VectorSpace(field
, n
)
816 mult_table
= [ [ V
.gen(i
)*(i
== j
) for j
in xrange(n
) ]
819 fdeja
= super(RealCartesianProductEJA
, self
)
820 return fdeja
.__init
__(field
, mult_table
, rank
=n
, **kwargs
)
822 def inner_product(self
, x
, y
):
824 Faster to reimplement than to use natural representations.
828 sage: from mjo.eja.eja_algebra import RealCartesianProductEJA
832 Ensure that this is the usual inner product for the algebras
835 sage: set_random_seed()
836 sage: J = RealCartesianProductEJA.random_instance()
837 sage: x,y = J.random_elements(2)
838 sage: X = x.natural_representation()
839 sage: Y = y.natural_representation()
840 sage: x.inner_product(y) == J.natural_inner_product(X,Y)
844 return x
.to_vector().inner_product(y
.to_vector())
847 def random_eja(field
=QQ
):
849 Return a "random" finite-dimensional Euclidean Jordan Algebra.
853 For now, we choose a random natural number ``n`` (greater than zero)
854 and then give you back one of the following:
856 * The cartesian product of the rational numbers ``n`` times; this is
857 ``QQ^n`` with the Hadamard product.
859 * The Jordan spin algebra on ``QQ^n``.
861 * The ``n``-by-``n`` rational symmetric matrices with the symmetric
864 * The ``n``-by-``n`` complex-rational Hermitian matrices embedded
865 in the space of ``2n``-by-``2n`` real symmetric matrices.
867 * The ``n``-by-``n`` quaternion-rational Hermitian matrices embedded
868 in the space of ``4n``-by-``4n`` real symmetric matrices.
870 Later this might be extended to return Cartesian products of the
875 sage: from mjo.eja.eja_algebra import random_eja
880 Euclidean Jordan algebra of dimension...
883 classname
= choice(KnownRankEJA
.__subclasses
__())
884 return classname
.random_instance(field
=field
)
891 class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra
):
893 def _max_test_case_size():
894 # Play it safe, since this will be squared and the underlying
895 # field can have dimension 4 (quaternions) too.
898 def __init__(self
, field
, basis
, rank
, normalize_basis
=True, **kwargs
):
900 Compared to the superclass constructor, we take a basis instead of
901 a multiplication table because the latter can be computed in terms
902 of the former when the product is known (like it is here).
904 # Used in this class's fast _charpoly_coeff() override.
905 self
._basis
_normalizers
= None
907 # We're going to loop through this a few times, so now's a good
908 # time to ensure that it isn't a generator expression.
911 if rank
> 1 and normalize_basis
:
912 # We'll need sqrt(2) to normalize the basis, and this
913 # winds up in the multiplication table, so the whole
914 # algebra needs to be over the field extension.
915 R
= PolynomialRing(field
, 'z')
918 if p
.is_irreducible():
919 field
= field
.extension(p
, 'sqrt2', embedding
=RLF(2).sqrt())
920 basis
= tuple( s
.change_ring(field
) for s
in basis
)
921 self
._basis
_normalizers
= tuple(
922 ~
(self
.natural_inner_product(s
,s
).sqrt()) for s
in basis
)
923 basis
= tuple(s
*c
for (s
,c
) in izip(basis
,self
._basis
_normalizers
))
925 Qs
= self
.multiplication_table_from_matrix_basis(basis
)
927 fdeja
= super(MatrixEuclideanJordanAlgebra
, self
)
928 return fdeja
.__init
__(field
,
936 def _charpoly_coeff(self
, i
):
938 Override the parent method with something that tries to compute
939 over a faster (non-extension) field.
941 if self
._basis
_normalizers
is None:
942 # We didn't normalize, so assume that the basis we started
943 # with had entries in a nice field.
944 return super(MatrixEuclideanJordanAlgebra
, self
)._charpoly
_coeff
(i
)
946 basis
= ( (b
/n
) for (b
,n
) in izip(self
.natural_basis(),
947 self
._basis
_normalizers
) )
949 # Do this over the rationals and convert back at the end.
950 J
= MatrixEuclideanJordanAlgebra(QQ
,
953 normalize_basis
=False)
954 (_
,x
,_
,_
) = J
._charpoly
_matrix
_system
()
955 p
= J
._charpoly
_coeff
(i
)
956 # p might be missing some vars, have to substitute "optionally"
957 pairs
= izip(x
.base_ring().gens(), self
._basis
_normalizers
)
958 substitutions
= { v: v*c for (v,c) in pairs }
959 result
= p
.subs(substitutions
)
961 # The result of "subs" can be either a coefficient-ring
962 # element or a polynomial. Gotta handle both cases.
964 return self
.base_ring()(result
)
966 return result
.change_ring(self
.base_ring())
970 def multiplication_table_from_matrix_basis(basis
):
972 At least three of the five simple Euclidean Jordan algebras have the
973 symmetric multiplication (A,B) |-> (AB + BA)/2, where the
974 multiplication on the right is matrix multiplication. Given a basis
975 for the underlying matrix space, this function returns a
976 multiplication table (obtained by looping through the basis
977 elements) for an algebra of those matrices.
979 # In S^2, for example, we nominally have four coordinates even
980 # though the space is of dimension three only. The vector space V
981 # is supposed to hold the entire long vector, and the subspace W
982 # of V will be spanned by the vectors that arise from symmetric
983 # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
984 field
= basis
[0].base_ring()
985 dimension
= basis
[0].nrows()
987 V
= VectorSpace(field
, dimension
**2)
988 W
= V
.span_of_basis( _mat2vec(s
) for s
in basis
)
990 mult_table
= [[W
.zero() for j
in xrange(n
)] for i
in xrange(n
)]
993 mat_entry
= (basis
[i
]*basis
[j
] + basis
[j
]*basis
[i
])/2
994 mult_table
[i
][j
] = W
.coordinate_vector(_mat2vec(mat_entry
))
1002 Embed the matrix ``M`` into a space of real matrices.
1004 The matrix ``M`` can have entries in any field at the moment:
1005 the real numbers, complex numbers, or quaternions. And although
1006 they are not a field, we can probably support octonions at some
1007 point, too. This function returns a real matrix that "acts like"
1008 the original with respect to matrix multiplication; i.e.
1010 real_embed(M*N) = real_embed(M)*real_embed(N)
1013 raise NotImplementedError
1017 def real_unembed(M
):
1019 The inverse of :meth:`real_embed`.
1021 raise NotImplementedError
1025 def natural_inner_product(cls
,X
,Y
):
1026 Xu
= cls
.real_unembed(X
)
1027 Yu
= cls
.real_unembed(Y
)
1028 tr
= (Xu
*Yu
).trace()
1031 # It's real already.
1034 # Otherwise, try the thing that works for complex numbers; and
1035 # if that doesn't work, the thing that works for quaternions.
1037 return tr
.vector()[0] # real part, imag part is index 1
1038 except AttributeError:
1039 # A quaternions doesn't have a vector() method, but does
1040 # have coefficient_tuple() method that returns the
1041 # coefficients of 1, i, j, and k -- in that order.
1042 return tr
.coefficient_tuple()[0]
1045 class RealMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1049 The identity function, for embedding real matrices into real
1055 def real_unembed(M
):
1057 The identity function, for unembedding real matrices from real
1063 class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra
, KnownRankEJA
):
1065 The rank-n simple EJA consisting of real symmetric n-by-n
1066 matrices, the usual symmetric Jordan product, and the trace inner
1067 product. It has dimension `(n^2 + n)/2` over the reals.
1071 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1075 sage: J = RealSymmetricEJA(2)
1076 sage: e0, e1, e2 = J.gens()
1084 In theory, our "field" can be any subfield of the reals::
1086 sage: RealSymmetricEJA(2, AA)
1087 Euclidean Jordan algebra of dimension 3 over Algebraic Real Field
1088 sage: RealSymmetricEJA(2, RR)
1089 Euclidean Jordan algebra of dimension 3 over Real Field with
1090 53 bits of precision
1094 The dimension of this algebra is `(n^2 + n) / 2`::
1096 sage: set_random_seed()
1097 sage: n_max = RealSymmetricEJA._max_test_case_size()
1098 sage: n = ZZ.random_element(1, n_max)
1099 sage: J = RealSymmetricEJA(n)
1100 sage: J.dimension() == (n^2 + n)/2
1103 The Jordan multiplication is what we think it is::
1105 sage: set_random_seed()
1106 sage: J = RealSymmetricEJA.random_instance()
1107 sage: x,y = J.random_elements(2)
1108 sage: actual = (x*y).natural_representation()
1109 sage: X = x.natural_representation()
1110 sage: Y = y.natural_representation()
1111 sage: expected = (X*Y + Y*X)/2
1112 sage: actual == expected
1114 sage: J(expected) == x*y
1117 We can change the generator prefix::
1119 sage: RealSymmetricEJA(3, prefix='q').gens()
1120 (q0, q1, q2, q3, q4, q5)
1122 Our natural basis is normalized with respect to the natural inner
1123 product unless we specify otherwise::
1125 sage: set_random_seed()
1126 sage: J = RealSymmetricEJA.random_instance()
1127 sage: all( b.norm() == 1 for b in J.gens() )
1130 Since our natural basis is normalized with respect to the natural
1131 inner product, and since we know that this algebra is an EJA, any
1132 left-multiplication operator's matrix will be symmetric because
1133 natural->EJA basis representation is an isometry and within the EJA
1134 the operator is self-adjoint by the Jordan axiom::
1136 sage: set_random_seed()
1137 sage: x = RealSymmetricEJA.random_instance().random_element()
1138 sage: x.operator().matrix().is_symmetric()
1143 def _denormalized_basis(cls
, n
, field
):
1145 Return a basis for the space of real symmetric n-by-n matrices.
1149 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1153 sage: set_random_seed()
1154 sage: n = ZZ.random_element(1,5)
1155 sage: B = RealSymmetricEJA._denormalized_basis(n,QQ)
1156 sage: all( M.is_symmetric() for M in B)
1160 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1164 for j
in xrange(i
+1):
1165 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
1169 Sij
= Eij
+ Eij
.transpose()
1175 def _max_test_case_size():
1176 return 4 # Dimension 10
1179 def __init__(self
, n
, field
=QQ
, **kwargs
):
1180 basis
= self
._denormalized
_basis
(n
, field
)
1181 super(RealSymmetricEJA
, self
).__init
__(field
, basis
, n
, **kwargs
)
1184 class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1188 Embed the n-by-n complex matrix ``M`` into the space of real
1189 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1190 bi` to the block matrix ``[[a,b],[-b,a]]``.
1194 sage: from mjo.eja.eja_algebra import \
1195 ....: ComplexMatrixEuclideanJordanAlgebra
1199 sage: F = QuadraticField(-1, 'i')
1200 sage: x1 = F(4 - 2*i)
1201 sage: x2 = F(1 + 2*i)
1204 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1205 sage: ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1214 Embedding is a homomorphism (isomorphism, in fact)::
1216 sage: set_random_seed()
1217 sage: n_max = ComplexMatrixEuclideanJordanAlgebra._max_test_case_size()
1218 sage: n = ZZ.random_element(n_max)
1219 sage: F = QuadraticField(-1, 'i')
1220 sage: X = random_matrix(F, n)
1221 sage: Y = random_matrix(F, n)
1222 sage: Xe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X)
1223 sage: Ye = ComplexMatrixEuclideanJordanAlgebra.real_embed(Y)
1224 sage: XYe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1231 raise ValueError("the matrix 'M' must be square")
1233 # We don't need any adjoined elements...
1234 field
= M
.base_ring().base_ring()
1238 a
= z
.list()[0] # real part, I guess
1239 b
= z
.list()[1] # imag part, I guess
1240 blocks
.append(matrix(field
, 2, [[a
,b
],[-b
,a
]]))
1242 return matrix
.block(field
, n
, blocks
)
1246 def real_unembed(M
):
1248 The inverse of _embed_complex_matrix().
1252 sage: from mjo.eja.eja_algebra import \
1253 ....: ComplexMatrixEuclideanJordanAlgebra
1257 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1258 ....: [-2, 1, -4, 3],
1259 ....: [ 9, 10, 11, 12],
1260 ....: [-10, 9, -12, 11] ])
1261 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(A)
1263 [ 10*i + 9 12*i + 11]
1267 Unembedding is the inverse of embedding::
1269 sage: set_random_seed()
1270 sage: F = QuadraticField(-1, 'i')
1271 sage: M = random_matrix(F, 3)
1272 sage: Me = ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1273 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1279 raise ValueError("the matrix 'M' must be square")
1280 if not n
.mod(2).is_zero():
1281 raise ValueError("the matrix 'M' must be a complex embedding")
1283 # If "M" was normalized, its base ring might have roots
1284 # adjoined and they can stick around after unembedding.
1285 field
= M
.base_ring()
1286 R
= PolynomialRing(field
, 'z')
1288 F
= field
.extension(z
**2 + 1, 'i', embedding
=CLF(-1).sqrt())
1291 # Go top-left to bottom-right (reading order), converting every
1292 # 2-by-2 block we see to a single complex element.
1294 for k
in xrange(n
/2):
1295 for j
in xrange(n
/2):
1296 submat
= M
[2*k
:2*k
+2,2*j
:2*j
+2]
1297 if submat
[0,0] != submat
[1,1]:
1298 raise ValueError('bad on-diagonal submatrix')
1299 if submat
[0,1] != -submat
[1,0]:
1300 raise ValueError('bad off-diagonal submatrix')
1301 z
= submat
[0,0] + submat
[0,1]*i
1304 return matrix(F
, n
/2, elements
)
1308 def natural_inner_product(cls
,X
,Y
):
1310 Compute a natural inner product in this algebra directly from
1315 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1319 This gives the same answer as the slow, default method implemented
1320 in :class:`MatrixEuclideanJordanAlgebra`::
1322 sage: set_random_seed()
1323 sage: J = ComplexHermitianEJA.random_instance()
1324 sage: x,y = J.random_elements(2)
1325 sage: Xe = x.natural_representation()
1326 sage: Ye = y.natural_representation()
1327 sage: X = ComplexHermitianEJA.real_unembed(Xe)
1328 sage: Y = ComplexHermitianEJA.real_unembed(Ye)
1329 sage: expected = (X*Y).trace().vector()[0]
1330 sage: actual = ComplexHermitianEJA.natural_inner_product(Xe,Ye)
1331 sage: actual == expected
1335 return RealMatrixEuclideanJordanAlgebra
.natural_inner_product(X
,Y
)/2
1338 class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra
, KnownRankEJA
):
1340 The rank-n simple EJA consisting of complex Hermitian n-by-n
1341 matrices over the real numbers, the usual symmetric Jordan product,
1342 and the real-part-of-trace inner product. It has dimension `n^2` over
1347 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1351 In theory, our "field" can be any subfield of the reals::
1353 sage: ComplexHermitianEJA(2, AA)
1354 Euclidean Jordan algebra of dimension 4 over Algebraic Real Field
1355 sage: ComplexHermitianEJA(2, RR)
1356 Euclidean Jordan algebra of dimension 4 over Real Field with
1357 53 bits of precision
1361 The dimension of this algebra is `n^2`::
1363 sage: set_random_seed()
1364 sage: n_max = ComplexHermitianEJA._max_test_case_size()
1365 sage: n = ZZ.random_element(1, n_max)
1366 sage: J = ComplexHermitianEJA(n)
1367 sage: J.dimension() == n^2
1370 The Jordan multiplication is what we think it is::
1372 sage: set_random_seed()
1373 sage: J = ComplexHermitianEJA.random_instance()
1374 sage: x,y = J.random_elements(2)
1375 sage: actual = (x*y).natural_representation()
1376 sage: X = x.natural_representation()
1377 sage: Y = y.natural_representation()
1378 sage: expected = (X*Y + Y*X)/2
1379 sage: actual == expected
1381 sage: J(expected) == x*y
1384 We can change the generator prefix::
1386 sage: ComplexHermitianEJA(2, prefix='z').gens()
1389 Our natural basis is normalized with respect to the natural inner
1390 product unless we specify otherwise::
1392 sage: set_random_seed()
1393 sage: J = ComplexHermitianEJA.random_instance()
1394 sage: all( b.norm() == 1 for b in J.gens() )
1397 Since our natural basis is normalized with respect to the natural
1398 inner product, and since we know that this algebra is an EJA, any
1399 left-multiplication operator's matrix will be symmetric because
1400 natural->EJA basis representation is an isometry and within the EJA
1401 the operator is self-adjoint by the Jordan axiom::
1403 sage: set_random_seed()
1404 sage: x = ComplexHermitianEJA.random_instance().random_element()
1405 sage: x.operator().matrix().is_symmetric()
1411 def _denormalized_basis(cls
, n
, field
):
1413 Returns a basis for the space of complex Hermitian n-by-n matrices.
1415 Why do we embed these? Basically, because all of numerical linear
1416 algebra assumes that you're working with vectors consisting of `n`
1417 entries from a field and scalars from the same field. There's no way
1418 to tell SageMath that (for example) the vectors contain complex
1419 numbers, while the scalar field is real.
1423 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1427 sage: set_random_seed()
1428 sage: n = ZZ.random_element(1,5)
1429 sage: field = QuadraticField(2, 'sqrt2')
1430 sage: B = ComplexHermitianEJA._denormalized_basis(n, field)
1431 sage: all( M.is_symmetric() for M in B)
1435 R
= PolynomialRing(field
, 'z')
1437 F
= field
.extension(z
**2 + 1, 'I')
1440 # This is like the symmetric case, but we need to be careful:
1442 # * We want conjugate-symmetry, not just symmetry.
1443 # * The diagonal will (as a result) be real.
1447 for j
in xrange(i
+1):
1448 Eij
= matrix(F
, n
, lambda k
,l
: k
==i
and l
==j
)
1450 Sij
= cls
.real_embed(Eij
)
1453 # The second one has a minus because it's conjugated.
1454 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
1456 Sij_imag
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
1459 # Since we embedded these, we can drop back to the "field" that we
1460 # started with instead of the complex extension "F".
1461 return ( s
.change_ring(field
) for s
in S
)
1464 def __init__(self
, n
, field
=QQ
, **kwargs
):
1465 basis
= self
._denormalized
_basis
(n
,field
)
1466 super(ComplexHermitianEJA
,self
).__init
__(field
, basis
, n
, **kwargs
)
1469 class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1473 Embed the n-by-n quaternion matrix ``M`` into the space of real
1474 matrices of size 4n-by-4n by first sending each quaternion entry `z
1475 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
1476 c+di],[-c + di, a-bi]]`, and then embedding those into a real
1481 sage: from mjo.eja.eja_algebra import \
1482 ....: QuaternionMatrixEuclideanJordanAlgebra
1486 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1487 sage: i,j,k = Q.gens()
1488 sage: x = 1 + 2*i + 3*j + 4*k
1489 sage: M = matrix(Q, 1, [[x]])
1490 sage: QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1496 Embedding is a homomorphism (isomorphism, in fact)::
1498 sage: set_random_seed()
1499 sage: n_max = QuaternionMatrixEuclideanJordanAlgebra._max_test_case_size()
1500 sage: n = ZZ.random_element(n_max)
1501 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1502 sage: X = random_matrix(Q, n)
1503 sage: Y = random_matrix(Q, n)
1504 sage: Xe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X)
1505 sage: Ye = QuaternionMatrixEuclideanJordanAlgebra.real_embed(Y)
1506 sage: XYe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1511 quaternions
= M
.base_ring()
1514 raise ValueError("the matrix 'M' must be square")
1516 F
= QuadraticField(-1, 'i')
1521 t
= z
.coefficient_tuple()
1526 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
1527 [-c
+ d
*i
, a
- b
*i
]])
1528 realM
= ComplexMatrixEuclideanJordanAlgebra
.real_embed(cplxM
)
1529 blocks
.append(realM
)
1531 # We should have real entries by now, so use the realest field
1532 # we've got for the return value.
1533 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
1538 def real_unembed(M
):
1540 The inverse of _embed_quaternion_matrix().
1544 sage: from mjo.eja.eja_algebra import \
1545 ....: QuaternionMatrixEuclideanJordanAlgebra
1549 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
1550 ....: [-2, 1, -4, 3],
1551 ....: [-3, 4, 1, -2],
1552 ....: [-4, -3, 2, 1]])
1553 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(M)
1554 [1 + 2*i + 3*j + 4*k]
1558 Unembedding is the inverse of embedding::
1560 sage: set_random_seed()
1561 sage: Q = QuaternionAlgebra(QQ, -1, -1)
1562 sage: M = random_matrix(Q, 3)
1563 sage: Me = QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1564 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1570 raise ValueError("the matrix 'M' must be square")
1571 if not n
.mod(4).is_zero():
1572 raise ValueError("the matrix 'M' must be a quaternion embedding")
1574 # Use the base ring of the matrix to ensure that its entries can be
1575 # multiplied by elements of the quaternion algebra.
1576 field
= M
.base_ring()
1577 Q
= QuaternionAlgebra(field
,-1,-1)
1580 # Go top-left to bottom-right (reading order), converting every
1581 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
1584 for l
in xrange(n
/4):
1585 for m
in xrange(n
/4):
1586 submat
= ComplexMatrixEuclideanJordanAlgebra
.real_unembed(
1587 M
[4*l
:4*l
+4,4*m
:4*m
+4] )
1588 if submat
[0,0] != submat
[1,1].conjugate():
1589 raise ValueError('bad on-diagonal submatrix')
1590 if submat
[0,1] != -submat
[1,0].conjugate():
1591 raise ValueError('bad off-diagonal submatrix')
1592 z
= submat
[0,0].vector()[0] # real part
1593 z
+= submat
[0,0].vector()[1]*i
# imag part
1594 z
+= submat
[0,1].vector()[0]*j
# real part
1595 z
+= submat
[0,1].vector()[1]*k
# imag part
1598 return matrix(Q
, n
/4, elements
)
1602 def natural_inner_product(cls
,X
,Y
):
1604 Compute a natural inner product in this algebra directly from
1609 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1613 This gives the same answer as the slow, default method implemented
1614 in :class:`MatrixEuclideanJordanAlgebra`::
1616 sage: set_random_seed()
1617 sage: J = QuaternionHermitianEJA.random_instance()
1618 sage: x,y = J.random_elements(2)
1619 sage: Xe = x.natural_representation()
1620 sage: Ye = y.natural_representation()
1621 sage: X = QuaternionHermitianEJA.real_unembed(Xe)
1622 sage: Y = QuaternionHermitianEJA.real_unembed(Ye)
1623 sage: expected = (X*Y).trace().coefficient_tuple()[0]
1624 sage: actual = QuaternionHermitianEJA.natural_inner_product(Xe,Ye)
1625 sage: actual == expected
1629 return RealMatrixEuclideanJordanAlgebra
.natural_inner_product(X
,Y
)/4
1632 class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra
,
1635 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
1636 matrices, the usual symmetric Jordan product, and the
1637 real-part-of-trace inner product. It has dimension `2n^2 - n` over
1642 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1646 In theory, our "field" can be any subfield of the reals::
1648 sage: QuaternionHermitianEJA(2, AA)
1649 Euclidean Jordan algebra of dimension 6 over Algebraic Real Field
1650 sage: QuaternionHermitianEJA(2, RR)
1651 Euclidean Jordan algebra of dimension 6 over Real Field with
1652 53 bits of precision
1656 The dimension of this algebra is `2*n^2 - n`::
1658 sage: set_random_seed()
1659 sage: n_max = QuaternionHermitianEJA._max_test_case_size()
1660 sage: n = ZZ.random_element(1, n_max)
1661 sage: J = QuaternionHermitianEJA(n)
1662 sage: J.dimension() == 2*(n^2) - n
1665 The Jordan multiplication is what we think it is::
1667 sage: set_random_seed()
1668 sage: J = QuaternionHermitianEJA.random_instance()
1669 sage: x,y = J.random_elements(2)
1670 sage: actual = (x*y).natural_representation()
1671 sage: X = x.natural_representation()
1672 sage: Y = y.natural_representation()
1673 sage: expected = (X*Y + Y*X)/2
1674 sage: actual == expected
1676 sage: J(expected) == x*y
1679 We can change the generator prefix::
1681 sage: QuaternionHermitianEJA(2, prefix='a').gens()
1682 (a0, a1, a2, a3, a4, a5)
1684 Our natural basis is normalized with respect to the natural inner
1685 product unless we specify otherwise::
1687 sage: set_random_seed()
1688 sage: J = QuaternionHermitianEJA.random_instance()
1689 sage: all( b.norm() == 1 for b in J.gens() )
1692 Since our natural basis is normalized with respect to the natural
1693 inner product, and since we know that this algebra is an EJA, any
1694 left-multiplication operator's matrix will be symmetric because
1695 natural->EJA basis representation is an isometry and within the EJA
1696 the operator is self-adjoint by the Jordan axiom::
1698 sage: set_random_seed()
1699 sage: x = QuaternionHermitianEJA.random_instance().random_element()
1700 sage: x.operator().matrix().is_symmetric()
1705 def _denormalized_basis(cls
, n
, field
):
1707 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
1709 Why do we embed these? Basically, because all of numerical
1710 linear algebra assumes that you're working with vectors consisting
1711 of `n` entries from a field and scalars from the same field. There's
1712 no way to tell SageMath that (for example) the vectors contain
1713 complex numbers, while the scalar field is real.
1717 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1721 sage: set_random_seed()
1722 sage: n = ZZ.random_element(1,5)
1723 sage: B = QuaternionHermitianEJA._denormalized_basis(n,QQ)
1724 sage: all( M.is_symmetric() for M in B )
1728 Q
= QuaternionAlgebra(QQ
,-1,-1)
1731 # This is like the symmetric case, but we need to be careful:
1733 # * We want conjugate-symmetry, not just symmetry.
1734 # * The diagonal will (as a result) be real.
1738 for j
in xrange(i
+1):
1739 Eij
= matrix(Q
, n
, lambda k
,l
: k
==i
and l
==j
)
1741 Sij
= cls
.real_embed(Eij
)
1744 # The second, third, and fourth ones have a minus
1745 # because they're conjugated.
1746 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
1748 Sij_I
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
1750 Sij_J
= cls
.real_embed(J
*Eij
- J
*Eij
.transpose())
1752 Sij_K
= cls
.real_embed(K
*Eij
- K
*Eij
.transpose())
1755 # Since we embedded these, we can drop back to the "field" that we
1756 # started with instead of the quaternion algebra "Q".
1757 return ( s
.change_ring(field
) for s
in S
)
1760 def __init__(self
, n
, field
=QQ
, **kwargs
):
1761 basis
= self
._denormalized
_basis
(n
,field
)
1762 super(QuaternionHermitianEJA
,self
).__init
__(field
, basis
, n
, **kwargs
)
1765 class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra
, KnownRankEJA
):
1767 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
1768 with the usual inner product and jordan product ``x*y =
1769 (<x_bar,y_bar>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
1774 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1778 This multiplication table can be verified by hand::
1780 sage: J = JordanSpinEJA(4)
1781 sage: e0,e1,e2,e3 = J.gens()
1797 We can change the generator prefix::
1799 sage: JordanSpinEJA(2, prefix='B').gens()
1803 def __init__(self
, n
, field
=QQ
, **kwargs
):
1804 V
= VectorSpace(field
, n
)
1805 mult_table
= [[V
.zero() for j
in xrange(n
)] for i
in xrange(n
)]
1815 z0
= x
.inner_product(y
)
1816 zbar
= y0
*xbar
+ x0
*ybar
1817 z
= V([z0
] + zbar
.list())
1818 mult_table
[i
][j
] = z
1820 # The rank of the spin algebra is two, unless we're in a
1821 # one-dimensional ambient space (because the rank is bounded by
1822 # the ambient dimension).
1823 fdeja
= super(JordanSpinEJA
, self
)
1824 return fdeja
.__init
__(field
, mult_table
, rank
=min(n
,2), **kwargs
)
1826 def inner_product(self
, x
, y
):
1828 Faster to reimplement than to use natural representations.
1832 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1836 Ensure that this is the usual inner product for the algebras
1839 sage: set_random_seed()
1840 sage: J = JordanSpinEJA.random_instance()
1841 sage: x,y = J.random_elements(2)
1842 sage: X = x.natural_representation()
1843 sage: Y = y.natural_representation()
1844 sage: x.inner_product(y) == J.natural_inner_product(X,Y)
1848 return x
.to_vector().inner_product(y
.to_vector())