2 Euclidean Jordan Algebras. These are formally-real Jordan Algebras;
3 specifically those where u^2 + v^2 = 0 implies that u = v = 0. They
4 are used in optimization, and have some additional nice methods beyond
5 what can be supported in a general Jordan Algebra.
8 from itertools
import izip
, repeat
10 from sage
.algebras
.quatalg
.quaternion_algebra
import QuaternionAlgebra
11 from sage
.categories
.magmatic_algebras
import MagmaticAlgebras
12 from sage
.combinat
.free_module
import CombinatorialFreeModule
13 from sage
.matrix
.constructor
import matrix
14 from sage
.matrix
.matrix_space
import MatrixSpace
15 from sage
.misc
.cachefunc
import cached_method
16 from sage
.misc
.prandom
import choice
17 from sage
.misc
.table
import table
18 from sage
.modules
.free_module
import FreeModule
, VectorSpace
19 from sage
.rings
.integer_ring
import ZZ
20 from sage
.rings
.number_field
.number_field
import NumberField
, QuadraticField
21 from sage
.rings
.polynomial
.polynomial_ring_constructor
import PolynomialRing
22 from sage
.rings
.rational_field
import QQ
23 from sage
.rings
.real_lazy
import CLF
, RLF
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
65 category
= MagmaticAlgebras(field
).FiniteDimensional()
66 category
= category
.WithBasis().Unital()
68 fda
= super(FiniteDimensionalEuclideanJordanAlgebra
, self
)
70 range(len(mult_table
)),
73 self
.print_options(bracket
='')
75 # The multiplication table we're given is necessarily in terms
76 # of vectors, because we don't have an algebra yet for
77 # anything to be an element of. However, it's faster in the
78 # long run to have the multiplication table be in terms of
79 # algebra elements. We do this after calling the superclass
80 # constructor so that from_vector() knows what to do.
81 self
._multiplication
_table
= [ map(lambda x
: self
.from_vector(x
), ls
)
82 for ls
in mult_table
]
85 def _element_constructor_(self
, elt
):
87 Construct an element of this algebra from its natural
90 This gets called only after the parent element _call_ method
91 fails to find a coercion for the argument.
95 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
96 ....: RealCartesianProductEJA,
97 ....: RealSymmetricEJA)
101 The identity in `S^n` is converted to the identity in the EJA::
103 sage: J = RealSymmetricEJA(3)
104 sage: I = matrix.identity(QQ,3)
105 sage: J(I) == J.one()
108 This skew-symmetric matrix can't be represented in the EJA::
110 sage: J = RealSymmetricEJA(3)
111 sage: A = matrix(QQ,3, lambda i,j: i-j)
113 Traceback (most recent call last):
115 ArithmeticError: vector is not in free module
119 Ensure that we can convert any element of the two non-matrix
120 simple algebras (whose natural representations are their usual
121 vector representations) back and forth faithfully::
123 sage: set_random_seed()
124 sage: J = RealCartesianProductEJA.random_instance()
125 sage: x = J.random_element()
126 sage: J(x.to_vector().column()) == x
128 sage: J = JordanSpinEJA.random_instance()
129 sage: x = J.random_element()
130 sage: J(x.to_vector().column()) == x
135 # The superclass implementation of random_element()
136 # needs to be able to coerce "0" into the algebra.
139 natural_basis
= self
.natural_basis()
140 basis_space
= natural_basis
[0].matrix_space()
141 if elt
not in basis_space
:
142 raise ValueError("not a naturally-represented algebra element")
144 # Thanks for nothing! Matrix spaces aren't vector spaces in
145 # Sage, so we have to figure out its natural-basis coordinates
146 # ourselves. We use the basis space's ring instead of the
147 # element's ring because the basis space might be an algebraic
148 # closure whereas the base ring of the 3-by-3 identity matrix
149 # could be QQ instead of QQbar.
150 V
= VectorSpace(basis_space
.base_ring(), elt
.nrows()*elt
.ncols())
151 W
= V
.span_of_basis( _mat2vec(s
) for s
in natural_basis
)
152 coords
= W
.coordinate_vector(_mat2vec(elt
))
153 return self
.from_vector(coords
)
158 Return a string representation of ``self``.
162 sage: from mjo.eja.eja_algebra import JordanSpinEJA
166 Ensure that it says what we think it says::
168 sage: JordanSpinEJA(2, field=QQ)
169 Euclidean Jordan algebra of dimension 2 over Rational Field
170 sage: JordanSpinEJA(3, field=RDF)
171 Euclidean Jordan algebra of dimension 3 over Real Double Field
174 fmt
= "Euclidean Jordan algebra of dimension {} over {}"
175 return fmt
.format(self
.dimension(), self
.base_ring())
177 def product_on_basis(self
, i
, j
):
178 return self
._multiplication
_table
[i
][j
]
180 def _a_regular_element(self
):
182 Guess a regular element. Needed to compute the basis for our
183 characteristic polynomial coefficients.
187 sage: from mjo.eja.eja_algebra import random_eja
191 Ensure that this hacky method succeeds for every algebra that we
192 know how to construct::
194 sage: set_random_seed()
195 sage: J = random_eja()
196 sage: J._a_regular_element().is_regular()
201 z
= self
.sum( (i
+1)*gs
[i
] for i
in range(len(gs
)) )
202 if not z
.is_regular():
203 raise ValueError("don't know a regular element")
208 def _charpoly_basis_space(self
):
210 Return the vector space spanned by the basis used in our
211 characteristic polynomial coefficients. This is used not only to
212 compute those coefficients, but also any time we need to
213 evaluate the coefficients (like when we compute the trace or
216 z
= self
._a
_regular
_element
()
217 # Don't use the parent vector space directly here in case this
218 # happens to be a subalgebra. In that case, we would be e.g.
219 # two-dimensional but span_of_basis() would expect three
221 V
= VectorSpace(self
.base_ring(), self
.vector_space().dimension())
222 basis
= [ (z
**k
).to_vector() for k
in range(self
.rank()) ]
223 V1
= V
.span_of_basis( basis
)
224 b
= (V1
.basis() + V1
.complement().basis())
225 return V
.span_of_basis(b
)
230 def _charpoly_coeff(self
, i
):
232 Return the coefficient polynomial "a_{i}" of this algebra's
233 general characteristic polynomial.
235 Having this be a separate cached method lets us compute and
236 store the trace/determinant (a_{r-1} and a_{0} respectively)
237 separate from the entire characteristic polynomial.
239 (A_of_x
, x
, xr
, detA
) = self
._charpoly
_matrix
_system
()
240 R
= A_of_x
.base_ring()
242 # Guaranteed by theory
245 # Danger: the in-place modification is done for performance
246 # reasons (reconstructing a matrix with huge polynomial
247 # entries is slow), but I don't know how cached_method works,
248 # so it's highly possible that we're modifying some global
249 # list variable by reference, here. In other words, you
250 # probably shouldn't call this method twice on the same
251 # algebra, at the same time, in two threads
252 Ai_orig
= A_of_x
.column(i
)
253 A_of_x
.set_column(i
,xr
)
254 numerator
= A_of_x
.det()
255 A_of_x
.set_column(i
,Ai_orig
)
257 # We're relying on the theory here to ensure that each a_i is
258 # indeed back in R, and the added negative signs are to make
259 # the whole charpoly expression sum to zero.
260 return R(-numerator
/detA
)
264 def _charpoly_matrix_system(self
):
266 Compute the matrix whose entries A_ij are polynomials in
267 X1,...,XN, the vector ``x`` of variables X1,...,XN, the vector
268 corresponding to `x^r` and the determinent of the matrix A =
269 [A_ij]. In other words, all of the fixed (cachable) data needed
270 to compute the coefficients of the characteristic polynomial.
275 # Turn my vector space into a module so that "vectors" can
276 # have multivatiate polynomial entries.
277 names
= tuple('X' + str(i
) for i
in range(1,n
+1))
278 R
= PolynomialRing(self
.base_ring(), names
)
280 # Using change_ring() on the parent's vector space doesn't work
281 # here because, in a subalgebra, that vector space has a basis
282 # and change_ring() tries to bring the basis along with it. And
283 # that doesn't work unless the new ring is a PID, which it usually
287 # Now let x = (X1,X2,...,Xn) be the vector whose entries are
291 # And figure out the "left multiplication by x" matrix in
294 monomial_matrices
= [ self
.monomial(i
).operator().matrix()
295 for i
in range(n
) ] # don't recompute these!
297 ek
= self
.monomial(k
).to_vector()
299 sum( x
[i
]*(monomial_matrices
[i
]*ek
)
300 for i
in range(n
) ) )
301 Lx
= matrix
.column(R
, lmbx_cols
)
303 # Now we can compute powers of x "symbolically"
304 x_powers
= [self
.one().to_vector(), x
]
305 for d
in range(2, r
+1):
306 x_powers
.append( Lx
*(x_powers
[-1]) )
308 idmat
= matrix
.identity(R
, n
)
310 W
= self
._charpoly
_basis
_space
()
311 W
= W
.change_ring(R
.fraction_field())
313 # Starting with the standard coordinates x = (X1,X2,...,Xn)
314 # and then converting the entries to W-coordinates allows us
315 # to pass in the standard coordinates to the charpoly and get
316 # back the right answer. Specifically, with x = (X1,X2,...,Xn),
319 # W.coordinates(x^2) eval'd at (standard z-coords)
323 # W-coords of (standard coords of x^2 eval'd at std-coords of z)
325 # We want the middle equivalent thing in our matrix, but use
326 # the first equivalent thing instead so that we can pass in
327 # standard coordinates.
328 x_powers
= [ W
.coordinate_vector(xp
) for xp
in x_powers
]
329 l2
= [idmat
.column(k
-1) for k
in range(r
+1, n
+1)]
330 A_of_x
= matrix
.column(R
, n
, (x_powers
[:r
] + l2
))
331 return (A_of_x
, x
, x_powers
[r
], A_of_x
.det())
335 def characteristic_polynomial(self
):
337 Return a characteristic polynomial that works for all elements
340 The resulting polynomial has `n+1` variables, where `n` is the
341 dimension of this algebra. The first `n` variables correspond to
342 the coordinates of an algebra element: when evaluated at the
343 coordinates of an algebra element with respect to a certain
344 basis, the result is a univariate polynomial (in the one
345 remaining variable ``t``), namely the characteristic polynomial
350 sage: from mjo.eja.eja_algebra import JordanSpinEJA
354 The characteristic polynomial in the spin algebra is given in
355 Alizadeh, Example 11.11::
357 sage: J = JordanSpinEJA(3)
358 sage: p = J.characteristic_polynomial(); p
359 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
360 sage: xvec = J.one().to_vector()
368 # The list of coefficient polynomials a_1, a_2, ..., a_n.
369 a
= [ self
._charpoly
_coeff
(i
) for i
in range(n
) ]
371 # We go to a bit of trouble here to reorder the
372 # indeterminates, so that it's easier to evaluate the
373 # characteristic polynomial at x's coordinates and get back
374 # something in terms of t, which is what we want.
376 S
= PolynomialRing(self
.base_ring(),'t')
378 S
= PolynomialRing(S
, R
.variable_names())
381 # Note: all entries past the rth should be zero. The
382 # coefficient of the highest power (x^r) is 1, but it doesn't
383 # appear in the solution vector which contains coefficients
384 # for the other powers (to make them sum to x^r).
386 a
[r
] = 1 # corresponds to x^r
388 # When the rank is equal to the dimension, trying to
389 # assign a[r] goes out-of-bounds.
390 a
.append(1) # corresponds to x^r
392 return sum( a
[k
]*(t
**k
) for k
in xrange(len(a
)) )
395 def inner_product(self
, x
, y
):
397 The inner product associated with this Euclidean Jordan algebra.
399 Defaults to the trace inner product, but can be overridden by
400 subclasses if they are sure that the necessary properties are
405 sage: from mjo.eja.eja_algebra import random_eja
409 Our inner product satisfies the Jordan axiom, which is also
410 referred to as "associativity" for a symmetric bilinear form::
412 sage: set_random_seed()
413 sage: J = random_eja()
414 sage: x,y,z = J.random_elements(3)
415 sage: (x*y).inner_product(z) == y.inner_product(x*z)
419 X
= x
.natural_representation()
420 Y
= y
.natural_representation()
421 return self
.natural_inner_product(X
,Y
)
424 def is_trivial(self
):
426 Return whether or not this algebra is trivial.
428 A trivial algebra contains only the zero element.
432 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
436 sage: J = ComplexHermitianEJA(3)
439 sage: A = J.zero().subalgebra_generated_by()
444 return self
.dimension() == 0
447 def multiplication_table(self
):
449 Return a visual representation of this algebra's multiplication
450 table (on basis elements).
454 sage: from mjo.eja.eja_algebra import JordanSpinEJA
458 sage: J = JordanSpinEJA(4)
459 sage: J.multiplication_table()
460 +----++----+----+----+----+
461 | * || e0 | e1 | e2 | e3 |
462 +====++====+====+====+====+
463 | e0 || e0 | e1 | e2 | e3 |
464 +----++----+----+----+----+
465 | e1 || e1 | e0 | 0 | 0 |
466 +----++----+----+----+----+
467 | e2 || e2 | 0 | e0 | 0 |
468 +----++----+----+----+----+
469 | e3 || e3 | 0 | 0 | e0 |
470 +----++----+----+----+----+
473 M
= list(self
._multiplication
_table
) # copy
474 for i
in xrange(len(M
)):
475 # M had better be "square"
476 M
[i
] = [self
.monomial(i
)] + M
[i
]
477 M
= [["*"] + list(self
.gens())] + M
478 return table(M
, header_row
=True, header_column
=True, frame
=True)
481 def natural_basis(self
):
483 Return a more-natural representation of this algebra's basis.
485 Every finite-dimensional Euclidean Jordan Algebra is a direct
486 sum of five simple algebras, four of which comprise Hermitian
487 matrices. This method returns the original "natural" basis
488 for our underlying vector space. (Typically, the natural basis
489 is used to construct the multiplication table in the first place.)
491 Note that this will always return a matrix. The standard basis
492 in `R^n` will be returned as `n`-by-`1` column matrices.
496 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
497 ....: RealSymmetricEJA)
501 sage: J = RealSymmetricEJA(2)
503 Finite family {0: e0, 1: e1, 2: e2}
504 sage: J.natural_basis()
506 [1 0] [ 0 1/2*sqrt2] [0 0]
507 [0 0], [1/2*sqrt2 0], [0 1]
512 sage: J = JordanSpinEJA(2)
514 Finite family {0: e0, 1: e1}
515 sage: J.natural_basis()
522 if self
._natural
_basis
is None:
523 M
= self
.natural_basis_space()
524 return tuple( M(b
.to_vector()) for b
in self
.basis() )
526 return self
._natural
_basis
529 def natural_basis_space(self
):
531 Return the matrix space in which this algebra's natural basis
534 if self
._natural
_basis
is None or len(self
._natural
_basis
) == 0:
535 return MatrixSpace(self
.base_ring(), self
.dimension(), 1)
537 return self
._natural
_basis
[0].matrix_space()
541 def natural_inner_product(X
,Y
):
543 Compute the inner product of two naturally-represented elements.
545 For example in the real symmetric matrix EJA, this will compute
546 the trace inner-product of two n-by-n symmetric matrices. The
547 default should work for the real cartesian product EJA, the
548 Jordan spin EJA, and the real symmetric matrices. The others
549 will have to be overridden.
551 return (X
.conjugate_transpose()*Y
).trace()
557 Return the unit element of this algebra.
561 sage: from mjo.eja.eja_algebra import (RealCartesianProductEJA,
566 sage: J = RealCartesianProductEJA(5)
568 e0 + e1 + e2 + e3 + e4
572 The identity element acts like the identity::
574 sage: set_random_seed()
575 sage: J = random_eja()
576 sage: x = J.random_element()
577 sage: J.one()*x == x and x*J.one() == x
580 The matrix of the unit element's operator is the identity::
582 sage: set_random_seed()
583 sage: J = random_eja()
584 sage: actual = J.one().operator().matrix()
585 sage: expected = matrix.identity(J.base_ring(), J.dimension())
586 sage: actual == expected
590 # We can brute-force compute the matrices of the operators
591 # that correspond to the basis elements of this algebra.
592 # If some linear combination of those basis elements is the
593 # algebra identity, then the same linear combination of
594 # their matrices has to be the identity matrix.
596 # Of course, matrices aren't vectors in sage, so we have to
597 # appeal to the "long vectors" isometry.
598 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
600 # Now we use basis linear algebra to find the coefficients,
601 # of the matrices-as-vectors-linear-combination, which should
602 # work for the original algebra basis too.
603 A
= matrix
.column(self
.base_ring(), oper_vecs
)
605 # We used the isometry on the left-hand side already, but we
606 # still need to do it for the right-hand side. Recall that we
607 # wanted something that summed to the identity matrix.
608 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
610 # Now if there's an identity element in the algebra, this should work.
611 coeffs
= A
.solve_right(b
)
612 return self
.linear_combination(zip(self
.gens(), coeffs
))
615 def random_element(self
):
616 # Temporary workaround for https://trac.sagemath.org/ticket/28327
617 if self
.is_trivial():
620 s
= super(FiniteDimensionalEuclideanJordanAlgebra
, self
)
621 return s
.random_element()
623 def random_elements(self
, count
):
625 Return ``count`` random elements as a tuple.
629 sage: from mjo.eja.eja_algebra import JordanSpinEJA
633 sage: J = JordanSpinEJA(3)
634 sage: x,y,z = J.random_elements(3)
635 sage: all( [ x in J, y in J, z in J ])
637 sage: len( J.random_elements(10) ) == 10
641 return tuple( self
.random_element() for idx
in xrange(count
) )
646 Return the rank of this EJA.
650 The author knows of no algorithm to compute the rank of an EJA
651 where only the multiplication table is known. In lieu of one, we
652 require the rank to be specified when the algebra is created,
653 and simply pass along that number here.
657 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
658 ....: RealSymmetricEJA,
659 ....: ComplexHermitianEJA,
660 ....: QuaternionHermitianEJA,
665 The rank of the Jordan spin algebra is always two::
667 sage: JordanSpinEJA(2).rank()
669 sage: JordanSpinEJA(3).rank()
671 sage: JordanSpinEJA(4).rank()
674 The rank of the `n`-by-`n` Hermitian real, complex, or
675 quaternion matrices is `n`::
677 sage: RealSymmetricEJA(4).rank()
679 sage: ComplexHermitianEJA(3).rank()
681 sage: QuaternionHermitianEJA(2).rank()
686 Ensure that every EJA that we know how to construct has a
687 positive integer rank::
689 sage: set_random_seed()
690 sage: r = random_eja().rank()
691 sage: r in ZZ and r > 0
698 def vector_space(self
):
700 Return the vector space that underlies this algebra.
704 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
708 sage: J = RealSymmetricEJA(2)
709 sage: J.vector_space()
710 Vector space of dimension 3 over...
713 return self
.zero().to_vector().parent().ambient_vector_space()
716 Element
= FiniteDimensionalEuclideanJordanAlgebraElement
719 class KnownRankEJA(object):
721 A class for algebras that we actually know we can construct. The
722 main issue is that, for most of our methods to make sense, we need
723 to know the rank of our algebra. Thus we can't simply generate a
724 "random" algebra, or even check that a given basis and product
725 satisfy the axioms; because even if everything looks OK, we wouldn't
726 know the rank we need to actuallty build the thing.
728 Not really a subclass of FDEJA because doing that causes method
729 resolution errors, e.g.
731 TypeError: Error when calling the metaclass bases
732 Cannot create a consistent method resolution
733 order (MRO) for bases FiniteDimensionalEuclideanJordanAlgebra,
738 def _max_test_case_size():
740 Return an integer "size" that is an upper bound on the size of
741 this algebra when it is used in a random test
742 case. Unfortunately, the term "size" is quite vague -- when
743 dealing with `R^n` under either the Hadamard or Jordan spin
744 product, the "size" refers to the dimension `n`. When dealing
745 with a matrix algebra (real symmetric or complex/quaternion
746 Hermitian), it refers to the size of the matrix, which is
747 far less than the dimension of the underlying vector space.
749 We default to five in this class, which is safe in `R^n`. The
750 matrix algebra subclasses (or any class where the "size" is
751 interpreted to be far less than the dimension) should override
752 with a smaller number.
757 def random_instance(cls
, field
=QQ
, **kwargs
):
759 Return a random instance of this type of algebra.
761 Beware, this will crash for "most instances" because the
762 constructor below looks wrong.
764 n
= ZZ
.random_element(cls
._max
_test
_case
_size
()) + 1
765 return cls(n
, field
, **kwargs
)
768 class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra
,
771 Return the Euclidean Jordan Algebra corresponding to the set
772 `R^n` under the Hadamard product.
774 Note: this is nothing more than the Cartesian product of ``n``
775 copies of the spin algebra. Once Cartesian product algebras
776 are implemented, this can go.
780 sage: from mjo.eja.eja_algebra import RealCartesianProductEJA
784 This multiplication table can be verified by hand::
786 sage: J = RealCartesianProductEJA(3)
787 sage: e0,e1,e2 = J.gens()
803 We can change the generator prefix::
805 sage: RealCartesianProductEJA(3, prefix='r').gens()
809 def __init__(self
, n
, field
=QQ
, **kwargs
):
810 V
= VectorSpace(field
, n
)
811 mult_table
= [ [ V
.gen(i
)*(i
== j
) for j
in xrange(n
) ]
814 fdeja
= super(RealCartesianProductEJA
, self
)
815 return fdeja
.__init
__(field
, mult_table
, rank
=n
, **kwargs
)
817 def inner_product(self
, x
, y
):
819 Faster to reimplement than to use natural representations.
823 sage: from mjo.eja.eja_algebra import RealCartesianProductEJA
827 Ensure that this is the usual inner product for the algebras
830 sage: set_random_seed()
831 sage: J = RealCartesianProductEJA.random_instance()
832 sage: x,y = J.random_elements(2)
833 sage: X = x.natural_representation()
834 sage: Y = y.natural_representation()
835 sage: x.inner_product(y) == J.natural_inner_product(X,Y)
839 return x
.to_vector().inner_product(y
.to_vector())
844 Return a "random" finite-dimensional Euclidean Jordan Algebra.
848 For now, we choose a random natural number ``n`` (greater than zero)
849 and then give you back one of the following:
851 * The cartesian product of the rational numbers ``n`` times; this is
852 ``QQ^n`` with the Hadamard product.
854 * The Jordan spin algebra on ``QQ^n``.
856 * The ``n``-by-``n`` rational symmetric matrices with the symmetric
859 * The ``n``-by-``n`` complex-rational Hermitian matrices embedded
860 in the space of ``2n``-by-``2n`` real symmetric matrices.
862 * The ``n``-by-``n`` quaternion-rational Hermitian matrices embedded
863 in the space of ``4n``-by-``4n`` real symmetric matrices.
865 Later this might be extended to return Cartesian products of the
870 sage: from mjo.eja.eja_algebra import random_eja
875 Euclidean Jordan algebra of dimension...
878 classname
= choice(KnownRankEJA
.__subclasses
__())
879 return classname
.random_instance()
886 class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra
):
888 def _max_test_case_size():
889 # Play it safe, since this will be squared and the underlying
890 # field can have dimension 4 (quaternions) too.
893 def __init__(self
, field
, basis
, rank
, normalize_basis
=True, **kwargs
):
895 Compared to the superclass constructor, we take a basis instead of
896 a multiplication table because the latter can be computed in terms
897 of the former when the product is known (like it is here).
899 # Used in this class's fast _charpoly_coeff() override.
900 self
._basis
_normalizers
= None
902 # We're going to loop through this a few times, so now's a good
903 # time to ensure that it isn't a generator expression.
906 if rank
> 1 and normalize_basis
:
907 # We'll need sqrt(2) to normalize the basis, and this
908 # winds up in the multiplication table, so the whole
909 # algebra needs to be over the field extension.
910 R
= PolynomialRing(field
, 'z')
913 if p
.is_irreducible():
914 field
= NumberField(p
, 'sqrt2', embedding
=RLF(2).sqrt())
915 basis
= tuple( s
.change_ring(field
) for s
in basis
)
916 self
._basis
_normalizers
= tuple(
917 ~
(self
.natural_inner_product(s
,s
).sqrt()) for s
in basis
)
918 basis
= tuple(s
*c
for (s
,c
) in izip(basis
,self
._basis
_normalizers
))
920 Qs
= self
.multiplication_table_from_matrix_basis(basis
)
922 fdeja
= super(MatrixEuclideanJordanAlgebra
, self
)
923 return fdeja
.__init
__(field
,
931 def _charpoly_coeff(self
, i
):
933 Override the parent method with something that tries to compute
934 over a faster (non-extension) field.
936 if self
._basis
_normalizers
is None:
937 # We didn't normalize, so assume that the basis we started
938 # with had entries in a nice field.
939 return super(MatrixEuclideanJordanAlgebra
, self
)._charpoly
_coeff
(i
)
941 basis
= ( (b
/n
) for (b
,n
) in izip(self
.natural_basis(),
942 self
._basis
_normalizers
) )
943 field
= self
.base_ring().base_ring() # yeeeaahhhhhhh
944 J
= MatrixEuclideanJordanAlgebra(field
,
947 normalize_basis
=False)
948 (_
,x
,_
,_
) = J
._charpoly
_matrix
_system
()
949 p
= J
._charpoly
_coeff
(i
)
950 # p might be missing some vars, have to substitute "optionally"
951 pairs
= izip(x
.base_ring().gens(), self
._basis
_normalizers
)
952 substitutions
= { v: v*c for (v,c) in pairs }
953 return p
.subs(substitutions
)
957 def multiplication_table_from_matrix_basis(basis
):
959 At least three of the five simple Euclidean Jordan algebras have the
960 symmetric multiplication (A,B) |-> (AB + BA)/2, where the
961 multiplication on the right is matrix multiplication. Given a basis
962 for the underlying matrix space, this function returns a
963 multiplication table (obtained by looping through the basis
964 elements) for an algebra of those matrices.
966 # In S^2, for example, we nominally have four coordinates even
967 # though the space is of dimension three only. The vector space V
968 # is supposed to hold the entire long vector, and the subspace W
969 # of V will be spanned by the vectors that arise from symmetric
970 # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
971 field
= basis
[0].base_ring()
972 dimension
= basis
[0].nrows()
974 V
= VectorSpace(field
, dimension
**2)
975 W
= V
.span_of_basis( _mat2vec(s
) for s
in basis
)
977 mult_table
= [[W
.zero() for j
in xrange(n
)] for i
in xrange(n
)]
980 mat_entry
= (basis
[i
]*basis
[j
] + basis
[j
]*basis
[i
])/2
981 mult_table
[i
][j
] = W
.coordinate_vector(_mat2vec(mat_entry
))
989 Embed the matrix ``M`` into a space of real matrices.
991 The matrix ``M`` can have entries in any field at the moment:
992 the real numbers, complex numbers, or quaternions. And although
993 they are not a field, we can probably support octonions at some
994 point, too. This function returns a real matrix that "acts like"
995 the original with respect to matrix multiplication; i.e.
997 real_embed(M*N) = real_embed(M)*real_embed(N)
1000 raise NotImplementedError
1004 def real_unembed(M
):
1006 The inverse of :meth:`real_embed`.
1008 raise NotImplementedError
1012 def natural_inner_product(cls
,X
,Y
):
1013 Xu
= cls
.real_unembed(X
)
1014 Yu
= cls
.real_unembed(Y
)
1015 tr
= (Xu
*Yu
).trace()
1017 # It's real already.
1020 # Otherwise, try the thing that works for complex numbers; and
1021 # if that doesn't work, the thing that works for quaternions.
1023 return tr
.vector()[0] # real part, imag part is index 1
1024 except AttributeError:
1025 # A quaternions doesn't have a vector() method, but does
1026 # have coefficient_tuple() method that returns the
1027 # coefficients of 1, i, j, and k -- in that order.
1028 return tr
.coefficient_tuple()[0]
1031 class RealMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1035 The identity function, for embedding real matrices into real
1041 def real_unembed(M
):
1043 The identity function, for unembedding real matrices from real
1049 class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra
, KnownRankEJA
):
1051 The rank-n simple EJA consisting of real symmetric n-by-n
1052 matrices, the usual symmetric Jordan product, and the trace inner
1053 product. It has dimension `(n^2 + n)/2` over the reals.
1057 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1061 sage: J = RealSymmetricEJA(2)
1062 sage: e0, e1, e2 = J.gens()
1072 The dimension of this algebra is `(n^2 + n) / 2`::
1074 sage: set_random_seed()
1075 sage: n_max = RealSymmetricEJA._max_test_case_size()
1076 sage: n = ZZ.random_element(1, n_max)
1077 sage: J = RealSymmetricEJA(n)
1078 sage: J.dimension() == (n^2 + n)/2
1081 The Jordan multiplication is what we think it is::
1083 sage: set_random_seed()
1084 sage: J = RealSymmetricEJA.random_instance()
1085 sage: x,y = J.random_elements(2)
1086 sage: actual = (x*y).natural_representation()
1087 sage: X = x.natural_representation()
1088 sage: Y = y.natural_representation()
1089 sage: expected = (X*Y + Y*X)/2
1090 sage: actual == expected
1092 sage: J(expected) == x*y
1095 We can change the generator prefix::
1097 sage: RealSymmetricEJA(3, prefix='q').gens()
1098 (q0, q1, q2, q3, q4, q5)
1100 Our natural basis is normalized with respect to the natural inner
1101 product unless we specify otherwise::
1103 sage: set_random_seed()
1104 sage: J = RealSymmetricEJA.random_instance()
1105 sage: all( b.norm() == 1 for b in J.gens() )
1108 Since our natural basis is normalized with respect to the natural
1109 inner product, and since we know that this algebra is an EJA, any
1110 left-multiplication operator's matrix will be symmetric because
1111 natural->EJA basis representation is an isometry and within the EJA
1112 the operator is self-adjoint by the Jordan axiom::
1114 sage: set_random_seed()
1115 sage: x = RealSymmetricEJA.random_instance().random_element()
1116 sage: x.operator().matrix().is_symmetric()
1121 def _denormalized_basis(cls
, n
, field
):
1123 Return a basis for the space of real symmetric n-by-n matrices.
1127 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1131 sage: set_random_seed()
1132 sage: n = ZZ.random_element(1,5)
1133 sage: B = RealSymmetricEJA._denormalized_basis(n,QQ)
1134 sage: all( M.is_symmetric() for M in B)
1138 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1142 for j
in xrange(i
+1):
1143 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
1147 Sij
= Eij
+ Eij
.transpose()
1153 def _max_test_case_size():
1154 return 4 # Dimension 10
1157 def __init__(self
, n
, field
=QQ
, **kwargs
):
1158 basis
= self
._denormalized
_basis
(n
, field
)
1159 super(RealSymmetricEJA
, self
).__init
__(field
, basis
, n
, **kwargs
)
1162 class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1166 Embed the n-by-n complex matrix ``M`` into the space of real
1167 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1168 bi` to the block matrix ``[[a,b],[-b,a]]``.
1172 sage: from mjo.eja.eja_algebra import \
1173 ....: ComplexMatrixEuclideanJordanAlgebra
1177 sage: F = QuadraticField(-1, 'i')
1178 sage: x1 = F(4 - 2*i)
1179 sage: x2 = F(1 + 2*i)
1182 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1183 sage: ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1192 Embedding is a homomorphism (isomorphism, in fact)::
1194 sage: set_random_seed()
1195 sage: n_max = ComplexMatrixEuclideanJordanAlgebra._max_test_case_size()
1196 sage: n = ZZ.random_element(n_max)
1197 sage: F = QuadraticField(-1, 'i')
1198 sage: X = random_matrix(F, n)
1199 sage: Y = random_matrix(F, n)
1200 sage: Xe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X)
1201 sage: Ye = ComplexMatrixEuclideanJordanAlgebra.real_embed(Y)
1202 sage: XYe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1209 raise ValueError("the matrix 'M' must be square")
1210 field
= M
.base_ring()
1213 a
= z
.vector()[0] # real part, I guess
1214 b
= z
.vector()[1] # imag part, I guess
1215 blocks
.append(matrix(field
, 2, [[a
,b
],[-b
,a
]]))
1217 # We can drop the imaginaries here.
1218 return matrix
.block(field
.base_ring(), n
, blocks
)
1222 def real_unembed(M
):
1224 The inverse of _embed_complex_matrix().
1228 sage: from mjo.eja.eja_algebra import \
1229 ....: ComplexMatrixEuclideanJordanAlgebra
1233 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1234 ....: [-2, 1, -4, 3],
1235 ....: [ 9, 10, 11, 12],
1236 ....: [-10, 9, -12, 11] ])
1237 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(A)
1239 [ 10*i + 9 12*i + 11]
1243 Unembedding is the inverse of embedding::
1245 sage: set_random_seed()
1246 sage: F = QuadraticField(-1, 'i')
1247 sage: M = random_matrix(F, 3)
1248 sage: Me = ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1249 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1255 raise ValueError("the matrix 'M' must be square")
1256 if not n
.mod(2).is_zero():
1257 raise ValueError("the matrix 'M' must be a complex embedding")
1259 field
= M
.base_ring() # This should already have sqrt2
1260 R
= PolynomialRing(field
, 'z')
1262 F
= NumberField(z
**2 + 1,'i', embedding
=CLF(-1).sqrt())
1265 # Go top-left to bottom-right (reading order), converting every
1266 # 2-by-2 block we see to a single complex element.
1268 for k
in xrange(n
/2):
1269 for j
in xrange(n
/2):
1270 submat
= M
[2*k
:2*k
+2,2*j
:2*j
+2]
1271 if submat
[0,0] != submat
[1,1]:
1272 raise ValueError('bad on-diagonal submatrix')
1273 if submat
[0,1] != -submat
[1,0]:
1274 raise ValueError('bad off-diagonal submatrix')
1275 z
= submat
[0,0] + submat
[0,1]*i
1278 return matrix(F
, n
/2, elements
)
1282 def natural_inner_product(cls
,X
,Y
):
1284 Compute a natural inner product in this algebra directly from
1289 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1293 This gives the same answer as the slow, default method implemented
1294 in :class:`MatrixEuclideanJordanAlgebra`::
1296 sage: set_random_seed()
1297 sage: J = ComplexHermitianEJA.random_instance()
1298 sage: x,y = J.random_elements(2)
1299 sage: Xe = x.natural_representation()
1300 sage: Ye = y.natural_representation()
1301 sage: X = ComplexHermitianEJA.real_unembed(Xe)
1302 sage: Y = ComplexHermitianEJA.real_unembed(Ye)
1303 sage: expected = (X*Y).trace().vector()[0]
1304 sage: actual = ComplexHermitianEJA.natural_inner_product(Xe,Ye)
1305 sage: actual == expected
1309 return RealMatrixEuclideanJordanAlgebra
.natural_inner_product(X
,Y
)/2
1312 class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra
, KnownRankEJA
):
1314 The rank-n simple EJA consisting of complex Hermitian n-by-n
1315 matrices over the real numbers, the usual symmetric Jordan product,
1316 and the real-part-of-trace inner product. It has dimension `n^2` over
1321 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1325 The dimension of this algebra is `n^2`::
1327 sage: set_random_seed()
1328 sage: n_max = ComplexHermitianEJA._max_test_case_size()
1329 sage: n = ZZ.random_element(1, n_max)
1330 sage: J = ComplexHermitianEJA(n)
1331 sage: J.dimension() == n^2
1334 The Jordan multiplication is what we think it is::
1336 sage: set_random_seed()
1337 sage: J = ComplexHermitianEJA.random_instance()
1338 sage: x,y = J.random_elements(2)
1339 sage: actual = (x*y).natural_representation()
1340 sage: X = x.natural_representation()
1341 sage: Y = y.natural_representation()
1342 sage: expected = (X*Y + Y*X)/2
1343 sage: actual == expected
1345 sage: J(expected) == x*y
1348 We can change the generator prefix::
1350 sage: ComplexHermitianEJA(2, prefix='z').gens()
1353 Our natural basis is normalized with respect to the natural inner
1354 product unless we specify otherwise::
1356 sage: set_random_seed()
1357 sage: J = ComplexHermitianEJA.random_instance()
1358 sage: all( b.norm() == 1 for b in J.gens() )
1361 Since our natural basis is normalized with respect to the natural
1362 inner product, and since we know that this algebra is an EJA, any
1363 left-multiplication operator's matrix will be symmetric because
1364 natural->EJA basis representation is an isometry and within the EJA
1365 the operator is self-adjoint by the Jordan axiom::
1367 sage: set_random_seed()
1368 sage: x = ComplexHermitianEJA.random_instance().random_element()
1369 sage: x.operator().matrix().is_symmetric()
1375 def _denormalized_basis(cls
, n
, field
):
1377 Returns a basis for the space of complex Hermitian n-by-n matrices.
1379 Why do we embed these? Basically, because all of numerical linear
1380 algebra assumes that you're working with vectors consisting of `n`
1381 entries from a field and scalars from the same field. There's no way
1382 to tell SageMath that (for example) the vectors contain complex
1383 numbers, while the scalar field is real.
1387 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1391 sage: set_random_seed()
1392 sage: n = ZZ.random_element(1,5)
1393 sage: field = QuadraticField(2, 'sqrt2')
1394 sage: B = ComplexHermitianEJA._denormalized_basis(n, field)
1395 sage: all( M.is_symmetric() for M in B)
1399 R
= PolynomialRing(field
, 'z')
1401 F
= NumberField(z
**2 + 1, 'I', embedding
=CLF(-1).sqrt())
1404 # This is like the symmetric case, but we need to be careful:
1406 # * We want conjugate-symmetry, not just symmetry.
1407 # * The diagonal will (as a result) be real.
1411 for j
in xrange(i
+1):
1412 Eij
= matrix(F
, n
, lambda k
,l
: k
==i
and l
==j
)
1414 Sij
= cls
.real_embed(Eij
)
1417 # The second one has a minus because it's conjugated.
1418 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
1420 Sij_imag
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
1423 # Since we embedded these, we can drop back to the "field" that we
1424 # started with instead of the complex extension "F".
1425 return ( s
.change_ring(field
) for s
in S
)
1428 def __init__(self
, n
, field
=QQ
, **kwargs
):
1429 basis
= self
._denormalized
_basis
(n
,field
)
1430 super(ComplexHermitianEJA
,self
).__init
__(field
, basis
, n
, **kwargs
)
1433 class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1437 Embed the n-by-n quaternion matrix ``M`` into the space of real
1438 matrices of size 4n-by-4n by first sending each quaternion entry `z
1439 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
1440 c+di],[-c + di, a-bi]]`, and then embedding those into a real
1445 sage: from mjo.eja.eja_algebra import \
1446 ....: QuaternionMatrixEuclideanJordanAlgebra
1450 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1451 sage: i,j,k = Q.gens()
1452 sage: x = 1 + 2*i + 3*j + 4*k
1453 sage: M = matrix(Q, 1, [[x]])
1454 sage: QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1460 Embedding is a homomorphism (isomorphism, in fact)::
1462 sage: set_random_seed()
1463 sage: n_max = QuaternionMatrixEuclideanJordanAlgebra._max_test_case_size()
1464 sage: n = ZZ.random_element(n_max)
1465 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1466 sage: X = random_matrix(Q, n)
1467 sage: Y = random_matrix(Q, n)
1468 sage: Xe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X)
1469 sage: Ye = QuaternionMatrixEuclideanJordanAlgebra.real_embed(Y)
1470 sage: XYe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1475 quaternions
= M
.base_ring()
1478 raise ValueError("the matrix 'M' must be square")
1480 F
= QuadraticField(-1, 'i')
1485 t
= z
.coefficient_tuple()
1490 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
1491 [-c
+ d
*i
, a
- b
*i
]])
1492 realM
= ComplexMatrixEuclideanJordanAlgebra
.real_embed(cplxM
)
1493 blocks
.append(realM
)
1495 # We should have real entries by now, so use the realest field
1496 # we've got for the return value.
1497 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
1502 def real_unembed(M
):
1504 The inverse of _embed_quaternion_matrix().
1508 sage: from mjo.eja.eja_algebra import \
1509 ....: QuaternionMatrixEuclideanJordanAlgebra
1513 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
1514 ....: [-2, 1, -4, 3],
1515 ....: [-3, 4, 1, -2],
1516 ....: [-4, -3, 2, 1]])
1517 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(M)
1518 [1 + 2*i + 3*j + 4*k]
1522 Unembedding is the inverse of embedding::
1524 sage: set_random_seed()
1525 sage: Q = QuaternionAlgebra(QQ, -1, -1)
1526 sage: M = random_matrix(Q, 3)
1527 sage: Me = QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1528 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1534 raise ValueError("the matrix 'M' must be square")
1535 if not n
.mod(4).is_zero():
1536 raise ValueError("the matrix 'M' must be a complex embedding")
1538 # Use the base ring of the matrix to ensure that its entries can be
1539 # multiplied by elements of the quaternion algebra.
1540 field
= M
.base_ring()
1541 Q
= QuaternionAlgebra(field
,-1,-1)
1544 # Go top-left to bottom-right (reading order), converting every
1545 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
1548 for l
in xrange(n
/4):
1549 for m
in xrange(n
/4):
1550 submat
= ComplexMatrixEuclideanJordanAlgebra
.real_unembed(
1551 M
[4*l
:4*l
+4,4*m
:4*m
+4] )
1552 if submat
[0,0] != submat
[1,1].conjugate():
1553 raise ValueError('bad on-diagonal submatrix')
1554 if submat
[0,1] != -submat
[1,0].conjugate():
1555 raise ValueError('bad off-diagonal submatrix')
1556 z
= submat
[0,0].vector()[0] # real part
1557 z
+= submat
[0,0].vector()[1]*i
# imag part
1558 z
+= submat
[0,1].vector()[0]*j
# real part
1559 z
+= submat
[0,1].vector()[1]*k
# imag part
1562 return matrix(Q
, n
/4, elements
)
1566 def natural_inner_product(cls
,X
,Y
):
1568 Compute a natural inner product in this algebra directly from
1573 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1577 This gives the same answer as the slow, default method implemented
1578 in :class:`MatrixEuclideanJordanAlgebra`::
1580 sage: set_random_seed()
1581 sage: J = QuaternionHermitianEJA.random_instance()
1582 sage: x,y = J.random_elements(2)
1583 sage: Xe = x.natural_representation()
1584 sage: Ye = y.natural_representation()
1585 sage: X = QuaternionHermitianEJA.real_unembed(Xe)
1586 sage: Y = QuaternionHermitianEJA.real_unembed(Ye)
1587 sage: expected = (X*Y).trace().coefficient_tuple()[0]
1588 sage: actual = QuaternionHermitianEJA.natural_inner_product(Xe,Ye)
1589 sage: actual == expected
1593 return RealMatrixEuclideanJordanAlgebra
.natural_inner_product(X
,Y
)/4
1596 class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra
,
1599 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
1600 matrices, the usual symmetric Jordan product, and the
1601 real-part-of-trace inner product. It has dimension `2n^2 - n` over
1606 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1610 The dimension of this algebra is `2*n^2 - n`::
1612 sage: set_random_seed()
1613 sage: n_max = QuaternionHermitianEJA._max_test_case_size()
1614 sage: n = ZZ.random_element(1, n_max)
1615 sage: J = QuaternionHermitianEJA(n)
1616 sage: J.dimension() == 2*(n^2) - n
1619 The Jordan multiplication is what we think it is::
1621 sage: set_random_seed()
1622 sage: J = QuaternionHermitianEJA.random_instance()
1623 sage: x,y = J.random_elements(2)
1624 sage: actual = (x*y).natural_representation()
1625 sage: X = x.natural_representation()
1626 sage: Y = y.natural_representation()
1627 sage: expected = (X*Y + Y*X)/2
1628 sage: actual == expected
1630 sage: J(expected) == x*y
1633 We can change the generator prefix::
1635 sage: QuaternionHermitianEJA(2, prefix='a').gens()
1636 (a0, a1, a2, a3, a4, a5)
1638 Our natural basis is normalized with respect to the natural inner
1639 product unless we specify otherwise::
1641 sage: set_random_seed()
1642 sage: J = QuaternionHermitianEJA.random_instance()
1643 sage: all( b.norm() == 1 for b in J.gens() )
1646 Since our natural basis is normalized with respect to the natural
1647 inner product, and since we know that this algebra is an EJA, any
1648 left-multiplication operator's matrix will be symmetric because
1649 natural->EJA basis representation is an isometry and within the EJA
1650 the operator is self-adjoint by the Jordan axiom::
1652 sage: set_random_seed()
1653 sage: x = QuaternionHermitianEJA.random_instance().random_element()
1654 sage: x.operator().matrix().is_symmetric()
1659 def _denormalized_basis(cls
, n
, field
):
1661 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
1663 Why do we embed these? Basically, because all of numerical
1664 linear algebra assumes that you're working with vectors consisting
1665 of `n` entries from a field and scalars from the same field. There's
1666 no way to tell SageMath that (for example) the vectors contain
1667 complex numbers, while the scalar field is real.
1671 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1675 sage: set_random_seed()
1676 sage: n = ZZ.random_element(1,5)
1677 sage: B = QuaternionHermitianEJA._denormalized_basis(n,QQ)
1678 sage: all( M.is_symmetric() for M in B )
1682 Q
= QuaternionAlgebra(QQ
,-1,-1)
1685 # This is like the symmetric case, but we need to be careful:
1687 # * We want conjugate-symmetry, not just symmetry.
1688 # * The diagonal will (as a result) be real.
1692 for j
in xrange(i
+1):
1693 Eij
= matrix(Q
, n
, lambda k
,l
: k
==i
and l
==j
)
1695 Sij
= cls
.real_embed(Eij
)
1698 # The second, third, and fourth ones have a minus
1699 # because they're conjugated.
1700 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
1702 Sij_I
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
1704 Sij_J
= cls
.real_embed(J
*Eij
- J
*Eij
.transpose())
1706 Sij_K
= cls
.real_embed(K
*Eij
- K
*Eij
.transpose())
1711 def __init__(self
, n
, field
=QQ
, **kwargs
):
1712 basis
= self
._denormalized
_basis
(n
,field
)
1713 super(QuaternionHermitianEJA
,self
).__init
__(field
, basis
, n
, **kwargs
)
1716 class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra
, KnownRankEJA
):
1718 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
1719 with the usual inner product and jordan product ``x*y =
1720 (<x_bar,y_bar>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
1725 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1729 This multiplication table can be verified by hand::
1731 sage: J = JordanSpinEJA(4)
1732 sage: e0,e1,e2,e3 = J.gens()
1748 We can change the generator prefix::
1750 sage: JordanSpinEJA(2, prefix='B').gens()
1754 def __init__(self
, n
, field
=QQ
, **kwargs
):
1755 V
= VectorSpace(field
, n
)
1756 mult_table
= [[V
.zero() for j
in xrange(n
)] for i
in xrange(n
)]
1766 z0
= x
.inner_product(y
)
1767 zbar
= y0
*xbar
+ x0
*ybar
1768 z
= V([z0
] + zbar
.list())
1769 mult_table
[i
][j
] = z
1771 # The rank of the spin algebra is two, unless we're in a
1772 # one-dimensional ambient space (because the rank is bounded by
1773 # the ambient dimension).
1774 fdeja
= super(JordanSpinEJA
, self
)
1775 return fdeja
.__init
__(field
, mult_table
, rank
=min(n
,2), **kwargs
)
1777 def inner_product(self
, x
, y
):
1779 Faster to reimplement than to use natural representations.
1783 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1787 Ensure that this is the usual inner product for the algebras
1790 sage: set_random_seed()
1791 sage: J = JordanSpinEJA.random_instance()
1792 sage: x,y = J.random_elements(2)
1793 sage: X = x.natural_representation()
1794 sage: Y = y.natural_representation()
1795 sage: x.inner_product(y) == J.natural_inner_product(X,Y)
1799 return x
.to_vector().inner_product(y
.to_vector())