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)
454 return self
.dimension() == 0
457 def multiplication_table(self
):
459 Return a visual representation of this algebra's multiplication
460 table (on basis elements).
464 sage: from mjo.eja.eja_algebra import JordanSpinEJA
468 sage: J = JordanSpinEJA(4)
469 sage: J.multiplication_table()
470 +----++----+----+----+----+
471 | * || e0 | e1 | e2 | e3 |
472 +====++====+====+====+====+
473 | e0 || e0 | e1 | e2 | e3 |
474 +----++----+----+----+----+
475 | e1 || e1 | e0 | 0 | 0 |
476 +----++----+----+----+----+
477 | e2 || e2 | 0 | e0 | 0 |
478 +----++----+----+----+----+
479 | e3 || e3 | 0 | 0 | e0 |
480 +----++----+----+----+----+
483 M
= list(self
._multiplication
_table
) # copy
484 for i
in xrange(len(M
)):
485 # M had better be "square"
486 M
[i
] = [self
.monomial(i
)] + M
[i
]
487 M
= [["*"] + list(self
.gens())] + M
488 return table(M
, header_row
=True, header_column
=True, frame
=True)
491 def natural_basis(self
):
493 Return a more-natural representation of this algebra's basis.
495 Every finite-dimensional Euclidean Jordan Algebra is a direct
496 sum of five simple algebras, four of which comprise Hermitian
497 matrices. This method returns the original "natural" basis
498 for our underlying vector space. (Typically, the natural basis
499 is used to construct the multiplication table in the first place.)
501 Note that this will always return a matrix. The standard basis
502 in `R^n` will be returned as `n`-by-`1` column matrices.
506 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
507 ....: RealSymmetricEJA)
511 sage: J = RealSymmetricEJA(2)
513 Finite family {0: e0, 1: e1, 2: e2}
514 sage: J.natural_basis()
516 [1 0] [ 0 1/2*sqrt2] [0 0]
517 [0 0], [1/2*sqrt2 0], [0 1]
522 sage: J = JordanSpinEJA(2)
524 Finite family {0: e0, 1: e1}
525 sage: J.natural_basis()
532 if self
._natural
_basis
is None:
533 M
= self
.natural_basis_space()
534 return tuple( M(b
.to_vector()) for b
in self
.basis() )
536 return self
._natural
_basis
539 def natural_basis_space(self
):
541 Return the matrix space in which this algebra's natural basis
544 if self
._natural
_basis
is None or len(self
._natural
_basis
) == 0:
545 return MatrixSpace(self
.base_ring(), self
.dimension(), 1)
547 return self
._natural
_basis
[0].matrix_space()
551 def natural_inner_product(X
,Y
):
553 Compute the inner product of two naturally-represented elements.
555 For example in the real symmetric matrix EJA, this will compute
556 the trace inner-product of two n-by-n symmetric matrices. The
557 default should work for the real cartesian product EJA, the
558 Jordan spin EJA, and the real symmetric matrices. The others
559 will have to be overridden.
561 return (X
.conjugate_transpose()*Y
).trace()
567 Return the unit element of this algebra.
571 sage: from mjo.eja.eja_algebra import (RealCartesianProductEJA,
576 sage: J = RealCartesianProductEJA(5)
578 e0 + e1 + e2 + e3 + e4
582 The identity element acts like the identity::
584 sage: set_random_seed()
585 sage: J = random_eja()
586 sage: x = J.random_element()
587 sage: J.one()*x == x and x*J.one() == x
590 The matrix of the unit element's operator is the identity::
592 sage: set_random_seed()
593 sage: J = random_eja()
594 sage: actual = J.one().operator().matrix()
595 sage: expected = matrix.identity(J.base_ring(), J.dimension())
596 sage: actual == expected
600 # We can brute-force compute the matrices of the operators
601 # that correspond to the basis elements of this algebra.
602 # If some linear combination of those basis elements is the
603 # algebra identity, then the same linear combination of
604 # their matrices has to be the identity matrix.
606 # Of course, matrices aren't vectors in sage, so we have to
607 # appeal to the "long vectors" isometry.
608 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
610 # Now we use basis linear algebra to find the coefficients,
611 # of the matrices-as-vectors-linear-combination, which should
612 # work for the original algebra basis too.
613 A
= matrix
.column(self
.base_ring(), oper_vecs
)
615 # We used the isometry on the left-hand side already, but we
616 # still need to do it for the right-hand side. Recall that we
617 # wanted something that summed to the identity matrix.
618 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
620 # Now if there's an identity element in the algebra, this should work.
621 coeffs
= A
.solve_right(b
)
622 return self
.linear_combination(zip(self
.gens(), coeffs
))
625 def random_elements(self
, count
):
627 Return ``count`` random elements as a tuple.
631 sage: from mjo.eja.eja_algebra import JordanSpinEJA
635 sage: J = JordanSpinEJA(3)
636 sage: x,y,z = J.random_elements(3)
637 sage: all( [ x in J, y in J, z in J ])
639 sage: len( J.random_elements(10) ) == 10
643 return tuple( self
.random_element() for idx
in xrange(count
) )
648 Return the rank of this EJA.
652 The author knows of no algorithm to compute the rank of an EJA
653 where only the multiplication table is known. In lieu of one, we
654 require the rank to be specified when the algebra is created,
655 and simply pass along that number here.
659 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
660 ....: RealSymmetricEJA,
661 ....: ComplexHermitianEJA,
662 ....: QuaternionHermitianEJA,
667 The rank of the Jordan spin algebra is always two::
669 sage: JordanSpinEJA(2).rank()
671 sage: JordanSpinEJA(3).rank()
673 sage: JordanSpinEJA(4).rank()
676 The rank of the `n`-by-`n` Hermitian real, complex, or
677 quaternion matrices is `n`::
679 sage: RealSymmetricEJA(4).rank()
681 sage: ComplexHermitianEJA(3).rank()
683 sage: QuaternionHermitianEJA(2).rank()
688 Ensure that every EJA that we know how to construct has a
689 positive integer rank::
691 sage: set_random_seed()
692 sage: r = random_eja().rank()
693 sage: r in ZZ and r > 0
700 def vector_space(self
):
702 Return the vector space that underlies this algebra.
706 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
710 sage: J = RealSymmetricEJA(2)
711 sage: J.vector_space()
712 Vector space of dimension 3 over...
715 return self
.zero().to_vector().parent().ambient_vector_space()
718 Element
= FiniteDimensionalEuclideanJordanAlgebraElement
721 class KnownRankEJA(object):
723 A class for algebras that we actually know we can construct. The
724 main issue is that, for most of our methods to make sense, we need
725 to know the rank of our algebra. Thus we can't simply generate a
726 "random" algebra, or even check that a given basis and product
727 satisfy the axioms; because even if everything looks OK, we wouldn't
728 know the rank we need to actuallty build the thing.
730 Not really a subclass of FDEJA because doing that causes method
731 resolution errors, e.g.
733 TypeError: Error when calling the metaclass bases
734 Cannot create a consistent method resolution
735 order (MRO) for bases FiniteDimensionalEuclideanJordanAlgebra,
740 def _max_test_case_size():
742 Return an integer "size" that is an upper bound on the size of
743 this algebra when it is used in a random test
744 case. Unfortunately, the term "size" is quite vague -- when
745 dealing with `R^n` under either the Hadamard or Jordan spin
746 product, the "size" refers to the dimension `n`. When dealing
747 with a matrix algebra (real symmetric or complex/quaternion
748 Hermitian), it refers to the size of the matrix, which is
749 far less than the dimension of the underlying vector space.
751 We default to five in this class, which is safe in `R^n`. The
752 matrix algebra subclasses (or any class where the "size" is
753 interpreted to be far less than the dimension) should override
754 with a smaller number.
759 def random_instance(cls
, field
=QQ
, **kwargs
):
761 Return a random instance of this type of algebra.
763 Beware, this will crash for "most instances" because the
764 constructor below looks wrong.
766 n
= ZZ
.random_element(cls
._max
_test
_case
_size
()) + 1
767 return cls(n
, field
, **kwargs
)
770 class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra
,
773 Return the Euclidean Jordan Algebra corresponding to the set
774 `R^n` under the Hadamard product.
776 Note: this is nothing more than the Cartesian product of ``n``
777 copies of the spin algebra. Once Cartesian product algebras
778 are implemented, this can go.
782 sage: from mjo.eja.eja_algebra import RealCartesianProductEJA
786 This multiplication table can be verified by hand::
788 sage: J = RealCartesianProductEJA(3)
789 sage: e0,e1,e2 = J.gens()
805 We can change the generator prefix::
807 sage: RealCartesianProductEJA(3, prefix='r').gens()
811 def __init__(self
, n
, field
=QQ
, **kwargs
):
812 V
= VectorSpace(field
, n
)
813 mult_table
= [ [ V
.gen(i
)*(i
== j
) for j
in xrange(n
) ]
816 fdeja
= super(RealCartesianProductEJA
, self
)
817 return fdeja
.__init
__(field
, mult_table
, rank
=n
, **kwargs
)
819 def inner_product(self
, x
, y
):
821 Faster to reimplement than to use natural representations.
825 sage: from mjo.eja.eja_algebra import RealCartesianProductEJA
829 Ensure that this is the usual inner product for the algebras
832 sage: set_random_seed()
833 sage: J = RealCartesianProductEJA.random_instance()
834 sage: x,y = J.random_elements(2)
835 sage: X = x.natural_representation()
836 sage: Y = y.natural_representation()
837 sage: x.inner_product(y) == J.natural_inner_product(X,Y)
841 return x
.to_vector().inner_product(y
.to_vector())
844 def random_eja(field
=QQ
):
846 Return a "random" finite-dimensional Euclidean Jordan Algebra.
850 For now, we choose a random natural number ``n`` (greater than zero)
851 and then give you back one of the following:
853 * The cartesian product of the rational numbers ``n`` times; this is
854 ``QQ^n`` with the Hadamard product.
856 * The Jordan spin algebra on ``QQ^n``.
858 * The ``n``-by-``n`` rational symmetric matrices with the symmetric
861 * The ``n``-by-``n`` complex-rational Hermitian matrices embedded
862 in the space of ``2n``-by-``2n`` real symmetric matrices.
864 * The ``n``-by-``n`` quaternion-rational Hermitian matrices embedded
865 in the space of ``4n``-by-``4n`` real symmetric matrices.
867 Later this might be extended to return Cartesian products of the
872 sage: from mjo.eja.eja_algebra import random_eja
877 Euclidean Jordan algebra of dimension...
880 classname
= choice(KnownRankEJA
.__subclasses
__())
881 return classname
.random_instance(field
=field
)
888 class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra
):
890 def _max_test_case_size():
891 # Play it safe, since this will be squared and the underlying
892 # field can have dimension 4 (quaternions) too.
895 def __init__(self
, field
, basis
, rank
, normalize_basis
=True, **kwargs
):
897 Compared to the superclass constructor, we take a basis instead of
898 a multiplication table because the latter can be computed in terms
899 of the former when the product is known (like it is here).
901 # Used in this class's fast _charpoly_coeff() override.
902 self
._basis
_normalizers
= None
904 # We're going to loop through this a few times, so now's a good
905 # time to ensure that it isn't a generator expression.
908 if rank
> 1 and normalize_basis
:
909 # We'll need sqrt(2) to normalize the basis, and this
910 # winds up in the multiplication table, so the whole
911 # algebra needs to be over the field extension.
912 R
= PolynomialRing(field
, 'z')
915 if p
.is_irreducible():
916 field
= field
.extension(p
, 'sqrt2', embedding
=RLF(2).sqrt())
917 basis
= tuple( s
.change_ring(field
) for s
in basis
)
918 self
._basis
_normalizers
= tuple(
919 ~
(self
.natural_inner_product(s
,s
).sqrt()) for s
in basis
)
920 basis
= tuple(s
*c
for (s
,c
) in izip(basis
,self
._basis
_normalizers
))
922 Qs
= self
.multiplication_table_from_matrix_basis(basis
)
924 fdeja
= super(MatrixEuclideanJordanAlgebra
, self
)
925 return fdeja
.__init
__(field
,
933 def _charpoly_coeff(self
, i
):
935 Override the parent method with something that tries to compute
936 over a faster (non-extension) field.
938 if self
._basis
_normalizers
is None:
939 # We didn't normalize, so assume that the basis we started
940 # with had entries in a nice field.
941 return super(MatrixEuclideanJordanAlgebra
, self
)._charpoly
_coeff
(i
)
943 basis
= ( (b
/n
) for (b
,n
) in izip(self
.natural_basis(),
944 self
._basis
_normalizers
) )
946 # Do this over the rationals and convert back at the end.
947 J
= MatrixEuclideanJordanAlgebra(QQ
,
950 normalize_basis
=False)
951 (_
,x
,_
,_
) = J
._charpoly
_matrix
_system
()
952 p
= J
._charpoly
_coeff
(i
)
953 # p might be missing some vars, have to substitute "optionally"
954 pairs
= izip(x
.base_ring().gens(), self
._basis
_normalizers
)
955 substitutions
= { v: v*c for (v,c) in pairs }
956 result
= p
.subs(substitutions
)
958 # The result of "subs" can be either a coefficient-ring
959 # element or a polynomial. Gotta handle both cases.
961 return self
.base_ring()(result
)
963 return result
.change_ring(self
.base_ring())
967 def multiplication_table_from_matrix_basis(basis
):
969 At least three of the five simple Euclidean Jordan algebras have the
970 symmetric multiplication (A,B) |-> (AB + BA)/2, where the
971 multiplication on the right is matrix multiplication. Given a basis
972 for the underlying matrix space, this function returns a
973 multiplication table (obtained by looping through the basis
974 elements) for an algebra of those matrices.
976 # In S^2, for example, we nominally have four coordinates even
977 # though the space is of dimension three only. The vector space V
978 # is supposed to hold the entire long vector, and the subspace W
979 # of V will be spanned by the vectors that arise from symmetric
980 # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
981 field
= basis
[0].base_ring()
982 dimension
= basis
[0].nrows()
984 V
= VectorSpace(field
, dimension
**2)
985 W
= V
.span_of_basis( _mat2vec(s
) for s
in basis
)
987 mult_table
= [[W
.zero() for j
in xrange(n
)] for i
in xrange(n
)]
990 mat_entry
= (basis
[i
]*basis
[j
] + basis
[j
]*basis
[i
])/2
991 mult_table
[i
][j
] = W
.coordinate_vector(_mat2vec(mat_entry
))
999 Embed the matrix ``M`` into a space of real matrices.
1001 The matrix ``M`` can have entries in any field at the moment:
1002 the real numbers, complex numbers, or quaternions. And although
1003 they are not a field, we can probably support octonions at some
1004 point, too. This function returns a real matrix that "acts like"
1005 the original with respect to matrix multiplication; i.e.
1007 real_embed(M*N) = real_embed(M)*real_embed(N)
1010 raise NotImplementedError
1014 def real_unembed(M
):
1016 The inverse of :meth:`real_embed`.
1018 raise NotImplementedError
1022 def natural_inner_product(cls
,X
,Y
):
1023 Xu
= cls
.real_unembed(X
)
1024 Yu
= cls
.real_unembed(Y
)
1025 tr
= (Xu
*Yu
).trace()
1028 # It's real already.
1031 # Otherwise, try the thing that works for complex numbers; and
1032 # if that doesn't work, the thing that works for quaternions.
1034 return tr
.vector()[0] # real part, imag part is index 1
1035 except AttributeError:
1036 # A quaternions doesn't have a vector() method, but does
1037 # have coefficient_tuple() method that returns the
1038 # coefficients of 1, i, j, and k -- in that order.
1039 return tr
.coefficient_tuple()[0]
1042 class RealMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1046 The identity function, for embedding real matrices into real
1052 def real_unembed(M
):
1054 The identity function, for unembedding real matrices from real
1060 class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra
, KnownRankEJA
):
1062 The rank-n simple EJA consisting of real symmetric n-by-n
1063 matrices, the usual symmetric Jordan product, and the trace inner
1064 product. It has dimension `(n^2 + n)/2` over the reals.
1068 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1072 sage: J = RealSymmetricEJA(2)
1073 sage: e0, e1, e2 = J.gens()
1081 In theory, our "field" can be any subfield of the reals::
1083 sage: RealSymmetricEJA(2, AA)
1084 Euclidean Jordan algebra of dimension 3 over Algebraic Real Field
1085 sage: RealSymmetricEJA(2, RR)
1086 Euclidean Jordan algebra of dimension 3 over Real Field with
1087 53 bits of precision
1091 The dimension of this algebra is `(n^2 + n) / 2`::
1093 sage: set_random_seed()
1094 sage: n_max = RealSymmetricEJA._max_test_case_size()
1095 sage: n = ZZ.random_element(1, n_max)
1096 sage: J = RealSymmetricEJA(n)
1097 sage: J.dimension() == (n^2 + n)/2
1100 The Jordan multiplication is what we think it is::
1102 sage: set_random_seed()
1103 sage: J = RealSymmetricEJA.random_instance()
1104 sage: x,y = J.random_elements(2)
1105 sage: actual = (x*y).natural_representation()
1106 sage: X = x.natural_representation()
1107 sage: Y = y.natural_representation()
1108 sage: expected = (X*Y + Y*X)/2
1109 sage: actual == expected
1111 sage: J(expected) == x*y
1114 We can change the generator prefix::
1116 sage: RealSymmetricEJA(3, prefix='q').gens()
1117 (q0, q1, q2, q3, q4, q5)
1119 Our natural basis is normalized with respect to the natural inner
1120 product unless we specify otherwise::
1122 sage: set_random_seed()
1123 sage: J = RealSymmetricEJA.random_instance()
1124 sage: all( b.norm() == 1 for b in J.gens() )
1127 Since our natural basis is normalized with respect to the natural
1128 inner product, and since we know that this algebra is an EJA, any
1129 left-multiplication operator's matrix will be symmetric because
1130 natural->EJA basis representation is an isometry and within the EJA
1131 the operator is self-adjoint by the Jordan axiom::
1133 sage: set_random_seed()
1134 sage: x = RealSymmetricEJA.random_instance().random_element()
1135 sage: x.operator().matrix().is_symmetric()
1140 def _denormalized_basis(cls
, n
, field
):
1142 Return a basis for the space of real symmetric n-by-n matrices.
1146 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1150 sage: set_random_seed()
1151 sage: n = ZZ.random_element(1,5)
1152 sage: B = RealSymmetricEJA._denormalized_basis(n,QQ)
1153 sage: all( M.is_symmetric() for M in B)
1157 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1161 for j
in xrange(i
+1):
1162 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
1166 Sij
= Eij
+ Eij
.transpose()
1172 def _max_test_case_size():
1173 return 4 # Dimension 10
1176 def __init__(self
, n
, field
=QQ
, **kwargs
):
1177 basis
= self
._denormalized
_basis
(n
, field
)
1178 super(RealSymmetricEJA
, self
).__init
__(field
, basis
, n
, **kwargs
)
1181 class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1185 Embed the n-by-n complex matrix ``M`` into the space of real
1186 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1187 bi` to the block matrix ``[[a,b],[-b,a]]``.
1191 sage: from mjo.eja.eja_algebra import \
1192 ....: ComplexMatrixEuclideanJordanAlgebra
1196 sage: F = QuadraticField(-1, 'i')
1197 sage: x1 = F(4 - 2*i)
1198 sage: x2 = F(1 + 2*i)
1201 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1202 sage: ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1211 Embedding is a homomorphism (isomorphism, in fact)::
1213 sage: set_random_seed()
1214 sage: n_max = ComplexMatrixEuclideanJordanAlgebra._max_test_case_size()
1215 sage: n = ZZ.random_element(n_max)
1216 sage: F = QuadraticField(-1, 'i')
1217 sage: X = random_matrix(F, n)
1218 sage: Y = random_matrix(F, n)
1219 sage: Xe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X)
1220 sage: Ye = ComplexMatrixEuclideanJordanAlgebra.real_embed(Y)
1221 sage: XYe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1228 raise ValueError("the matrix 'M' must be square")
1230 # We don't need any adjoined elements...
1231 field
= M
.base_ring().base_ring()
1235 a
= z
.list()[0] # real part, I guess
1236 b
= z
.list()[1] # imag part, I guess
1237 blocks
.append(matrix(field
, 2, [[a
,b
],[-b
,a
]]))
1239 return matrix
.block(field
, n
, blocks
)
1243 def real_unembed(M
):
1245 The inverse of _embed_complex_matrix().
1249 sage: from mjo.eja.eja_algebra import \
1250 ....: ComplexMatrixEuclideanJordanAlgebra
1254 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1255 ....: [-2, 1, -4, 3],
1256 ....: [ 9, 10, 11, 12],
1257 ....: [-10, 9, -12, 11] ])
1258 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(A)
1260 [ 10*i + 9 12*i + 11]
1264 Unembedding is the inverse of embedding::
1266 sage: set_random_seed()
1267 sage: F = QuadraticField(-1, 'i')
1268 sage: M = random_matrix(F, 3)
1269 sage: Me = ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1270 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1276 raise ValueError("the matrix 'M' must be square")
1277 if not n
.mod(2).is_zero():
1278 raise ValueError("the matrix 'M' must be a complex embedding")
1280 # If "M" was normalized, its base ring might have roots
1281 # adjoined and they can stick around after unembedding.
1282 field
= M
.base_ring()
1283 R
= PolynomialRing(field
, 'z')
1285 F
= field
.extension(z
**2 + 1, 'i', embedding
=CLF(-1).sqrt())
1288 # Go top-left to bottom-right (reading order), converting every
1289 # 2-by-2 block we see to a single complex element.
1291 for k
in xrange(n
/2):
1292 for j
in xrange(n
/2):
1293 submat
= M
[2*k
:2*k
+2,2*j
:2*j
+2]
1294 if submat
[0,0] != submat
[1,1]:
1295 raise ValueError('bad on-diagonal submatrix')
1296 if submat
[0,1] != -submat
[1,0]:
1297 raise ValueError('bad off-diagonal submatrix')
1298 z
= submat
[0,0] + submat
[0,1]*i
1301 return matrix(F
, n
/2, elements
)
1305 def natural_inner_product(cls
,X
,Y
):
1307 Compute a natural inner product in this algebra directly from
1312 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1316 This gives the same answer as the slow, default method implemented
1317 in :class:`MatrixEuclideanJordanAlgebra`::
1319 sage: set_random_seed()
1320 sage: J = ComplexHermitianEJA.random_instance()
1321 sage: x,y = J.random_elements(2)
1322 sage: Xe = x.natural_representation()
1323 sage: Ye = y.natural_representation()
1324 sage: X = ComplexHermitianEJA.real_unembed(Xe)
1325 sage: Y = ComplexHermitianEJA.real_unembed(Ye)
1326 sage: expected = (X*Y).trace().vector()[0]
1327 sage: actual = ComplexHermitianEJA.natural_inner_product(Xe,Ye)
1328 sage: actual == expected
1332 return RealMatrixEuclideanJordanAlgebra
.natural_inner_product(X
,Y
)/2
1335 class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra
, KnownRankEJA
):
1337 The rank-n simple EJA consisting of complex Hermitian n-by-n
1338 matrices over the real numbers, the usual symmetric Jordan product,
1339 and the real-part-of-trace inner product. It has dimension `n^2` over
1344 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1348 In theory, our "field" can be any subfield of the reals::
1350 sage: ComplexHermitianEJA(2, AA)
1351 Euclidean Jordan algebra of dimension 4 over Algebraic Real Field
1352 sage: ComplexHermitianEJA(2, RR)
1353 Euclidean Jordan algebra of dimension 4 over Real Field with
1354 53 bits of precision
1358 The dimension of this algebra is `n^2`::
1360 sage: set_random_seed()
1361 sage: n_max = ComplexHermitianEJA._max_test_case_size()
1362 sage: n = ZZ.random_element(1, n_max)
1363 sage: J = ComplexHermitianEJA(n)
1364 sage: J.dimension() == n^2
1367 The Jordan multiplication is what we think it is::
1369 sage: set_random_seed()
1370 sage: J = ComplexHermitianEJA.random_instance()
1371 sage: x,y = J.random_elements(2)
1372 sage: actual = (x*y).natural_representation()
1373 sage: X = x.natural_representation()
1374 sage: Y = y.natural_representation()
1375 sage: expected = (X*Y + Y*X)/2
1376 sage: actual == expected
1378 sage: J(expected) == x*y
1381 We can change the generator prefix::
1383 sage: ComplexHermitianEJA(2, prefix='z').gens()
1386 Our natural basis is normalized with respect to the natural inner
1387 product unless we specify otherwise::
1389 sage: set_random_seed()
1390 sage: J = ComplexHermitianEJA.random_instance()
1391 sage: all( b.norm() == 1 for b in J.gens() )
1394 Since our natural basis is normalized with respect to the natural
1395 inner product, and since we know that this algebra is an EJA, any
1396 left-multiplication operator's matrix will be symmetric because
1397 natural->EJA basis representation is an isometry and within the EJA
1398 the operator is self-adjoint by the Jordan axiom::
1400 sage: set_random_seed()
1401 sage: x = ComplexHermitianEJA.random_instance().random_element()
1402 sage: x.operator().matrix().is_symmetric()
1408 def _denormalized_basis(cls
, n
, field
):
1410 Returns a basis for the space of complex Hermitian n-by-n matrices.
1412 Why do we embed these? Basically, because all of numerical linear
1413 algebra assumes that you're working with vectors consisting of `n`
1414 entries from a field and scalars from the same field. There's no way
1415 to tell SageMath that (for example) the vectors contain complex
1416 numbers, while the scalar field is real.
1420 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1424 sage: set_random_seed()
1425 sage: n = ZZ.random_element(1,5)
1426 sage: field = QuadraticField(2, 'sqrt2')
1427 sage: B = ComplexHermitianEJA._denormalized_basis(n, field)
1428 sage: all( M.is_symmetric() for M in B)
1432 R
= PolynomialRing(field
, 'z')
1434 F
= field
.extension(z
**2 + 1, 'I')
1437 # This is like the symmetric case, but we need to be careful:
1439 # * We want conjugate-symmetry, not just symmetry.
1440 # * The diagonal will (as a result) be real.
1444 for j
in xrange(i
+1):
1445 Eij
= matrix(F
, n
, lambda k
,l
: k
==i
and l
==j
)
1447 Sij
= cls
.real_embed(Eij
)
1450 # The second one has a minus because it's conjugated.
1451 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
1453 Sij_imag
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
1456 # Since we embedded these, we can drop back to the "field" that we
1457 # started with instead of the complex extension "F".
1458 return ( s
.change_ring(field
) for s
in S
)
1461 def __init__(self
, n
, field
=QQ
, **kwargs
):
1462 basis
= self
._denormalized
_basis
(n
,field
)
1463 super(ComplexHermitianEJA
,self
).__init
__(field
, basis
, n
, **kwargs
)
1466 class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1470 Embed the n-by-n quaternion matrix ``M`` into the space of real
1471 matrices of size 4n-by-4n by first sending each quaternion entry `z
1472 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
1473 c+di],[-c + di, a-bi]]`, and then embedding those into a real
1478 sage: from mjo.eja.eja_algebra import \
1479 ....: QuaternionMatrixEuclideanJordanAlgebra
1483 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1484 sage: i,j,k = Q.gens()
1485 sage: x = 1 + 2*i + 3*j + 4*k
1486 sage: M = matrix(Q, 1, [[x]])
1487 sage: QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1493 Embedding is a homomorphism (isomorphism, in fact)::
1495 sage: set_random_seed()
1496 sage: n_max = QuaternionMatrixEuclideanJordanAlgebra._max_test_case_size()
1497 sage: n = ZZ.random_element(n_max)
1498 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1499 sage: X = random_matrix(Q, n)
1500 sage: Y = random_matrix(Q, n)
1501 sage: Xe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X)
1502 sage: Ye = QuaternionMatrixEuclideanJordanAlgebra.real_embed(Y)
1503 sage: XYe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1508 quaternions
= M
.base_ring()
1511 raise ValueError("the matrix 'M' must be square")
1513 F
= QuadraticField(-1, 'i')
1518 t
= z
.coefficient_tuple()
1523 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
1524 [-c
+ d
*i
, a
- b
*i
]])
1525 realM
= ComplexMatrixEuclideanJordanAlgebra
.real_embed(cplxM
)
1526 blocks
.append(realM
)
1528 # We should have real entries by now, so use the realest field
1529 # we've got for the return value.
1530 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
1535 def real_unembed(M
):
1537 The inverse of _embed_quaternion_matrix().
1541 sage: from mjo.eja.eja_algebra import \
1542 ....: QuaternionMatrixEuclideanJordanAlgebra
1546 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
1547 ....: [-2, 1, -4, 3],
1548 ....: [-3, 4, 1, -2],
1549 ....: [-4, -3, 2, 1]])
1550 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(M)
1551 [1 + 2*i + 3*j + 4*k]
1555 Unembedding is the inverse of embedding::
1557 sage: set_random_seed()
1558 sage: Q = QuaternionAlgebra(QQ, -1, -1)
1559 sage: M = random_matrix(Q, 3)
1560 sage: Me = QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1561 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1567 raise ValueError("the matrix 'M' must be square")
1568 if not n
.mod(4).is_zero():
1569 raise ValueError("the matrix 'M' must be a quaternion embedding")
1571 # Use the base ring of the matrix to ensure that its entries can be
1572 # multiplied by elements of the quaternion algebra.
1573 field
= M
.base_ring()
1574 Q
= QuaternionAlgebra(field
,-1,-1)
1577 # Go top-left to bottom-right (reading order), converting every
1578 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
1581 for l
in xrange(n
/4):
1582 for m
in xrange(n
/4):
1583 submat
= ComplexMatrixEuclideanJordanAlgebra
.real_unembed(
1584 M
[4*l
:4*l
+4,4*m
:4*m
+4] )
1585 if submat
[0,0] != submat
[1,1].conjugate():
1586 raise ValueError('bad on-diagonal submatrix')
1587 if submat
[0,1] != -submat
[1,0].conjugate():
1588 raise ValueError('bad off-diagonal submatrix')
1589 z
= submat
[0,0].vector()[0] # real part
1590 z
+= submat
[0,0].vector()[1]*i
# imag part
1591 z
+= submat
[0,1].vector()[0]*j
# real part
1592 z
+= submat
[0,1].vector()[1]*k
# imag part
1595 return matrix(Q
, n
/4, elements
)
1599 def natural_inner_product(cls
,X
,Y
):
1601 Compute a natural inner product in this algebra directly from
1606 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1610 This gives the same answer as the slow, default method implemented
1611 in :class:`MatrixEuclideanJordanAlgebra`::
1613 sage: set_random_seed()
1614 sage: J = QuaternionHermitianEJA.random_instance()
1615 sage: x,y = J.random_elements(2)
1616 sage: Xe = x.natural_representation()
1617 sage: Ye = y.natural_representation()
1618 sage: X = QuaternionHermitianEJA.real_unembed(Xe)
1619 sage: Y = QuaternionHermitianEJA.real_unembed(Ye)
1620 sage: expected = (X*Y).trace().coefficient_tuple()[0]
1621 sage: actual = QuaternionHermitianEJA.natural_inner_product(Xe,Ye)
1622 sage: actual == expected
1626 return RealMatrixEuclideanJordanAlgebra
.natural_inner_product(X
,Y
)/4
1629 class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra
,
1632 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
1633 matrices, the usual symmetric Jordan product, and the
1634 real-part-of-trace inner product. It has dimension `2n^2 - n` over
1639 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1643 In theory, our "field" can be any subfield of the reals::
1645 sage: QuaternionHermitianEJA(2, AA)
1646 Euclidean Jordan algebra of dimension 6 over Algebraic Real Field
1647 sage: QuaternionHermitianEJA(2, RR)
1648 Euclidean Jordan algebra of dimension 6 over Real Field with
1649 53 bits of precision
1653 The dimension of this algebra is `2*n^2 - n`::
1655 sage: set_random_seed()
1656 sage: n_max = QuaternionHermitianEJA._max_test_case_size()
1657 sage: n = ZZ.random_element(1, n_max)
1658 sage: J = QuaternionHermitianEJA(n)
1659 sage: J.dimension() == 2*(n^2) - n
1662 The Jordan multiplication is what we think it is::
1664 sage: set_random_seed()
1665 sage: J = QuaternionHermitianEJA.random_instance()
1666 sage: x,y = J.random_elements(2)
1667 sage: actual = (x*y).natural_representation()
1668 sage: X = x.natural_representation()
1669 sage: Y = y.natural_representation()
1670 sage: expected = (X*Y + Y*X)/2
1671 sage: actual == expected
1673 sage: J(expected) == x*y
1676 We can change the generator prefix::
1678 sage: QuaternionHermitianEJA(2, prefix='a').gens()
1679 (a0, a1, a2, a3, a4, a5)
1681 Our natural basis is normalized with respect to the natural inner
1682 product unless we specify otherwise::
1684 sage: set_random_seed()
1685 sage: J = QuaternionHermitianEJA.random_instance()
1686 sage: all( b.norm() == 1 for b in J.gens() )
1689 Since our natural basis is normalized with respect to the natural
1690 inner product, and since we know that this algebra is an EJA, any
1691 left-multiplication operator's matrix will be symmetric because
1692 natural->EJA basis representation is an isometry and within the EJA
1693 the operator is self-adjoint by the Jordan axiom::
1695 sage: set_random_seed()
1696 sage: x = QuaternionHermitianEJA.random_instance().random_element()
1697 sage: x.operator().matrix().is_symmetric()
1702 def _denormalized_basis(cls
, n
, field
):
1704 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
1706 Why do we embed these? Basically, because all of numerical
1707 linear algebra assumes that you're working with vectors consisting
1708 of `n` entries from a field and scalars from the same field. There's
1709 no way to tell SageMath that (for example) the vectors contain
1710 complex numbers, while the scalar field is real.
1714 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1718 sage: set_random_seed()
1719 sage: n = ZZ.random_element(1,5)
1720 sage: B = QuaternionHermitianEJA._denormalized_basis(n,QQ)
1721 sage: all( M.is_symmetric() for M in B )
1725 Q
= QuaternionAlgebra(QQ
,-1,-1)
1728 # This is like the symmetric case, but we need to be careful:
1730 # * We want conjugate-symmetry, not just symmetry.
1731 # * The diagonal will (as a result) be real.
1735 for j
in xrange(i
+1):
1736 Eij
= matrix(Q
, n
, lambda k
,l
: k
==i
and l
==j
)
1738 Sij
= cls
.real_embed(Eij
)
1741 # The second, third, and fourth ones have a minus
1742 # because they're conjugated.
1743 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
1745 Sij_I
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
1747 Sij_J
= cls
.real_embed(J
*Eij
- J
*Eij
.transpose())
1749 Sij_K
= cls
.real_embed(K
*Eij
- K
*Eij
.transpose())
1752 # Since we embedded these, we can drop back to the "field" that we
1753 # started with instead of the quaternion algebra "Q".
1754 return ( s
.change_ring(field
) for s
in S
)
1757 def __init__(self
, n
, field
=QQ
, **kwargs
):
1758 basis
= self
._denormalized
_basis
(n
,field
)
1759 super(QuaternionHermitianEJA
,self
).__init
__(field
, basis
, n
, **kwargs
)
1762 class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra
, KnownRankEJA
):
1764 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
1765 with the usual inner product and jordan product ``x*y =
1766 (<x_bar,y_bar>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
1771 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1775 This multiplication table can be verified by hand::
1777 sage: J = JordanSpinEJA(4)
1778 sage: e0,e1,e2,e3 = J.gens()
1794 We can change the generator prefix::
1796 sage: JordanSpinEJA(2, prefix='B').gens()
1800 def __init__(self
, n
, field
=QQ
, **kwargs
):
1801 V
= VectorSpace(field
, n
)
1802 mult_table
= [[V
.zero() for j
in xrange(n
)] for i
in xrange(n
)]
1812 z0
= x
.inner_product(y
)
1813 zbar
= y0
*xbar
+ x0
*ybar
1814 z
= V([z0
] + zbar
.list())
1815 mult_table
[i
][j
] = z
1817 # The rank of the spin algebra is two, unless we're in a
1818 # one-dimensional ambient space (because the rank is bounded by
1819 # the ambient dimension).
1820 fdeja
= super(JordanSpinEJA
, self
)
1821 return fdeja
.__init
__(field
, mult_table
, rank
=min(n
,2), **kwargs
)
1823 def inner_product(self
, x
, y
):
1825 Faster to reimplement than to use natural representations.
1829 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1833 Ensure that this is the usual inner product for the algebras
1836 sage: set_random_seed()
1837 sage: J = JordanSpinEJA.random_instance()
1838 sage: x,y = J.random_elements(2)
1839 sage: X = x.natural_representation()
1840 sage: Y = y.natural_representation()
1841 sage: x.inner_product(y) == J.natural_inner_product(X,Y)
1845 return x
.to_vector().inner_product(y
.to_vector())