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, unless the algebra is trivial in
690 which case its rank will be zero::
692 sage: set_random_seed()
693 sage: J = random_eja()
697 sage: r > 0 or (r == 0 and J.is_trivial())
704 def vector_space(self
):
706 Return the vector space that underlies this algebra.
710 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
714 sage: J = RealSymmetricEJA(2)
715 sage: J.vector_space()
716 Vector space of dimension 3 over...
719 return self
.zero().to_vector().parent().ambient_vector_space()
722 Element
= FiniteDimensionalEuclideanJordanAlgebraElement
725 class KnownRankEJA(object):
727 A class for algebras that we actually know we can construct. The
728 main issue is that, for most of our methods to make sense, we need
729 to know the rank of our algebra. Thus we can't simply generate a
730 "random" algebra, or even check that a given basis and product
731 satisfy the axioms; because even if everything looks OK, we wouldn't
732 know the rank we need to actuallty build the thing.
734 Not really a subclass of FDEJA because doing that causes method
735 resolution errors, e.g.
737 TypeError: Error when calling the metaclass bases
738 Cannot create a consistent method resolution
739 order (MRO) for bases FiniteDimensionalEuclideanJordanAlgebra,
744 def _max_test_case_size():
746 Return an integer "size" that is an upper bound on the size of
747 this algebra when it is used in a random test
748 case. Unfortunately, the term "size" is quite vague -- when
749 dealing with `R^n` under either the Hadamard or Jordan spin
750 product, the "size" refers to the dimension `n`. When dealing
751 with a matrix algebra (real symmetric or complex/quaternion
752 Hermitian), it refers to the size of the matrix, which is
753 far less than the dimension of the underlying vector space.
755 We default to five in this class, which is safe in `R^n`. The
756 matrix algebra subclasses (or any class where the "size" is
757 interpreted to be far less than the dimension) should override
758 with a smaller number.
763 def random_instance(cls
, field
=QQ
, **kwargs
):
765 Return a random instance of this type of algebra.
767 Beware, this will crash for "most instances" because the
768 constructor below looks wrong.
770 n
= ZZ
.random_element(cls
._max
_test
_case
_size
()) + 1
771 return cls(n
, field
, **kwargs
)
774 class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra
,
777 Return the Euclidean Jordan Algebra corresponding to the set
778 `R^n` under the Hadamard product.
780 Note: this is nothing more than the Cartesian product of ``n``
781 copies of the spin algebra. Once Cartesian product algebras
782 are implemented, this can go.
786 sage: from mjo.eja.eja_algebra import RealCartesianProductEJA
790 This multiplication table can be verified by hand::
792 sage: J = RealCartesianProductEJA(3)
793 sage: e0,e1,e2 = J.gens()
809 We can change the generator prefix::
811 sage: RealCartesianProductEJA(3, prefix='r').gens()
815 def __init__(self
, n
, field
=QQ
, **kwargs
):
816 V
= VectorSpace(field
, n
)
817 mult_table
= [ [ V
.gen(i
)*(i
== j
) for j
in xrange(n
) ]
820 fdeja
= super(RealCartesianProductEJA
, self
)
821 return fdeja
.__init
__(field
, mult_table
, rank
=n
, **kwargs
)
823 def inner_product(self
, x
, y
):
825 Faster to reimplement than to use natural representations.
829 sage: from mjo.eja.eja_algebra import RealCartesianProductEJA
833 Ensure that this is the usual inner product for the algebras
836 sage: set_random_seed()
837 sage: J = RealCartesianProductEJA.random_instance()
838 sage: x,y = J.random_elements(2)
839 sage: X = x.natural_representation()
840 sage: Y = y.natural_representation()
841 sage: x.inner_product(y) == J.natural_inner_product(X,Y)
845 return x
.to_vector().inner_product(y
.to_vector())
848 def random_eja(field
=QQ
):
850 Return a "random" finite-dimensional Euclidean Jordan Algebra.
854 For now, we choose a random natural number ``n`` (greater than zero)
855 and then give you back one of the following:
857 * The cartesian product of the rational numbers ``n`` times; this is
858 ``QQ^n`` with the Hadamard product.
860 * The Jordan spin algebra on ``QQ^n``.
862 * The ``n``-by-``n`` rational symmetric matrices with the symmetric
865 * The ``n``-by-``n`` complex-rational Hermitian matrices embedded
866 in the space of ``2n``-by-``2n`` real symmetric matrices.
868 * The ``n``-by-``n`` quaternion-rational Hermitian matrices embedded
869 in the space of ``4n``-by-``4n`` real symmetric matrices.
871 Later this might be extended to return Cartesian products of the
876 sage: from mjo.eja.eja_algebra import random_eja
881 Euclidean Jordan algebra of dimension...
884 classname
= choice(KnownRankEJA
.__subclasses
__())
885 return classname
.random_instance(field
=field
)
892 class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra
):
894 def _max_test_case_size():
895 # Play it safe, since this will be squared and the underlying
896 # field can have dimension 4 (quaternions) too.
899 def __init__(self
, field
, basis
, rank
, normalize_basis
=True, **kwargs
):
901 Compared to the superclass constructor, we take a basis instead of
902 a multiplication table because the latter can be computed in terms
903 of the former when the product is known (like it is here).
905 # Used in this class's fast _charpoly_coeff() override.
906 self
._basis
_normalizers
= None
908 # We're going to loop through this a few times, so now's a good
909 # time to ensure that it isn't a generator expression.
912 if rank
> 1 and normalize_basis
:
913 # We'll need sqrt(2) to normalize the basis, and this
914 # winds up in the multiplication table, so the whole
915 # algebra needs to be over the field extension.
916 R
= PolynomialRing(field
, 'z')
919 if p
.is_irreducible():
920 field
= field
.extension(p
, 'sqrt2', embedding
=RLF(2).sqrt())
921 basis
= tuple( s
.change_ring(field
) for s
in basis
)
922 self
._basis
_normalizers
= tuple(
923 ~
(self
.natural_inner_product(s
,s
).sqrt()) for s
in basis
)
924 basis
= tuple(s
*c
for (s
,c
) in izip(basis
,self
._basis
_normalizers
))
926 Qs
= self
.multiplication_table_from_matrix_basis(basis
)
928 fdeja
= super(MatrixEuclideanJordanAlgebra
, self
)
929 return fdeja
.__init
__(field
,
937 def _charpoly_coeff(self
, i
):
939 Override the parent method with something that tries to compute
940 over a faster (non-extension) field.
942 if self
._basis
_normalizers
is None:
943 # We didn't normalize, so assume that the basis we started
944 # with had entries in a nice field.
945 return super(MatrixEuclideanJordanAlgebra
, self
)._charpoly
_coeff
(i
)
947 basis
= ( (b
/n
) for (b
,n
) in izip(self
.natural_basis(),
948 self
._basis
_normalizers
) )
950 # Do this over the rationals and convert back at the end.
951 J
= MatrixEuclideanJordanAlgebra(QQ
,
954 normalize_basis
=False)
955 (_
,x
,_
,_
) = J
._charpoly
_matrix
_system
()
956 p
= J
._charpoly
_coeff
(i
)
957 # p might be missing some vars, have to substitute "optionally"
958 pairs
= izip(x
.base_ring().gens(), self
._basis
_normalizers
)
959 substitutions
= { v: v*c for (v,c) in pairs }
960 result
= p
.subs(substitutions
)
962 # The result of "subs" can be either a coefficient-ring
963 # element or a polynomial. Gotta handle both cases.
965 return self
.base_ring()(result
)
967 return result
.change_ring(self
.base_ring())
971 def multiplication_table_from_matrix_basis(basis
):
973 At least three of the five simple Euclidean Jordan algebras have the
974 symmetric multiplication (A,B) |-> (AB + BA)/2, where the
975 multiplication on the right is matrix multiplication. Given a basis
976 for the underlying matrix space, this function returns a
977 multiplication table (obtained by looping through the basis
978 elements) for an algebra of those matrices.
980 # In S^2, for example, we nominally have four coordinates even
981 # though the space is of dimension three only. The vector space V
982 # is supposed to hold the entire long vector, and the subspace W
983 # of V will be spanned by the vectors that arise from symmetric
984 # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
985 field
= basis
[0].base_ring()
986 dimension
= basis
[0].nrows()
988 V
= VectorSpace(field
, dimension
**2)
989 W
= V
.span_of_basis( _mat2vec(s
) for s
in basis
)
991 mult_table
= [[W
.zero() for j
in xrange(n
)] for i
in xrange(n
)]
994 mat_entry
= (basis
[i
]*basis
[j
] + basis
[j
]*basis
[i
])/2
995 mult_table
[i
][j
] = W
.coordinate_vector(_mat2vec(mat_entry
))
1003 Embed the matrix ``M`` into a space of real matrices.
1005 The matrix ``M`` can have entries in any field at the moment:
1006 the real numbers, complex numbers, or quaternions. And although
1007 they are not a field, we can probably support octonions at some
1008 point, too. This function returns a real matrix that "acts like"
1009 the original with respect to matrix multiplication; i.e.
1011 real_embed(M*N) = real_embed(M)*real_embed(N)
1014 raise NotImplementedError
1018 def real_unembed(M
):
1020 The inverse of :meth:`real_embed`.
1022 raise NotImplementedError
1026 def natural_inner_product(cls
,X
,Y
):
1027 Xu
= cls
.real_unembed(X
)
1028 Yu
= cls
.real_unembed(Y
)
1029 tr
= (Xu
*Yu
).trace()
1032 # It's real already.
1035 # Otherwise, try the thing that works for complex numbers; and
1036 # if that doesn't work, the thing that works for quaternions.
1038 return tr
.vector()[0] # real part, imag part is index 1
1039 except AttributeError:
1040 # A quaternions doesn't have a vector() method, but does
1041 # have coefficient_tuple() method that returns the
1042 # coefficients of 1, i, j, and k -- in that order.
1043 return tr
.coefficient_tuple()[0]
1046 class RealMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1050 The identity function, for embedding real matrices into real
1056 def real_unembed(M
):
1058 The identity function, for unembedding real matrices from real
1064 class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra
, KnownRankEJA
):
1066 The rank-n simple EJA consisting of real symmetric n-by-n
1067 matrices, the usual symmetric Jordan product, and the trace inner
1068 product. It has dimension `(n^2 + n)/2` over the reals.
1072 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1076 sage: J = RealSymmetricEJA(2)
1077 sage: e0, e1, e2 = J.gens()
1085 In theory, our "field" can be any subfield of the reals::
1087 sage: RealSymmetricEJA(2, AA)
1088 Euclidean Jordan algebra of dimension 3 over Algebraic Real Field
1089 sage: RealSymmetricEJA(2, RR)
1090 Euclidean Jordan algebra of dimension 3 over Real Field with
1091 53 bits of precision
1095 The dimension of this algebra is `(n^2 + n) / 2`::
1097 sage: set_random_seed()
1098 sage: n_max = RealSymmetricEJA._max_test_case_size()
1099 sage: n = ZZ.random_element(1, n_max)
1100 sage: J = RealSymmetricEJA(n)
1101 sage: J.dimension() == (n^2 + n)/2
1104 The Jordan multiplication is what we think it is::
1106 sage: set_random_seed()
1107 sage: J = RealSymmetricEJA.random_instance()
1108 sage: x,y = J.random_elements(2)
1109 sage: actual = (x*y).natural_representation()
1110 sage: X = x.natural_representation()
1111 sage: Y = y.natural_representation()
1112 sage: expected = (X*Y + Y*X)/2
1113 sage: actual == expected
1115 sage: J(expected) == x*y
1118 We can change the generator prefix::
1120 sage: RealSymmetricEJA(3, prefix='q').gens()
1121 (q0, q1, q2, q3, q4, q5)
1123 Our natural basis is normalized with respect to the natural inner
1124 product unless we specify otherwise::
1126 sage: set_random_seed()
1127 sage: J = RealSymmetricEJA.random_instance()
1128 sage: all( b.norm() == 1 for b in J.gens() )
1131 Since our natural basis is normalized with respect to the natural
1132 inner product, and since we know that this algebra is an EJA, any
1133 left-multiplication operator's matrix will be symmetric because
1134 natural->EJA basis representation is an isometry and within the EJA
1135 the operator is self-adjoint by the Jordan axiom::
1137 sage: set_random_seed()
1138 sage: x = RealSymmetricEJA.random_instance().random_element()
1139 sage: x.operator().matrix().is_symmetric()
1144 def _denormalized_basis(cls
, n
, field
):
1146 Return a basis for the space of real symmetric n-by-n matrices.
1150 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1154 sage: set_random_seed()
1155 sage: n = ZZ.random_element(1,5)
1156 sage: B = RealSymmetricEJA._denormalized_basis(n,QQ)
1157 sage: all( M.is_symmetric() for M in B)
1161 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1165 for j
in xrange(i
+1):
1166 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
1170 Sij
= Eij
+ Eij
.transpose()
1176 def _max_test_case_size():
1177 return 4 # Dimension 10
1180 def __init__(self
, n
, field
=QQ
, **kwargs
):
1181 basis
= self
._denormalized
_basis
(n
, field
)
1182 super(RealSymmetricEJA
, self
).__init
__(field
, basis
, n
, **kwargs
)
1185 class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1189 Embed the n-by-n complex matrix ``M`` into the space of real
1190 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1191 bi` to the block matrix ``[[a,b],[-b,a]]``.
1195 sage: from mjo.eja.eja_algebra import \
1196 ....: ComplexMatrixEuclideanJordanAlgebra
1200 sage: F = QuadraticField(-1, 'i')
1201 sage: x1 = F(4 - 2*i)
1202 sage: x2 = F(1 + 2*i)
1205 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1206 sage: ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1215 Embedding is a homomorphism (isomorphism, in fact)::
1217 sage: set_random_seed()
1218 sage: n_max = ComplexMatrixEuclideanJordanAlgebra._max_test_case_size()
1219 sage: n = ZZ.random_element(n_max)
1220 sage: F = QuadraticField(-1, 'i')
1221 sage: X = random_matrix(F, n)
1222 sage: Y = random_matrix(F, n)
1223 sage: Xe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X)
1224 sage: Ye = ComplexMatrixEuclideanJordanAlgebra.real_embed(Y)
1225 sage: XYe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1232 raise ValueError("the matrix 'M' must be square")
1234 # We don't need any adjoined elements...
1235 field
= M
.base_ring().base_ring()
1239 a
= z
.list()[0] # real part, I guess
1240 b
= z
.list()[1] # imag part, I guess
1241 blocks
.append(matrix(field
, 2, [[a
,b
],[-b
,a
]]))
1243 return matrix
.block(field
, n
, blocks
)
1247 def real_unembed(M
):
1249 The inverse of _embed_complex_matrix().
1253 sage: from mjo.eja.eja_algebra import \
1254 ....: ComplexMatrixEuclideanJordanAlgebra
1258 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1259 ....: [-2, 1, -4, 3],
1260 ....: [ 9, 10, 11, 12],
1261 ....: [-10, 9, -12, 11] ])
1262 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(A)
1264 [ 10*i + 9 12*i + 11]
1268 Unembedding is the inverse of embedding::
1270 sage: set_random_seed()
1271 sage: F = QuadraticField(-1, 'i')
1272 sage: M = random_matrix(F, 3)
1273 sage: Me = ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1274 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1280 raise ValueError("the matrix 'M' must be square")
1281 if not n
.mod(2).is_zero():
1282 raise ValueError("the matrix 'M' must be a complex embedding")
1284 # If "M" was normalized, its base ring might have roots
1285 # adjoined and they can stick around after unembedding.
1286 field
= M
.base_ring()
1287 R
= PolynomialRing(field
, 'z')
1289 F
= field
.extension(z
**2 + 1, 'i', embedding
=CLF(-1).sqrt())
1292 # Go top-left to bottom-right (reading order), converting every
1293 # 2-by-2 block we see to a single complex element.
1295 for k
in xrange(n
/2):
1296 for j
in xrange(n
/2):
1297 submat
= M
[2*k
:2*k
+2,2*j
:2*j
+2]
1298 if submat
[0,0] != submat
[1,1]:
1299 raise ValueError('bad on-diagonal submatrix')
1300 if submat
[0,1] != -submat
[1,0]:
1301 raise ValueError('bad off-diagonal submatrix')
1302 z
= submat
[0,0] + submat
[0,1]*i
1305 return matrix(F
, n
/2, elements
)
1309 def natural_inner_product(cls
,X
,Y
):
1311 Compute a natural inner product in this algebra directly from
1316 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1320 This gives the same answer as the slow, default method implemented
1321 in :class:`MatrixEuclideanJordanAlgebra`::
1323 sage: set_random_seed()
1324 sage: J = ComplexHermitianEJA.random_instance()
1325 sage: x,y = J.random_elements(2)
1326 sage: Xe = x.natural_representation()
1327 sage: Ye = y.natural_representation()
1328 sage: X = ComplexHermitianEJA.real_unembed(Xe)
1329 sage: Y = ComplexHermitianEJA.real_unembed(Ye)
1330 sage: expected = (X*Y).trace().vector()[0]
1331 sage: actual = ComplexHermitianEJA.natural_inner_product(Xe,Ye)
1332 sage: actual == expected
1336 return RealMatrixEuclideanJordanAlgebra
.natural_inner_product(X
,Y
)/2
1339 class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra
, KnownRankEJA
):
1341 The rank-n simple EJA consisting of complex Hermitian n-by-n
1342 matrices over the real numbers, the usual symmetric Jordan product,
1343 and the real-part-of-trace inner product. It has dimension `n^2` over
1348 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1352 In theory, our "field" can be any subfield of the reals::
1354 sage: ComplexHermitianEJA(2, AA)
1355 Euclidean Jordan algebra of dimension 4 over Algebraic Real Field
1356 sage: ComplexHermitianEJA(2, RR)
1357 Euclidean Jordan algebra of dimension 4 over Real Field with
1358 53 bits of precision
1362 The dimension of this algebra is `n^2`::
1364 sage: set_random_seed()
1365 sage: n_max = ComplexHermitianEJA._max_test_case_size()
1366 sage: n = ZZ.random_element(1, n_max)
1367 sage: J = ComplexHermitianEJA(n)
1368 sage: J.dimension() == n^2
1371 The Jordan multiplication is what we think it is::
1373 sage: set_random_seed()
1374 sage: J = ComplexHermitianEJA.random_instance()
1375 sage: x,y = J.random_elements(2)
1376 sage: actual = (x*y).natural_representation()
1377 sage: X = x.natural_representation()
1378 sage: Y = y.natural_representation()
1379 sage: expected = (X*Y + Y*X)/2
1380 sage: actual == expected
1382 sage: J(expected) == x*y
1385 We can change the generator prefix::
1387 sage: ComplexHermitianEJA(2, prefix='z').gens()
1390 Our natural basis is normalized with respect to the natural inner
1391 product unless we specify otherwise::
1393 sage: set_random_seed()
1394 sage: J = ComplexHermitianEJA.random_instance()
1395 sage: all( b.norm() == 1 for b in J.gens() )
1398 Since our natural basis is normalized with respect to the natural
1399 inner product, and since we know that this algebra is an EJA, any
1400 left-multiplication operator's matrix will be symmetric because
1401 natural->EJA basis representation is an isometry and within the EJA
1402 the operator is self-adjoint by the Jordan axiom::
1404 sage: set_random_seed()
1405 sage: x = ComplexHermitianEJA.random_instance().random_element()
1406 sage: x.operator().matrix().is_symmetric()
1412 def _denormalized_basis(cls
, n
, field
):
1414 Returns a basis for the space of complex Hermitian n-by-n matrices.
1416 Why do we embed these? Basically, because all of numerical linear
1417 algebra assumes that you're working with vectors consisting of `n`
1418 entries from a field and scalars from the same field. There's no way
1419 to tell SageMath that (for example) the vectors contain complex
1420 numbers, while the scalar field is real.
1424 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1428 sage: set_random_seed()
1429 sage: n = ZZ.random_element(1,5)
1430 sage: field = QuadraticField(2, 'sqrt2')
1431 sage: B = ComplexHermitianEJA._denormalized_basis(n, field)
1432 sage: all( M.is_symmetric() for M in B)
1436 R
= PolynomialRing(field
, 'z')
1438 F
= field
.extension(z
**2 + 1, 'I')
1441 # This is like the symmetric case, but we need to be careful:
1443 # * We want conjugate-symmetry, not just symmetry.
1444 # * The diagonal will (as a result) be real.
1448 for j
in xrange(i
+1):
1449 Eij
= matrix(F
, n
, lambda k
,l
: k
==i
and l
==j
)
1451 Sij
= cls
.real_embed(Eij
)
1454 # The second one has a minus because it's conjugated.
1455 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
1457 Sij_imag
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
1460 # Since we embedded these, we can drop back to the "field" that we
1461 # started with instead of the complex extension "F".
1462 return ( s
.change_ring(field
) for s
in S
)
1465 def __init__(self
, n
, field
=QQ
, **kwargs
):
1466 basis
= self
._denormalized
_basis
(n
,field
)
1467 super(ComplexHermitianEJA
,self
).__init
__(field
, basis
, n
, **kwargs
)
1470 class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1474 Embed the n-by-n quaternion matrix ``M`` into the space of real
1475 matrices of size 4n-by-4n by first sending each quaternion entry `z
1476 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
1477 c+di],[-c + di, a-bi]]`, and then embedding those into a real
1482 sage: from mjo.eja.eja_algebra import \
1483 ....: QuaternionMatrixEuclideanJordanAlgebra
1487 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1488 sage: i,j,k = Q.gens()
1489 sage: x = 1 + 2*i + 3*j + 4*k
1490 sage: M = matrix(Q, 1, [[x]])
1491 sage: QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1497 Embedding is a homomorphism (isomorphism, in fact)::
1499 sage: set_random_seed()
1500 sage: n_max = QuaternionMatrixEuclideanJordanAlgebra._max_test_case_size()
1501 sage: n = ZZ.random_element(n_max)
1502 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1503 sage: X = random_matrix(Q, n)
1504 sage: Y = random_matrix(Q, n)
1505 sage: Xe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X)
1506 sage: Ye = QuaternionMatrixEuclideanJordanAlgebra.real_embed(Y)
1507 sage: XYe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1512 quaternions
= M
.base_ring()
1515 raise ValueError("the matrix 'M' must be square")
1517 F
= QuadraticField(-1, 'i')
1522 t
= z
.coefficient_tuple()
1527 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
1528 [-c
+ d
*i
, a
- b
*i
]])
1529 realM
= ComplexMatrixEuclideanJordanAlgebra
.real_embed(cplxM
)
1530 blocks
.append(realM
)
1532 # We should have real entries by now, so use the realest field
1533 # we've got for the return value.
1534 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
1539 def real_unembed(M
):
1541 The inverse of _embed_quaternion_matrix().
1545 sage: from mjo.eja.eja_algebra import \
1546 ....: QuaternionMatrixEuclideanJordanAlgebra
1550 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
1551 ....: [-2, 1, -4, 3],
1552 ....: [-3, 4, 1, -2],
1553 ....: [-4, -3, 2, 1]])
1554 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(M)
1555 [1 + 2*i + 3*j + 4*k]
1559 Unembedding is the inverse of embedding::
1561 sage: set_random_seed()
1562 sage: Q = QuaternionAlgebra(QQ, -1, -1)
1563 sage: M = random_matrix(Q, 3)
1564 sage: Me = QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1565 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1571 raise ValueError("the matrix 'M' must be square")
1572 if not n
.mod(4).is_zero():
1573 raise ValueError("the matrix 'M' must be a quaternion embedding")
1575 # Use the base ring of the matrix to ensure that its entries can be
1576 # multiplied by elements of the quaternion algebra.
1577 field
= M
.base_ring()
1578 Q
= QuaternionAlgebra(field
,-1,-1)
1581 # Go top-left to bottom-right (reading order), converting every
1582 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
1585 for l
in xrange(n
/4):
1586 for m
in xrange(n
/4):
1587 submat
= ComplexMatrixEuclideanJordanAlgebra
.real_unembed(
1588 M
[4*l
:4*l
+4,4*m
:4*m
+4] )
1589 if submat
[0,0] != submat
[1,1].conjugate():
1590 raise ValueError('bad on-diagonal submatrix')
1591 if submat
[0,1] != -submat
[1,0].conjugate():
1592 raise ValueError('bad off-diagonal submatrix')
1593 z
= submat
[0,0].vector()[0] # real part
1594 z
+= submat
[0,0].vector()[1]*i
# imag part
1595 z
+= submat
[0,1].vector()[0]*j
# real part
1596 z
+= submat
[0,1].vector()[1]*k
# imag part
1599 return matrix(Q
, n
/4, elements
)
1603 def natural_inner_product(cls
,X
,Y
):
1605 Compute a natural inner product in this algebra directly from
1610 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1614 This gives the same answer as the slow, default method implemented
1615 in :class:`MatrixEuclideanJordanAlgebra`::
1617 sage: set_random_seed()
1618 sage: J = QuaternionHermitianEJA.random_instance()
1619 sage: x,y = J.random_elements(2)
1620 sage: Xe = x.natural_representation()
1621 sage: Ye = y.natural_representation()
1622 sage: X = QuaternionHermitianEJA.real_unembed(Xe)
1623 sage: Y = QuaternionHermitianEJA.real_unembed(Ye)
1624 sage: expected = (X*Y).trace().coefficient_tuple()[0]
1625 sage: actual = QuaternionHermitianEJA.natural_inner_product(Xe,Ye)
1626 sage: actual == expected
1630 return RealMatrixEuclideanJordanAlgebra
.natural_inner_product(X
,Y
)/4
1633 class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra
,
1636 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
1637 matrices, the usual symmetric Jordan product, and the
1638 real-part-of-trace inner product. It has dimension `2n^2 - n` over
1643 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1647 In theory, our "field" can be any subfield of the reals::
1649 sage: QuaternionHermitianEJA(2, AA)
1650 Euclidean Jordan algebra of dimension 6 over Algebraic Real Field
1651 sage: QuaternionHermitianEJA(2, RR)
1652 Euclidean Jordan algebra of dimension 6 over Real Field with
1653 53 bits of precision
1657 The dimension of this algebra is `2*n^2 - n`::
1659 sage: set_random_seed()
1660 sage: n_max = QuaternionHermitianEJA._max_test_case_size()
1661 sage: n = ZZ.random_element(1, n_max)
1662 sage: J = QuaternionHermitianEJA(n)
1663 sage: J.dimension() == 2*(n^2) - n
1666 The Jordan multiplication is what we think it is::
1668 sage: set_random_seed()
1669 sage: J = QuaternionHermitianEJA.random_instance()
1670 sage: x,y = J.random_elements(2)
1671 sage: actual = (x*y).natural_representation()
1672 sage: X = x.natural_representation()
1673 sage: Y = y.natural_representation()
1674 sage: expected = (X*Y + Y*X)/2
1675 sage: actual == expected
1677 sage: J(expected) == x*y
1680 We can change the generator prefix::
1682 sage: QuaternionHermitianEJA(2, prefix='a').gens()
1683 (a0, a1, a2, a3, a4, a5)
1685 Our natural basis is normalized with respect to the natural inner
1686 product unless we specify otherwise::
1688 sage: set_random_seed()
1689 sage: J = QuaternionHermitianEJA.random_instance()
1690 sage: all( b.norm() == 1 for b in J.gens() )
1693 Since our natural basis is normalized with respect to the natural
1694 inner product, and since we know that this algebra is an EJA, any
1695 left-multiplication operator's matrix will be symmetric because
1696 natural->EJA basis representation is an isometry and within the EJA
1697 the operator is self-adjoint by the Jordan axiom::
1699 sage: set_random_seed()
1700 sage: x = QuaternionHermitianEJA.random_instance().random_element()
1701 sage: x.operator().matrix().is_symmetric()
1706 def _denormalized_basis(cls
, n
, field
):
1708 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
1710 Why do we embed these? Basically, because all of numerical
1711 linear algebra assumes that you're working with vectors consisting
1712 of `n` entries from a field and scalars from the same field. There's
1713 no way to tell SageMath that (for example) the vectors contain
1714 complex numbers, while the scalar field is real.
1718 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1722 sage: set_random_seed()
1723 sage: n = ZZ.random_element(1,5)
1724 sage: B = QuaternionHermitianEJA._denormalized_basis(n,QQ)
1725 sage: all( M.is_symmetric() for M in B )
1729 Q
= QuaternionAlgebra(QQ
,-1,-1)
1732 # This is like the symmetric case, but we need to be careful:
1734 # * We want conjugate-symmetry, not just symmetry.
1735 # * The diagonal will (as a result) be real.
1739 for j
in xrange(i
+1):
1740 Eij
= matrix(Q
, n
, lambda k
,l
: k
==i
and l
==j
)
1742 Sij
= cls
.real_embed(Eij
)
1745 # The second, third, and fourth ones have a minus
1746 # because they're conjugated.
1747 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
1749 Sij_I
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
1751 Sij_J
= cls
.real_embed(J
*Eij
- J
*Eij
.transpose())
1753 Sij_K
= cls
.real_embed(K
*Eij
- K
*Eij
.transpose())
1756 # Since we embedded these, we can drop back to the "field" that we
1757 # started with instead of the quaternion algebra "Q".
1758 return ( s
.change_ring(field
) for s
in S
)
1761 def __init__(self
, n
, field
=QQ
, **kwargs
):
1762 basis
= self
._denormalized
_basis
(n
,field
)
1763 super(QuaternionHermitianEJA
,self
).__init
__(field
, basis
, n
, **kwargs
)
1766 class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra
, KnownRankEJA
):
1768 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
1769 with the usual inner product and jordan product ``x*y =
1770 (<x_bar,y_bar>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
1775 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1779 This multiplication table can be verified by hand::
1781 sage: J = JordanSpinEJA(4)
1782 sage: e0,e1,e2,e3 = J.gens()
1798 We can change the generator prefix::
1800 sage: JordanSpinEJA(2, prefix='B').gens()
1804 def __init__(self
, n
, field
=QQ
, **kwargs
):
1805 V
= VectorSpace(field
, n
)
1806 mult_table
= [[V
.zero() for j
in xrange(n
)] for i
in xrange(n
)]
1816 z0
= x
.inner_product(y
)
1817 zbar
= y0
*xbar
+ x0
*ybar
1818 z
= V([z0
] + zbar
.list())
1819 mult_table
[i
][j
] = z
1821 # The rank of the spin algebra is two, unless we're in a
1822 # one-dimensional ambient space (because the rank is bounded by
1823 # the ambient dimension).
1824 fdeja
= super(JordanSpinEJA
, self
)
1825 return fdeja
.__init
__(field
, mult_table
, rank
=min(n
,2), **kwargs
)
1827 def inner_product(self
, x
, y
):
1829 Faster to reimplement than to use natural representations.
1833 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1837 Ensure that this is the usual inner product for the algebras
1840 sage: set_random_seed()
1841 sage: J = JordanSpinEJA.random_instance()
1842 sage: x,y = J.random_elements(2)
1843 sage: X = x.natural_representation()
1844 sage: Y = y.natural_representation()
1845 sage: x.inner_product(y) == J.natural_inner_product(X,Y)
1849 return x
.to_vector().inner_product(y
.to_vector())