2 Euclidean Jordan Algebras. These are formally-real Jordan Algebras;
3 specifically those where u^2 + v^2 = 0 implies that u = v = 0. They
4 are used in optimization, and have some additional nice methods beyond
5 what can be supported in a general Jordan Algebra.
8 from itertools
import repeat
10 from sage
.algebras
.quatalg
.quaternion_algebra
import QuaternionAlgebra
11 from sage
.categories
.magmatic_algebras
import MagmaticAlgebras
12 from sage
.combinat
.free_module
import CombinatorialFreeModule
13 from sage
.matrix
.constructor
import matrix
14 from sage
.matrix
.matrix_space
import MatrixSpace
15 from sage
.misc
.cachefunc
import cached_method
16 from sage
.misc
.prandom
import choice
17 from sage
.misc
.table
import table
18 from sage
.modules
.free_module
import FreeModule
, VectorSpace
19 from sage
.rings
.integer_ring
import ZZ
20 from sage
.rings
.number_field
.number_field
import NumberField
, QuadraticField
21 from sage
.rings
.polynomial
.polynomial_ring_constructor
import PolynomialRing
22 from sage
.rings
.rational_field
import QQ
23 from sage
.rings
.real_lazy
import CLF
, RLF
24 from sage
.structure
.element
import is_Matrix
26 from mjo
.eja
.eja_element
import FiniteDimensionalEuclideanJordanAlgebraElement
27 from mjo
.eja
.eja_utils
import _mat2vec
29 class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule
):
30 # This is an ugly hack needed to prevent the category framework
31 # from implementing a coercion from our base ring (e.g. the
32 # rationals) into the algebra. First of all -- such a coercion is
33 # nonsense to begin with. But more importantly, it tries to do so
34 # in the category of rings, and since our algebras aren't
35 # associative they generally won't be rings.
36 _no_generic_basering_coercion
= True
48 sage: from mjo.eja.eja_algebra import random_eja
52 By definition, Jordan multiplication commutes::
54 sage: set_random_seed()
55 sage: J = random_eja()
56 sage: x,y = J.random_elements(2)
62 self
._natural
_basis
= natural_basis
64 # TODO: HACK for the charpoly.. needs redesign badly.
65 self
._basis
_normalizers
= None
68 category
= MagmaticAlgebras(field
).FiniteDimensional()
69 category
= category
.WithBasis().Unital()
71 fda
= super(FiniteDimensionalEuclideanJordanAlgebra
, self
)
73 range(len(mult_table
)),
76 self
.print_options(bracket
='')
78 # The multiplication table we're given is necessarily in terms
79 # of vectors, because we don't have an algebra yet for
80 # anything to be an element of. However, it's faster in the
81 # long run to have the multiplication table be in terms of
82 # algebra elements. We do this after calling the superclass
83 # constructor so that from_vector() knows what to do.
84 self
._multiplication
_table
= [ map(lambda x
: self
.from_vector(x
), ls
)
85 for ls
in mult_table
]
88 def _element_constructor_(self
, elt
):
90 Construct an element of this algebra from its natural
93 This gets called only after the parent element _call_ method
94 fails to find a coercion for the argument.
98 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
99 ....: RealCartesianProductEJA,
100 ....: RealSymmetricEJA)
104 The identity in `S^n` is converted to the identity in the EJA::
106 sage: J = RealSymmetricEJA(3)
107 sage: I = matrix.identity(QQ,3)
108 sage: J(I) == J.one()
111 This skew-symmetric matrix can't be represented in the EJA::
113 sage: J = RealSymmetricEJA(3)
114 sage: A = matrix(QQ,3, lambda i,j: i-j)
116 Traceback (most recent call last):
118 ArithmeticError: vector is not in free module
122 Ensure that we can convert any element of the two non-matrix
123 simple algebras (whose natural representations are their usual
124 vector representations) back and forth faithfully::
126 sage: set_random_seed()
127 sage: J = RealCartesianProductEJA.random_instance()
128 sage: x = J.random_element()
129 sage: J(x.to_vector().column()) == x
131 sage: J = JordanSpinEJA.random_instance()
132 sage: x = J.random_element()
133 sage: J(x.to_vector().column()) == x
138 # The superclass implementation of random_element()
139 # needs to be able to coerce "0" into the algebra.
142 natural_basis
= self
.natural_basis()
143 basis_space
= natural_basis
[0].matrix_space()
144 if elt
not in basis_space
:
145 raise ValueError("not a naturally-represented algebra element")
147 # Thanks for nothing! Matrix spaces aren't vector spaces in
148 # Sage, so we have to figure out its natural-basis coordinates
149 # ourselves. We use the basis space's ring instead of the
150 # element's ring because the basis space might be an algebraic
151 # closure whereas the base ring of the 3-by-3 identity matrix
152 # could be QQ instead of QQbar.
153 V
= VectorSpace(basis_space
.base_ring(), elt
.nrows()*elt
.ncols())
154 W
= V
.span_of_basis( _mat2vec(s
) for s
in natural_basis
)
155 coords
= W
.coordinate_vector(_mat2vec(elt
))
156 return self
.from_vector(coords
)
160 def _max_test_case_size():
162 Return an integer "size" that is an upper bound on the size of
163 this algebra when it is used in a random test
164 case. Unfortunately, the term "size" is quite vague -- when
165 dealing with `R^n` under either the Hadamard or Jordan spin
166 product, the "size" refers to the dimension `n`. When dealing
167 with a matrix algebra (real symmetric or complex/quaternion
168 Hermitian), it refers to the size of the matrix, which is
169 far less than the dimension of the underlying vector space.
171 We default to five in this class, which is safe in `R^n`. The
172 matrix algebra subclasses (or any class where the "size" is
173 interpreted to be far less than the dimension) should override
174 with a smaller number.
181 Return a string representation of ``self``.
185 sage: from mjo.eja.eja_algebra import JordanSpinEJA
189 Ensure that it says what we think it says::
191 sage: JordanSpinEJA(2, field=QQ)
192 Euclidean Jordan algebra of dimension 2 over Rational Field
193 sage: JordanSpinEJA(3, field=RDF)
194 Euclidean Jordan algebra of dimension 3 over Real Double Field
197 fmt
= "Euclidean Jordan algebra of dimension {} over {}"
198 return fmt
.format(self
.dimension(), self
.base_ring())
200 def product_on_basis(self
, i
, j
):
201 return self
._multiplication
_table
[i
][j
]
203 def _a_regular_element(self
):
205 Guess a regular element. Needed to compute the basis for our
206 characteristic polynomial coefficients.
210 sage: from mjo.eja.eja_algebra import random_eja
214 Ensure that this hacky method succeeds for every algebra that we
215 know how to construct::
217 sage: set_random_seed()
218 sage: J = random_eja()
219 sage: J._a_regular_element().is_regular()
224 z
= self
.sum( (i
+1)*gs
[i
] for i
in range(len(gs
)) )
225 if not z
.is_regular():
226 raise ValueError("don't know a regular element")
231 def _charpoly_basis_space(self
):
233 Return the vector space spanned by the basis used in our
234 characteristic polynomial coefficients. This is used not only to
235 compute those coefficients, but also any time we need to
236 evaluate the coefficients (like when we compute the trace or
239 z
= self
._a
_regular
_element
()
240 # Don't use the parent vector space directly here in case this
241 # happens to be a subalgebra. In that case, we would be e.g.
242 # two-dimensional but span_of_basis() would expect three
244 V
= VectorSpace(self
.base_ring(), self
.vector_space().dimension())
245 basis
= [ (z
**k
).to_vector() for k
in range(self
.rank()) ]
246 V1
= V
.span_of_basis( basis
)
247 b
= (V1
.basis() + V1
.complement().basis())
248 return V
.span_of_basis(b
)
253 def _charpoly_coeff(self
, i
):
255 Return the coefficient polynomial "a_{i}" of this algebra's
256 general characteristic polynomial.
258 Having this be a separate cached method lets us compute and
259 store the trace/determinant (a_{r-1} and a_{0} respectively)
260 separate from the entire characteristic polynomial.
262 if self
._basis
_normalizers
is not None:
263 # Must be a matrix class?
264 # WARNING/TODO: this whole mess is mis-designed.
265 n
= self
.natural_basis_space().nrows()
266 field
= self
.base_ring().base_ring() # yeeeeaaaahhh
267 J
= self
.__class
__(n
, field
, False)
268 (_
,x
,_
,_
) = J
._charpoly
_matrix
_system
()
269 p
= J
._charpoly
_coeff
(i
)
270 # p might be missing some vars, have to substitute "optionally"
271 pairs
= zip(x
.base_ring().gens(), self
._basis
_normalizers
)
272 substitutions
= { v: v*c for (v,c) in pairs }
273 return p
.subs(substitutions
)
275 (A_of_x
, x
, xr
, detA
) = self
._charpoly
_matrix
_system
()
276 R
= A_of_x
.base_ring()
278 # Guaranteed by theory
281 # Danger: the in-place modification is done for performance
282 # reasons (reconstructing a matrix with huge polynomial
283 # entries is slow), but I don't know how cached_method works,
284 # so it's highly possible that we're modifying some global
285 # list variable by reference, here. In other words, you
286 # probably shouldn't call this method twice on the same
287 # algebra, at the same time, in two threads
288 Ai_orig
= A_of_x
.column(i
)
289 A_of_x
.set_column(i
,xr
)
290 numerator
= A_of_x
.det()
291 A_of_x
.set_column(i
,Ai_orig
)
293 # We're relying on the theory here to ensure that each a_i is
294 # indeed back in R, and the added negative signs are to make
295 # the whole charpoly expression sum to zero.
296 return R(-numerator
/detA
)
300 def _charpoly_matrix_system(self
):
302 Compute the matrix whose entries A_ij are polynomials in
303 X1,...,XN, the vector ``x`` of variables X1,...,XN, the vector
304 corresponding to `x^r` and the determinent of the matrix A =
305 [A_ij]. In other words, all of the fixed (cachable) data needed
306 to compute the coefficients of the characteristic polynomial.
311 # Turn my vector space into a module so that "vectors" can
312 # have multivatiate polynomial entries.
313 names
= tuple('X' + str(i
) for i
in range(1,n
+1))
314 R
= PolynomialRing(self
.base_ring(), names
)
316 # Using change_ring() on the parent's vector space doesn't work
317 # here because, in a subalgebra, that vector space has a basis
318 # and change_ring() tries to bring the basis along with it. And
319 # that doesn't work unless the new ring is a PID, which it usually
323 # Now let x = (X1,X2,...,Xn) be the vector whose entries are
327 # And figure out the "left multiplication by x" matrix in
330 monomial_matrices
= [ self
.monomial(i
).operator().matrix()
331 for i
in range(n
) ] # don't recompute these!
333 ek
= self
.monomial(k
).to_vector()
335 sum( x
[i
]*(monomial_matrices
[i
]*ek
)
336 for i
in range(n
) ) )
337 Lx
= matrix
.column(R
, lmbx_cols
)
339 # Now we can compute powers of x "symbolically"
340 x_powers
= [self
.one().to_vector(), x
]
341 for d
in range(2, r
+1):
342 x_powers
.append( Lx
*(x_powers
[-1]) )
344 idmat
= matrix
.identity(R
, n
)
346 W
= self
._charpoly
_basis
_space
()
347 W
= W
.change_ring(R
.fraction_field())
349 # Starting with the standard coordinates x = (X1,X2,...,Xn)
350 # and then converting the entries to W-coordinates allows us
351 # to pass in the standard coordinates to the charpoly and get
352 # back the right answer. Specifically, with x = (X1,X2,...,Xn),
355 # W.coordinates(x^2) eval'd at (standard z-coords)
359 # W-coords of (standard coords of x^2 eval'd at std-coords of z)
361 # We want the middle equivalent thing in our matrix, but use
362 # the first equivalent thing instead so that we can pass in
363 # standard coordinates.
364 x_powers
= [ W
.coordinate_vector(xp
) for xp
in x_powers
]
365 l2
= [idmat
.column(k
-1) for k
in range(r
+1, n
+1)]
366 A_of_x
= matrix
.column(R
, n
, (x_powers
[:r
] + l2
))
367 return (A_of_x
, x
, x_powers
[r
], A_of_x
.det())
371 def characteristic_polynomial(self
):
373 Return a characteristic polynomial that works for all elements
376 The resulting polynomial has `n+1` variables, where `n` is the
377 dimension of this algebra. The first `n` variables correspond to
378 the coordinates of an algebra element: when evaluated at the
379 coordinates of an algebra element with respect to a certain
380 basis, the result is a univariate polynomial (in the one
381 remaining variable ``t``), namely the characteristic polynomial
386 sage: from mjo.eja.eja_algebra import JordanSpinEJA
390 The characteristic polynomial in the spin algebra is given in
391 Alizadeh, Example 11.11::
393 sage: J = JordanSpinEJA(3)
394 sage: p = J.characteristic_polynomial(); p
395 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
396 sage: xvec = J.one().to_vector()
404 # The list of coefficient polynomials a_1, a_2, ..., a_n.
405 a
= [ self
._charpoly
_coeff
(i
) for i
in range(n
) ]
407 # We go to a bit of trouble here to reorder the
408 # indeterminates, so that it's easier to evaluate the
409 # characteristic polynomial at x's coordinates and get back
410 # something in terms of t, which is what we want.
412 S
= PolynomialRing(self
.base_ring(),'t')
414 S
= PolynomialRing(S
, R
.variable_names())
417 # Note: all entries past the rth should be zero. The
418 # coefficient of the highest power (x^r) is 1, but it doesn't
419 # appear in the solution vector which contains coefficients
420 # for the other powers (to make them sum to x^r).
422 a
[r
] = 1 # corresponds to x^r
424 # When the rank is equal to the dimension, trying to
425 # assign a[r] goes out-of-bounds.
426 a
.append(1) # corresponds to x^r
428 return sum( a
[k
]*(t
**k
) for k
in range(len(a
)) )
431 def inner_product(self
, x
, y
):
433 The inner product associated with this Euclidean Jordan algebra.
435 Defaults to the trace inner product, but can be overridden by
436 subclasses if they are sure that the necessary properties are
441 sage: from mjo.eja.eja_algebra import random_eja
445 Our inner product satisfies the Jordan axiom, which is also
446 referred to as "associativity" for a symmetric bilinear form::
448 sage: set_random_seed()
449 sage: J = random_eja()
450 sage: x,y,z = J.random_elements(3)
451 sage: (x*y).inner_product(z) == y.inner_product(x*z)
455 X
= x
.natural_representation()
456 Y
= y
.natural_representation()
457 return self
.natural_inner_product(X
,Y
)
460 def is_trivial(self
):
462 Return whether or not this algebra is trivial.
464 A trivial algebra contains only the zero element.
468 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
472 sage: J = ComplexHermitianEJA(3)
475 sage: A = J.zero().subalgebra_generated_by()
480 return self
.dimension() == 0
483 def multiplication_table(self
):
485 Return a visual representation of this algebra's multiplication
486 table (on basis elements).
490 sage: from mjo.eja.eja_algebra import JordanSpinEJA
494 sage: J = JordanSpinEJA(4)
495 sage: J.multiplication_table()
496 +----++----+----+----+----+
497 | * || e0 | e1 | e2 | e3 |
498 +====++====+====+====+====+
499 | e0 || e0 | e1 | e2 | e3 |
500 +----++----+----+----+----+
501 | e1 || e1 | e0 | 0 | 0 |
502 +----++----+----+----+----+
503 | e2 || e2 | 0 | e0 | 0 |
504 +----++----+----+----+----+
505 | e3 || e3 | 0 | 0 | e0 |
506 +----++----+----+----+----+
509 M
= list(self
._multiplication
_table
) # copy
510 for i
in range(len(M
)):
511 # M had better be "square"
512 M
[i
] = [self
.monomial(i
)] + M
[i
]
513 M
= [["*"] + list(self
.gens())] + M
514 return table(M
, header_row
=True, header_column
=True, frame
=True)
517 def natural_basis(self
):
519 Return a more-natural representation of this algebra's basis.
521 Every finite-dimensional Euclidean Jordan Algebra is a direct
522 sum of five simple algebras, four of which comprise Hermitian
523 matrices. This method returns the original "natural" basis
524 for our underlying vector space. (Typically, the natural basis
525 is used to construct the multiplication table in the first place.)
527 Note that this will always return a matrix. The standard basis
528 in `R^n` will be returned as `n`-by-`1` column matrices.
532 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
533 ....: RealSymmetricEJA)
537 sage: J = RealSymmetricEJA(2)
539 Finite family {0: e0, 1: e1, 2: e2}
540 sage: J.natural_basis()
542 [1 0] [ 0 1/2*sqrt2] [0 0]
543 [0 0], [1/2*sqrt2 0], [0 1]
548 sage: J = JordanSpinEJA(2)
550 Finite family {0: e0, 1: e1}
551 sage: J.natural_basis()
558 if self
._natural
_basis
is None:
559 M
= self
.natural_basis_space()
560 return tuple( M(b
.to_vector()) for b
in self
.basis() )
562 return self
._natural
_basis
565 def natural_basis_space(self
):
567 Return the matrix space in which this algebra's natural basis
570 if self
._natural
_basis
is None or len(self
._natural
_basis
) == 0:
571 return MatrixSpace(self
.base_ring(), self
.dimension(), 1)
573 return self
._natural
_basis
[0].matrix_space()
577 def natural_inner_product(X
,Y
):
579 Compute the inner product of two naturally-represented elements.
581 For example in the real symmetric matrix EJA, this will compute
582 the trace inner-product of two n-by-n symmetric matrices. The
583 default should work for the real cartesian product EJA, the
584 Jordan spin EJA, and the real symmetric matrices. The others
585 will have to be overridden.
587 return (X
.conjugate_transpose()*Y
).trace()
593 Return the unit element of this algebra.
597 sage: from mjo.eja.eja_algebra import (RealCartesianProductEJA,
602 sage: J = RealCartesianProductEJA(5)
604 e0 + e1 + e2 + e3 + e4
608 The identity element acts like the identity::
610 sage: set_random_seed()
611 sage: J = random_eja()
612 sage: x = J.random_element()
613 sage: J.one()*x == x and x*J.one() == x
616 The matrix of the unit element's operator is the identity::
618 sage: set_random_seed()
619 sage: J = random_eja()
620 sage: actual = J.one().operator().matrix()
621 sage: expected = matrix.identity(J.base_ring(), J.dimension())
622 sage: actual == expected
626 # We can brute-force compute the matrices of the operators
627 # that correspond to the basis elements of this algebra.
628 # If some linear combination of those basis elements is the
629 # algebra identity, then the same linear combination of
630 # their matrices has to be the identity matrix.
632 # Of course, matrices aren't vectors in sage, so we have to
633 # appeal to the "long vectors" isometry.
634 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
636 # Now we use basis linear algebra to find the coefficients,
637 # of the matrices-as-vectors-linear-combination, which should
638 # work for the original algebra basis too.
639 A
= matrix
.column(self
.base_ring(), oper_vecs
)
641 # We used the isometry on the left-hand side already, but we
642 # still need to do it for the right-hand side. Recall that we
643 # wanted something that summed to the identity matrix.
644 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
646 # Now if there's an identity element in the algebra, this should work.
647 coeffs
= A
.solve_right(b
)
648 return self
.linear_combination(zip(self
.gens(), coeffs
))
651 def random_element(self
):
652 # Temporary workaround for https://trac.sagemath.org/ticket/28327
653 if self
.is_trivial():
656 s
= super(FiniteDimensionalEuclideanJordanAlgebra
, self
)
657 return s
.random_element()
659 def random_elements(self
, count
):
661 Return ``count`` random elements as a tuple.
665 sage: from mjo.eja.eja_algebra import JordanSpinEJA
669 sage: J = JordanSpinEJA(3)
670 sage: x,y,z = J.random_elements(3)
671 sage: all( [ x in J, y in J, z in J ])
673 sage: len( J.random_elements(10) ) == 10
677 return tuple( self
.random_element() for idx
in xrange(count
) )
680 def random_instance(cls
, field
=QQ
, **kwargs
):
682 Return a random instance of this type of algebra.
684 In subclasses for algebras that we know how to construct, this
685 is a shortcut for constructing test cases and examples.
687 if cls
is FiniteDimensionalEuclideanJordanAlgebra
:
688 # Red flag! But in theory we could do this I guess. The
689 # only finite-dimensional exceptional EJA is the
690 # octononions. So, we could just create an EJA from an
691 # associative matrix algebra (generated by a subset of
692 # elements) with the symmetric product. Or, we could punt
693 # to random_eja() here, override it in our subclasses, and
694 # not worry about it.
695 raise NotImplementedError
697 n
= ZZ
.random_element(1, cls
._max
_test
_case
_size
())
698 return cls(n
, field
, **kwargs
)
703 Return the rank of this EJA.
707 The author knows of no algorithm to compute the rank of an EJA
708 where only the multiplication table is known. In lieu of one, we
709 require the rank to be specified when the algebra is created,
710 and simply pass along that number here.
714 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
715 ....: RealSymmetricEJA,
716 ....: ComplexHermitianEJA,
717 ....: QuaternionHermitianEJA,
722 The rank of the Jordan spin algebra is always two::
724 sage: JordanSpinEJA(2).rank()
726 sage: JordanSpinEJA(3).rank()
728 sage: JordanSpinEJA(4).rank()
731 The rank of the `n`-by-`n` Hermitian real, complex, or
732 quaternion matrices is `n`::
734 sage: RealSymmetricEJA(4).rank()
736 sage: ComplexHermitianEJA(3).rank()
738 sage: QuaternionHermitianEJA(2).rank()
743 Ensure that every EJA that we know how to construct has a
744 positive integer rank::
746 sage: set_random_seed()
747 sage: r = random_eja().rank()
748 sage: r in ZZ and r > 0
755 def vector_space(self
):
757 Return the vector space that underlies this algebra.
761 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
765 sage: J = RealSymmetricEJA(2)
766 sage: J.vector_space()
767 Vector space of dimension 3 over...
770 return self
.zero().to_vector().parent().ambient_vector_space()
773 Element
= FiniteDimensionalEuclideanJordanAlgebraElement
776 class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra
):
778 Return the Euclidean Jordan Algebra corresponding to the set
779 `R^n` under the Hadamard product.
781 Note: this is nothing more than the Cartesian product of ``n``
782 copies of the spin algebra. Once Cartesian product algebras
783 are implemented, this can go.
787 sage: from mjo.eja.eja_algebra import RealCartesianProductEJA
791 This multiplication table can be verified by hand::
793 sage: J = RealCartesianProductEJA(3)
794 sage: e0,e1,e2 = J.gens()
810 We can change the generator prefix::
812 sage: RealCartesianProductEJA(3, prefix='r').gens()
816 def __init__(self
, n
, field
=QQ
, **kwargs
):
817 V
= VectorSpace(field
, n
)
818 mult_table
= [ [ V
.gen(i
)*(i
== j
) for j
in range(n
) ]
821 fdeja
= super(RealCartesianProductEJA
, self
)
822 return fdeja
.__init
__(field
, mult_table
, rank
=n
, **kwargs
)
824 def inner_product(self
, x
, y
):
826 Faster to reimplement than to use natural representations.
830 sage: from mjo.eja.eja_algebra import RealCartesianProductEJA
834 Ensure that this is the usual inner product for the algebras
837 sage: set_random_seed()
838 sage: J = RealCartesianProductEJA.random_instance()
839 sage: x,y = J.random_elements(2)
840 sage: X = x.natural_representation()
841 sage: Y = y.natural_representation()
842 sage: x.inner_product(y) == J.natural_inner_product(X,Y)
846 return x
.to_vector().inner_product(y
.to_vector())
851 Return a "random" finite-dimensional Euclidean Jordan Algebra.
855 For now, we choose a random natural number ``n`` (greater than zero)
856 and then give you back one of the following:
858 * The cartesian product of the rational numbers ``n`` times; this is
859 ``QQ^n`` with the Hadamard product.
861 * The Jordan spin algebra on ``QQ^n``.
863 * The ``n``-by-``n`` rational symmetric matrices with the symmetric
866 * The ``n``-by-``n`` complex-rational Hermitian matrices embedded
867 in the space of ``2n``-by-``2n`` real symmetric matrices.
869 * The ``n``-by-``n`` quaternion-rational Hermitian matrices embedded
870 in the space of ``4n``-by-``4n`` real symmetric matrices.
872 Later this might be extended to return Cartesian products of the
877 sage: from mjo.eja.eja_algebra import random_eja
882 Euclidean Jordan algebra of dimension...
885 classname
= choice([RealCartesianProductEJA
,
889 QuaternionHermitianEJA
])
890 return classname
.random_instance()
897 class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra
):
899 def _max_test_case_size():
900 # Play it safe, since this will be squared and the underlying
901 # field can have dimension 4 (quaternions) too.
905 def _denormalized_basis(cls
, n
, field
):
906 raise NotImplementedError
908 def __init__(self
, n
, field
=QQ
, normalize_basis
=True, **kwargs
):
909 S
= self
._denormalized
_basis
(n
, field
)
911 if n
> 1 and normalize_basis
:
912 # We'll need sqrt(2) to normalize the basis, and this
913 # winds up in the multiplication table, so the whole
914 # algebra needs to be over the field extension.
915 R
= PolynomialRing(field
, 'z')
918 if p
.is_irreducible():
919 field
= NumberField(p
, 'sqrt2', embedding
=RLF(2).sqrt())
920 S
= [ s
.change_ring(field
) for s
in S
]
921 self
._basis
_normalizers
= tuple(
922 ~
(self
.natural_inner_product(s
,s
).sqrt()) for s
in S
)
923 S
= tuple( s
*c
for (s
,c
) in zip(S
,self
._basis
_normalizers
) )
925 Qs
= self
.multiplication_table_from_matrix_basis(S
)
927 fdeja
= super(MatrixEuclideanJordanAlgebra
, self
)
928 return fdeja
.__init
__(field
,
936 def multiplication_table_from_matrix_basis(basis
):
938 At least three of the five simple Euclidean Jordan algebras have the
939 symmetric multiplication (A,B) |-> (AB + BA)/2, where the
940 multiplication on the right is matrix multiplication. Given a basis
941 for the underlying matrix space, this function returns a
942 multiplication table (obtained by looping through the basis
943 elements) for an algebra of those matrices.
945 # In S^2, for example, we nominally have four coordinates even
946 # though the space is of dimension three only. The vector space V
947 # is supposed to hold the entire long vector, and the subspace W
948 # of V will be spanned by the vectors that arise from symmetric
949 # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
950 field
= basis
[0].base_ring()
951 dimension
= basis
[0].nrows()
953 V
= VectorSpace(field
, dimension
**2)
954 W
= V
.span_of_basis( _mat2vec(s
) for s
in basis
)
956 mult_table
= [[W
.zero() for j
in range(n
)] for i
in range(n
)]
959 mat_entry
= (basis
[i
]*basis
[j
] + basis
[j
]*basis
[i
])/2
960 mult_table
[i
][j
] = W
.coordinate_vector(_mat2vec(mat_entry
))
968 Embed the matrix ``M`` into a space of real matrices.
970 The matrix ``M`` can have entries in any field at the moment:
971 the real numbers, complex numbers, or quaternions. And although
972 they are not a field, we can probably support octonions at some
973 point, too. This function returns a real matrix that "acts like"
974 the original with respect to matrix multiplication; i.e.
976 real_embed(M*N) = real_embed(M)*real_embed(N)
979 raise NotImplementedError
985 The inverse of :meth:`real_embed`.
987 raise NotImplementedError
991 def natural_inner_product(cls
,X
,Y
):
992 Xu
= cls
.real_unembed(X
)
993 Yu
= cls
.real_unembed(Y
)
999 # Otherwise, try the thing that works for complex numbers; and
1000 # if that doesn't work, the thing that works for quaternions.
1002 return tr
.vector()[0] # real part, imag part is index 1
1003 except AttributeError:
1004 # A quaternions doesn't have a vector() method, but does
1005 # have coefficient_tuple() method that returns the
1006 # coefficients of 1, i, j, and k -- in that order.
1007 return tr
.coefficient_tuple()[0]
1010 class RealMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1014 Embed the matrix ``M`` into a space of real matrices.
1016 The matrix ``M`` can have entries in any field at the moment:
1017 the real numbers, complex numbers, or quaternions. And although
1018 they are not a field, we can probably support octonions at some
1019 point, too. This function returns a real matrix that "acts like"
1020 the original with respect to matrix multiplication; i.e.
1022 real_embed(M*N) = real_embed(M)*real_embed(N)
1029 def real_unembed(M
):
1031 The inverse of :meth:`real_embed`.
1036 class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra
):
1038 The rank-n simple EJA consisting of real symmetric n-by-n
1039 matrices, the usual symmetric Jordan product, and the trace inner
1040 product. It has dimension `(n^2 + n)/2` over the reals.
1044 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1048 sage: J = RealSymmetricEJA(2)
1049 sage: e0, e1, e2 = J.gens()
1059 The dimension of this algebra is `(n^2 + n) / 2`::
1061 sage: set_random_seed()
1062 sage: n_max = RealSymmetricEJA._max_test_case_size()
1063 sage: n = ZZ.random_element(1, n_max)
1064 sage: J = RealSymmetricEJA(n)
1065 sage: J.dimension() == (n^2 + n)/2
1068 The Jordan multiplication is what we think it is::
1070 sage: set_random_seed()
1071 sage: J = RealSymmetricEJA.random_instance()
1072 sage: x,y = J.random_elements(2)
1073 sage: actual = (x*y).natural_representation()
1074 sage: X = x.natural_representation()
1075 sage: Y = y.natural_representation()
1076 sage: expected = (X*Y + Y*X)/2
1077 sage: actual == expected
1079 sage: J(expected) == x*y
1082 We can change the generator prefix::
1084 sage: RealSymmetricEJA(3, prefix='q').gens()
1085 (q0, q1, q2, q3, q4, q5)
1087 Our natural basis is normalized with respect to the natural inner
1088 product unless we specify otherwise::
1090 sage: set_random_seed()
1091 sage: J = RealSymmetricEJA.random_instance()
1092 sage: all( b.norm() == 1 for b in J.gens() )
1095 Since our natural basis is normalized with respect to the natural
1096 inner product, and since we know that this algebra is an EJA, any
1097 left-multiplication operator's matrix will be symmetric because
1098 natural->EJA basis representation is an isometry and within the EJA
1099 the operator is self-adjoint by the Jordan axiom::
1101 sage: set_random_seed()
1102 sage: x = RealSymmetricEJA.random_instance().random_element()
1103 sage: x.operator().matrix().is_symmetric()
1108 def _denormalized_basis(cls
, n
, field
):
1110 Return a basis for the space of real symmetric n-by-n matrices.
1114 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1118 sage: set_random_seed()
1119 sage: n = ZZ.random_element(1,5)
1120 sage: B = RealSymmetricEJA._denormalized_basis(n,QQ)
1121 sage: all( M.is_symmetric() for M in B)
1125 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1129 for j
in xrange(i
+1):
1130 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
1134 Sij
= Eij
+ Eij
.transpose()
1140 def _max_test_case_size():
1141 return 5 # Dimension 10
1145 class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1149 Embed the n-by-n complex matrix ``M`` into the space of real
1150 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1151 bi` to the block matrix ``[[a,b],[-b,a]]``.
1155 sage: from mjo.eja.eja_algebra import \
1156 ....: ComplexMatrixEuclideanJordanAlgebra
1160 sage: F = QuadraticField(-1, 'i')
1161 sage: x1 = F(4 - 2*i)
1162 sage: x2 = F(1 + 2*i)
1165 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1166 sage: ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1175 Embedding is a homomorphism (isomorphism, in fact)::
1177 sage: set_random_seed()
1178 sage: n_max = ComplexMatrixEuclideanJordanAlgebra._max_test_case_size()
1179 sage: n = ZZ.random_element(n_max)
1180 sage: F = QuadraticField(-1, 'i')
1181 sage: X = random_matrix(F, n)
1182 sage: Y = random_matrix(F, n)
1183 sage: Xe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X)
1184 sage: Ye = ComplexMatrixEuclideanJordanAlgebra.real_embed(Y)
1185 sage: XYe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1192 raise ValueError("the matrix 'M' must be square")
1193 field
= M
.base_ring()
1196 a
= z
.vector()[0] # real part, I guess
1197 b
= z
.vector()[1] # imag part, I guess
1198 blocks
.append(matrix(field
, 2, [[a
,b
],[-b
,a
]]))
1200 # We can drop the imaginaries here.
1201 return matrix
.block(field
.base_ring(), n
, blocks
)
1205 def real_unembed(M
):
1207 The inverse of _embed_complex_matrix().
1211 sage: from mjo.eja.eja_algebra import \
1212 ....: ComplexMatrixEuclideanJordanAlgebra
1216 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1217 ....: [-2, 1, -4, 3],
1218 ....: [ 9, 10, 11, 12],
1219 ....: [-10, 9, -12, 11] ])
1220 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(A)
1222 [ 10*i + 9 12*i + 11]
1226 Unembedding is the inverse of embedding::
1228 sage: set_random_seed()
1229 sage: F = QuadraticField(-1, 'i')
1230 sage: M = random_matrix(F, 3)
1231 sage: Me = ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1232 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1238 raise ValueError("the matrix 'M' must be square")
1239 if not n
.mod(2).is_zero():
1240 raise ValueError("the matrix 'M' must be a complex embedding")
1242 field
= M
.base_ring() # This should already have sqrt2
1243 R
= PolynomialRing(field
, 'z')
1245 F
= NumberField(z
**2 + 1,'i', embedding
=CLF(-1).sqrt())
1248 # Go top-left to bottom-right (reading order), converting every
1249 # 2-by-2 block we see to a single complex element.
1251 for k
in xrange(n
/2):
1252 for j
in xrange(n
/2):
1253 submat
= M
[2*k
:2*k
+2,2*j
:2*j
+2]
1254 if submat
[0,0] != submat
[1,1]:
1255 raise ValueError('bad on-diagonal submatrix')
1256 if submat
[0,1] != -submat
[1,0]:
1257 raise ValueError('bad off-diagonal submatrix')
1258 z
= submat
[0,0] + submat
[0,1]*i
1261 return matrix(F
, n
/2, elements
)
1265 def natural_inner_product(cls
,X
,Y
):
1267 Compute a natural inner product in this algebra directly from
1272 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1276 This gives the same answer as the slow, default method implemented
1277 in :class:`MatrixEuclideanJordanAlgebra`::
1279 sage: set_random_seed()
1280 sage: J = ComplexHermitianEJA.random_instance()
1281 sage: x,y = J.random_elements(2)
1282 sage: Xe = x.natural_representation()
1283 sage: Ye = y.natural_representation()
1284 sage: X = ComplexHermitianEJA.real_unembed(Xe)
1285 sage: Y = ComplexHermitianEJA.real_unembed(Ye)
1286 sage: expected = (X*Y).trace().vector()[0]
1287 sage: actual = ComplexHermitianEJA.natural_inner_product(Xe,Ye)
1288 sage: actual == expected
1292 return RealMatrixEuclideanJordanAlgebra
.natural_inner_product(X
,Y
)/2
1295 class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra
):
1297 The rank-n simple EJA consisting of complex Hermitian n-by-n
1298 matrices over the real numbers, the usual symmetric Jordan product,
1299 and the real-part-of-trace inner product. It has dimension `n^2` over
1304 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1308 The dimension of this algebra is `n^2`::
1310 sage: set_random_seed()
1311 sage: n_max = ComplexHermitianEJA._max_test_case_size()
1312 sage: n = ZZ.random_element(1, n_max)
1313 sage: J = ComplexHermitianEJA(n)
1314 sage: J.dimension() == n^2
1317 The Jordan multiplication is what we think it is::
1319 sage: set_random_seed()
1320 sage: J = ComplexHermitianEJA.random_instance()
1321 sage: x,y = J.random_elements(2)
1322 sage: actual = (x*y).natural_representation()
1323 sage: X = x.natural_representation()
1324 sage: Y = y.natural_representation()
1325 sage: expected = (X*Y + Y*X)/2
1326 sage: actual == expected
1328 sage: J(expected) == x*y
1331 We can change the generator prefix::
1333 sage: ComplexHermitianEJA(2, prefix='z').gens()
1336 Our natural basis is normalized with respect to the natural inner
1337 product unless we specify otherwise::
1339 sage: set_random_seed()
1340 sage: J = ComplexHermitianEJA.random_instance()
1341 sage: all( b.norm() == 1 for b in J.gens() )
1344 Since our natural basis is normalized with respect to the natural
1345 inner product, and since we know that this algebra is an EJA, any
1346 left-multiplication operator's matrix will be symmetric because
1347 natural->EJA basis representation is an isometry and within the EJA
1348 the operator is self-adjoint by the Jordan axiom::
1350 sage: set_random_seed()
1351 sage: x = ComplexHermitianEJA.random_instance().random_element()
1352 sage: x.operator().matrix().is_symmetric()
1357 def _denormalized_basis(cls
, n
, field
):
1359 Returns a basis for the space of complex Hermitian n-by-n matrices.
1361 Why do we embed these? Basically, because all of numerical linear
1362 algebra assumes that you're working with vectors consisting of `n`
1363 entries from a field and scalars from the same field. There's no way
1364 to tell SageMath that (for example) the vectors contain complex
1365 numbers, while the scalar field is real.
1369 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1373 sage: set_random_seed()
1374 sage: n = ZZ.random_element(1,5)
1375 sage: field = QuadraticField(2, 'sqrt2')
1376 sage: B = ComplexHermitianEJA._denormalized_basis(n, field)
1377 sage: all( M.is_symmetric() for M in B)
1381 R
= PolynomialRing(field
, 'z')
1383 F
= NumberField(z
**2 + 1, 'I', embedding
=CLF(-1).sqrt())
1386 # This is like the symmetric case, but we need to be careful:
1388 # * We want conjugate-symmetry, not just symmetry.
1389 # * The diagonal will (as a result) be real.
1393 for j
in xrange(i
+1):
1394 Eij
= matrix(F
, n
, lambda k
,l
: k
==i
and l
==j
)
1396 Sij
= cls
.real_embed(Eij
)
1399 # The second one has a minus because it's conjugated.
1400 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
1402 Sij_imag
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
1405 # Since we embedded these, we can drop back to the "field" that we
1406 # started with instead of the complex extension "F".
1407 return tuple( s
.change_ring(field
) for s
in S
)
1411 class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1415 Embed the n-by-n quaternion matrix ``M`` into the space of real
1416 matrices of size 4n-by-4n by first sending each quaternion entry `z
1417 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
1418 c+di],[-c + di, a-bi]]`, and then embedding those into a real
1423 sage: from mjo.eja.eja_algebra import \
1424 ....: QuaternionMatrixEuclideanJordanAlgebra
1428 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1429 sage: i,j,k = Q.gens()
1430 sage: x = 1 + 2*i + 3*j + 4*k
1431 sage: M = matrix(Q, 1, [[x]])
1432 sage: QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1438 Embedding is a homomorphism (isomorphism, in fact)::
1440 sage: set_random_seed()
1441 sage: n_max = QuaternionMatrixEuclideanJordanAlgebra._max_test_case_size()
1442 sage: n = ZZ.random_element(n_max)
1443 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1444 sage: X = random_matrix(Q, n)
1445 sage: Y = random_matrix(Q, n)
1446 sage: Xe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X)
1447 sage: Ye = QuaternionMatrixEuclideanJordanAlgebra.real_embed(Y)
1448 sage: XYe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1453 quaternions
= M
.base_ring()
1456 raise ValueError("the matrix 'M' must be square")
1458 F
= QuadraticField(-1, 'i')
1463 t
= z
.coefficient_tuple()
1468 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
1469 [-c
+ d
*i
, a
- b
*i
]])
1470 realM
= ComplexMatrixEuclideanJordanAlgebra
.real_embed(cplxM
)
1471 blocks
.append(realM
)
1473 # We should have real entries by now, so use the realest field
1474 # we've got for the return value.
1475 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
1480 def real_unembed(M
):
1482 The inverse of _embed_quaternion_matrix().
1486 sage: from mjo.eja.eja_algebra import \
1487 ....: QuaternionMatrixEuclideanJordanAlgebra
1491 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
1492 ....: [-2, 1, -4, 3],
1493 ....: [-3, 4, 1, -2],
1494 ....: [-4, -3, 2, 1]])
1495 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(M)
1496 [1 + 2*i + 3*j + 4*k]
1500 Unembedding is the inverse of embedding::
1502 sage: set_random_seed()
1503 sage: Q = QuaternionAlgebra(QQ, -1, -1)
1504 sage: M = random_matrix(Q, 3)
1505 sage: Me = QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1506 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1512 raise ValueError("the matrix 'M' must be square")
1513 if not n
.mod(4).is_zero():
1514 raise ValueError("the matrix 'M' must be a complex embedding")
1516 # Use the base ring of the matrix to ensure that its entries can be
1517 # multiplied by elements of the quaternion algebra.
1518 field
= M
.base_ring()
1519 Q
= QuaternionAlgebra(field
,-1,-1)
1522 # Go top-left to bottom-right (reading order), converting every
1523 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
1526 for l
in xrange(n
/4):
1527 for m
in xrange(n
/4):
1528 submat
= ComplexMatrixEuclideanJordanAlgebra
.real_unembed(
1529 M
[4*l
:4*l
+4,4*m
:4*m
+4] )
1530 if submat
[0,0] != submat
[1,1].conjugate():
1531 raise ValueError('bad on-diagonal submatrix')
1532 if submat
[0,1] != -submat
[1,0].conjugate():
1533 raise ValueError('bad off-diagonal submatrix')
1534 z
= submat
[0,0].vector()[0] # real part
1535 z
+= submat
[0,0].vector()[1]*i
# imag part
1536 z
+= submat
[0,1].vector()[0]*j
# real part
1537 z
+= submat
[0,1].vector()[1]*k
# imag part
1540 return matrix(Q
, n
/4, elements
)
1544 def natural_inner_product(cls
,X
,Y
):
1546 Compute a natural inner product in this algebra directly from
1551 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1555 This gives the same answer as the slow, default method implemented
1556 in :class:`MatrixEuclideanJordanAlgebra`::
1558 sage: set_random_seed()
1559 sage: J = QuaternionHermitianEJA.random_instance()
1560 sage: x,y = J.random_elements(2)
1561 sage: Xe = x.natural_representation()
1562 sage: Ye = y.natural_representation()
1563 sage: X = QuaternionHermitianEJA.real_unembed(Xe)
1564 sage: Y = QuaternionHermitianEJA.real_unembed(Ye)
1565 sage: expected = (X*Y).trace().coefficient_tuple()[0]
1566 sage: actual = QuaternionHermitianEJA.natural_inner_product(Xe,Ye)
1567 sage: actual == expected
1571 return RealMatrixEuclideanJordanAlgebra
.natural_inner_product(X
,Y
)/4
1574 class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra
):
1576 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
1577 matrices, the usual symmetric Jordan product, and the
1578 real-part-of-trace inner product. It has dimension `2n^2 - n` over
1583 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1587 The dimension of this algebra is `2*n^2 - n`::
1589 sage: set_random_seed()
1590 sage: n_max = QuaternionHermitianEJA._max_test_case_size()
1591 sage: n = ZZ.random_element(1, n_max)
1592 sage: J = QuaternionHermitianEJA(n)
1593 sage: J.dimension() == 2*(n^2) - n
1596 The Jordan multiplication is what we think it is::
1598 sage: set_random_seed()
1599 sage: J = QuaternionHermitianEJA.random_instance()
1600 sage: x,y = J.random_elements(2)
1601 sage: actual = (x*y).natural_representation()
1602 sage: X = x.natural_representation()
1603 sage: Y = y.natural_representation()
1604 sage: expected = (X*Y + Y*X)/2
1605 sage: actual == expected
1607 sage: J(expected) == x*y
1610 We can change the generator prefix::
1612 sage: QuaternionHermitianEJA(2, prefix='a').gens()
1613 (a0, a1, a2, a3, a4, a5)
1615 Our natural basis is normalized with respect to the natural inner
1616 product unless we specify otherwise::
1618 sage: set_random_seed()
1619 sage: J = QuaternionHermitianEJA.random_instance()
1620 sage: all( b.norm() == 1 for b in J.gens() )
1623 Since our natural basis is normalized with respect to the natural
1624 inner product, and since we know that this algebra is an EJA, any
1625 left-multiplication operator's matrix will be symmetric because
1626 natural->EJA basis representation is an isometry and within the EJA
1627 the operator is self-adjoint by the Jordan axiom::
1629 sage: set_random_seed()
1630 sage: x = QuaternionHermitianEJA.random_instance().random_element()
1631 sage: x.operator().matrix().is_symmetric()
1636 def _denormalized_basis(cls
, n
, field
):
1638 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
1640 Why do we embed these? Basically, because all of numerical
1641 linear algebra assumes that you're working with vectors consisting
1642 of `n` entries from a field and scalars from the same field. There's
1643 no way to tell SageMath that (for example) the vectors contain
1644 complex numbers, while the scalar field is real.
1648 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1652 sage: set_random_seed()
1653 sage: n = ZZ.random_element(1,5)
1654 sage: B = QuaternionHermitianEJA._denormalized_basis(n,QQ)
1655 sage: all( M.is_symmetric() for M in B )
1659 Q
= QuaternionAlgebra(QQ
,-1,-1)
1662 # This is like the symmetric case, but we need to be careful:
1664 # * We want conjugate-symmetry, not just symmetry.
1665 # * The diagonal will (as a result) be real.
1669 for j
in xrange(i
+1):
1670 Eij
= matrix(Q
, n
, lambda k
,l
: k
==i
and l
==j
)
1672 Sij
= cls
.real_embed(Eij
)
1675 # The second, third, and fourth ones have a minus
1676 # because they're conjugated.
1677 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
1679 Sij_I
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
1681 Sij_J
= cls
.real_embed(J
*Eij
- J
*Eij
.transpose())
1683 Sij_K
= cls
.real_embed(K
*Eij
- K
*Eij
.transpose())
1689 class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1691 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
1692 with the usual inner product and jordan product ``x*y =
1693 (<x_bar,y_bar>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
1698 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1702 This multiplication table can be verified by hand::
1704 sage: J = JordanSpinEJA(4)
1705 sage: e0,e1,e2,e3 = J.gens()
1721 We can change the generator prefix::
1723 sage: JordanSpinEJA(2, prefix='B').gens()
1727 def __init__(self
, n
, field
=QQ
, **kwargs
):
1728 V
= VectorSpace(field
, n
)
1729 mult_table
= [[V
.zero() for j
in range(n
)] for i
in range(n
)]
1739 z0
= x
.inner_product(y
)
1740 zbar
= y0
*xbar
+ x0
*ybar
1741 z
= V([z0
] + zbar
.list())
1742 mult_table
[i
][j
] = z
1744 # The rank of the spin algebra is two, unless we're in a
1745 # one-dimensional ambient space (because the rank is bounded by
1746 # the ambient dimension).
1747 fdeja
= super(JordanSpinEJA
, self
)
1748 return fdeja
.__init
__(field
, mult_table
, rank
=min(n
,2), **kwargs
)
1750 def inner_product(self
, x
, y
):
1752 Faster to reimplement than to use natural representations.
1756 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1760 Ensure that this is the usual inner product for the algebras
1763 sage: set_random_seed()
1764 sage: J = JordanSpinEJA.random_instance()
1765 sage: x,y = J.random_elements(2)
1766 sage: X = x.natural_representation()
1767 sage: Y = y.natural_representation()
1768 sage: x.inner_product(y) == J.natural_inner_product(X,Y)
1772 return x
.to_vector().inner_product(y
.to_vector())