2 Euclidean Jordan Algebras. These are formally-real Jordan Algebras;
3 specifically those where u^2 + v^2 = 0 implies that u = v = 0. They
4 are used in optimization, and have some additional nice methods beyond
5 what can be supported in a general Jordan Algebra.
8 from itertools
import izip
, repeat
10 from sage
.algebras
.quatalg
.quaternion_algebra
import QuaternionAlgebra
11 from sage
.categories
.magmatic_algebras
import MagmaticAlgebras
12 from sage
.combinat
.free_module
import CombinatorialFreeModule
13 from sage
.matrix
.constructor
import matrix
14 from sage
.matrix
.matrix_space
import MatrixSpace
15 from sage
.misc
.cachefunc
import cached_method
16 from sage
.misc
.prandom
import choice
17 from sage
.misc
.table
import table
18 from sage
.modules
.free_module
import FreeModule
, VectorSpace
19 from sage
.rings
.all
import (ZZ
, QQ
, RR
, RLF
, CLF
,
22 from mjo
.eja
.eja_element
import FiniteDimensionalEuclideanJordanAlgebraElement
23 from mjo
.eja
.eja_utils
import _mat2vec
25 class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule
):
26 # This is an ugly hack needed to prevent the category framework
27 # from implementing a coercion from our base ring (e.g. the
28 # rationals) into the algebra. First of all -- such a coercion is
29 # nonsense to begin with. But more importantly, it tries to do so
30 # in the category of rings, and since our algebras aren't
31 # associative they generally won't be rings.
32 _no_generic_basering_coercion
= True
45 sage: from mjo.eja.eja_algebra import (JordanSpinEJA, random_eja)
49 By definition, Jordan multiplication commutes::
51 sage: set_random_seed()
52 sage: J = random_eja()
53 sage: x,y = J.random_elements(2)
59 The ``field`` we're given must be real::
61 sage: JordanSpinEJA(2,QQbar)
62 Traceback (most recent call last):
64 ValueError: field is not real
68 if not field
.is_subring(RR
):
69 # Note: this does return true for the real algebraic
70 # field, and any quadratic field where we've specified
72 raise ValueError('field is not real')
75 self
._natural
_basis
= natural_basis
78 category
= MagmaticAlgebras(field
).FiniteDimensional()
79 category
= category
.WithBasis().Unital()
81 fda
= super(FiniteDimensionalEuclideanJordanAlgebra
, self
)
83 range(len(mult_table
)),
86 self
.print_options(bracket
='')
88 # The multiplication table we're given is necessarily in terms
89 # of vectors, because we don't have an algebra yet for
90 # anything to be an element of. However, it's faster in the
91 # long run to have the multiplication table be in terms of
92 # algebra elements. We do this after calling the superclass
93 # constructor so that from_vector() knows what to do.
94 self
._multiplication
_table
= [ map(lambda x
: self
.from_vector(x
), ls
)
95 for ls
in mult_table
]
98 def _element_constructor_(self
, elt
):
100 Construct an element of this algebra from its natural
103 This gets called only after the parent element _call_ method
104 fails to find a coercion for the argument.
108 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
109 ....: RealCartesianProductEJA,
110 ....: RealSymmetricEJA)
114 The identity in `S^n` is converted to the identity in the EJA::
116 sage: J = RealSymmetricEJA(3)
117 sage: I = matrix.identity(QQ,3)
118 sage: J(I) == J.one()
121 This skew-symmetric matrix can't be represented in the EJA::
123 sage: J = RealSymmetricEJA(3)
124 sage: A = matrix(QQ,3, lambda i,j: i-j)
126 Traceback (most recent call last):
128 ArithmeticError: vector is not in free module
132 Ensure that we can convert any element of the two non-matrix
133 simple algebras (whose natural representations are their usual
134 vector representations) back and forth faithfully::
136 sage: set_random_seed()
137 sage: J = RealCartesianProductEJA.random_instance()
138 sage: x = J.random_element()
139 sage: J(x.to_vector().column()) == x
141 sage: J = JordanSpinEJA.random_instance()
142 sage: x = J.random_element()
143 sage: J(x.to_vector().column()) == x
148 # The superclass implementation of random_element()
149 # needs to be able to coerce "0" into the algebra.
152 natural_basis
= self
.natural_basis()
153 basis_space
= natural_basis
[0].matrix_space()
154 if elt
not in basis_space
:
155 raise ValueError("not a naturally-represented algebra element")
157 # Thanks for nothing! Matrix spaces aren't vector spaces in
158 # Sage, so we have to figure out its natural-basis coordinates
159 # ourselves. We use the basis space's ring instead of the
160 # element's ring because the basis space might be an algebraic
161 # closure whereas the base ring of the 3-by-3 identity matrix
162 # could be QQ instead of QQbar.
163 V
= VectorSpace(basis_space
.base_ring(), elt
.nrows()*elt
.ncols())
164 W
= V
.span_of_basis( _mat2vec(s
) for s
in natural_basis
)
165 coords
= W
.coordinate_vector(_mat2vec(elt
))
166 return self
.from_vector(coords
)
171 Return a string representation of ``self``.
175 sage: from mjo.eja.eja_algebra import JordanSpinEJA
179 Ensure that it says what we think it says::
181 sage: JordanSpinEJA(2, field=QQ)
182 Euclidean Jordan algebra of dimension 2 over Rational Field
183 sage: JordanSpinEJA(3, field=RDF)
184 Euclidean Jordan algebra of dimension 3 over Real Double Field
187 fmt
= "Euclidean Jordan algebra of dimension {} over {}"
188 return fmt
.format(self
.dimension(), self
.base_ring())
190 def product_on_basis(self
, i
, j
):
191 return self
._multiplication
_table
[i
][j
]
193 def _a_regular_element(self
):
195 Guess a regular element. Needed to compute the basis for our
196 characteristic polynomial coefficients.
200 sage: from mjo.eja.eja_algebra import random_eja
204 Ensure that this hacky method succeeds for every algebra that we
205 know how to construct::
207 sage: set_random_seed()
208 sage: J = random_eja()
209 sage: J._a_regular_element().is_regular()
214 z
= self
.sum( (i
+1)*gs
[i
] for i
in range(len(gs
)) )
215 if not z
.is_regular():
216 raise ValueError("don't know a regular element")
221 def _charpoly_basis_space(self
):
223 Return the vector space spanned by the basis used in our
224 characteristic polynomial coefficients. This is used not only to
225 compute those coefficients, but also any time we need to
226 evaluate the coefficients (like when we compute the trace or
229 z
= self
._a
_regular
_element
()
230 # Don't use the parent vector space directly here in case this
231 # happens to be a subalgebra. In that case, we would be e.g.
232 # two-dimensional but span_of_basis() would expect three
234 V
= VectorSpace(self
.base_ring(), self
.vector_space().dimension())
235 basis
= [ (z
**k
).to_vector() for k
in range(self
.rank()) ]
236 V1
= V
.span_of_basis( basis
)
237 b
= (V1
.basis() + V1
.complement().basis())
238 return V
.span_of_basis(b
)
243 def _charpoly_coeff(self
, i
):
245 Return the coefficient polynomial "a_{i}" of this algebra's
246 general characteristic polynomial.
248 Having this be a separate cached method lets us compute and
249 store the trace/determinant (a_{r-1} and a_{0} respectively)
250 separate from the entire characteristic polynomial.
252 (A_of_x
, x
, xr
, detA
) = self
._charpoly
_matrix
_system
()
253 R
= A_of_x
.base_ring()
258 # Guaranteed by theory
261 # Danger: the in-place modification is done for performance
262 # reasons (reconstructing a matrix with huge polynomial
263 # entries is slow), but I don't know how cached_method works,
264 # so it's highly possible that we're modifying some global
265 # list variable by reference, here. In other words, you
266 # probably shouldn't call this method twice on the same
267 # algebra, at the same time, in two threads
268 Ai_orig
= A_of_x
.column(i
)
269 A_of_x
.set_column(i
,xr
)
270 numerator
= A_of_x
.det()
271 A_of_x
.set_column(i
,Ai_orig
)
273 # We're relying on the theory here to ensure that each a_i is
274 # indeed back in R, and the added negative signs are to make
275 # the whole charpoly expression sum to zero.
276 return R(-numerator
/detA
)
280 def _charpoly_matrix_system(self
):
282 Compute the matrix whose entries A_ij are polynomials in
283 X1,...,XN, the vector ``x`` of variables X1,...,XN, the vector
284 corresponding to `x^r` and the determinent of the matrix A =
285 [A_ij]. In other words, all of the fixed (cachable) data needed
286 to compute the coefficients of the characteristic polynomial.
291 # Turn my vector space into a module so that "vectors" can
292 # have multivatiate polynomial entries.
293 names
= tuple('X' + str(i
) for i
in range(1,n
+1))
294 R
= PolynomialRing(self
.base_ring(), names
)
296 # Using change_ring() on the parent's vector space doesn't work
297 # here because, in a subalgebra, that vector space has a basis
298 # and change_ring() tries to bring the basis along with it. And
299 # that doesn't work unless the new ring is a PID, which it usually
303 # Now let x = (X1,X2,...,Xn) be the vector whose entries are
307 # And figure out the "left multiplication by x" matrix in
310 monomial_matrices
= [ self
.monomial(i
).operator().matrix()
311 for i
in range(n
) ] # don't recompute these!
313 ek
= self
.monomial(k
).to_vector()
315 sum( x
[i
]*(monomial_matrices
[i
]*ek
)
316 for i
in range(n
) ) )
317 Lx
= matrix
.column(R
, lmbx_cols
)
319 # Now we can compute powers of x "symbolically"
320 x_powers
= [self
.one().to_vector(), x
]
321 for d
in range(2, r
+1):
322 x_powers
.append( Lx
*(x_powers
[-1]) )
324 idmat
= matrix
.identity(R
, n
)
326 W
= self
._charpoly
_basis
_space
()
327 W
= W
.change_ring(R
.fraction_field())
329 # Starting with the standard coordinates x = (X1,X2,...,Xn)
330 # and then converting the entries to W-coordinates allows us
331 # to pass in the standard coordinates to the charpoly and get
332 # back the right answer. Specifically, with x = (X1,X2,...,Xn),
335 # W.coordinates(x^2) eval'd at (standard z-coords)
339 # W-coords of (standard coords of x^2 eval'd at std-coords of z)
341 # We want the middle equivalent thing in our matrix, but use
342 # the first equivalent thing instead so that we can pass in
343 # standard coordinates.
344 x_powers
= [ W
.coordinate_vector(xp
) for xp
in x_powers
]
345 l2
= [idmat
.column(k
-1) for k
in range(r
+1, n
+1)]
346 A_of_x
= matrix
.column(R
, n
, (x_powers
[:r
] + l2
))
347 return (A_of_x
, x
, x_powers
[r
], A_of_x
.det())
351 def characteristic_polynomial(self
):
353 Return a characteristic polynomial that works for all elements
356 The resulting polynomial has `n+1` variables, where `n` is the
357 dimension of this algebra. The first `n` variables correspond to
358 the coordinates of an algebra element: when evaluated at the
359 coordinates of an algebra element with respect to a certain
360 basis, the result is a univariate polynomial (in the one
361 remaining variable ``t``), namely the characteristic polynomial
366 sage: from mjo.eja.eja_algebra import JordanSpinEJA
370 The characteristic polynomial in the spin algebra is given in
371 Alizadeh, Example 11.11::
373 sage: J = JordanSpinEJA(3)
374 sage: p = J.characteristic_polynomial(); p
375 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
376 sage: xvec = J.one().to_vector()
384 # The list of coefficient polynomials a_0, a_1, a_2, ..., a_n.
385 a
= [ self
._charpoly
_coeff
(i
) for i
in range(r
+1) ]
387 # We go to a bit of trouble here to reorder the
388 # indeterminates, so that it's easier to evaluate the
389 # characteristic polynomial at x's coordinates and get back
390 # something in terms of t, which is what we want.
392 S
= PolynomialRing(self
.base_ring(),'t')
394 S
= PolynomialRing(S
, R
.variable_names())
397 return sum( a
[k
]*(t
**k
) for k
in xrange(len(a
)) )
400 def inner_product(self
, x
, y
):
402 The inner product associated with this Euclidean Jordan algebra.
404 Defaults to the trace inner product, but can be overridden by
405 subclasses if they are sure that the necessary properties are
410 sage: from mjo.eja.eja_algebra import random_eja
414 Our inner product is "associative," which means the following for
415 a symmetric bilinear form::
417 sage: set_random_seed()
418 sage: J = random_eja()
419 sage: x,y,z = J.random_elements(3)
420 sage: (x*y).inner_product(z) == y.inner_product(x*z)
424 X
= x
.natural_representation()
425 Y
= y
.natural_representation()
426 return self
.natural_inner_product(X
,Y
)
429 def is_trivial(self
):
431 Return whether or not this algebra is trivial.
433 A trivial algebra contains only the zero element.
437 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
441 sage: J = ComplexHermitianEJA(3)
446 return self
.dimension() == 0
449 def multiplication_table(self
):
451 Return a visual representation of this algebra's multiplication
452 table (on basis elements).
456 sage: from mjo.eja.eja_algebra import JordanSpinEJA
460 sage: J = JordanSpinEJA(4)
461 sage: J.multiplication_table()
462 +----++----+----+----+----+
463 | * || e0 | e1 | e2 | e3 |
464 +====++====+====+====+====+
465 | e0 || e0 | e1 | e2 | e3 |
466 +----++----+----+----+----+
467 | e1 || e1 | e0 | 0 | 0 |
468 +----++----+----+----+----+
469 | e2 || e2 | 0 | e0 | 0 |
470 +----++----+----+----+----+
471 | e3 || e3 | 0 | 0 | e0 |
472 +----++----+----+----+----+
475 M
= list(self
._multiplication
_table
) # copy
476 for i
in xrange(len(M
)):
477 # M had better be "square"
478 M
[i
] = [self
.monomial(i
)] + M
[i
]
479 M
= [["*"] + list(self
.gens())] + M
480 return table(M
, header_row
=True, header_column
=True, frame
=True)
483 def natural_basis(self
):
485 Return a more-natural representation of this algebra's basis.
487 Every finite-dimensional Euclidean Jordan Algebra is a direct
488 sum of five simple algebras, four of which comprise Hermitian
489 matrices. This method returns the original "natural" basis
490 for our underlying vector space. (Typically, the natural basis
491 is used to construct the multiplication table in the first place.)
493 Note that this will always return a matrix. The standard basis
494 in `R^n` will be returned as `n`-by-`1` column matrices.
498 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
499 ....: RealSymmetricEJA)
503 sage: J = RealSymmetricEJA(2)
505 Finite family {0: e0, 1: e1, 2: e2}
506 sage: J.natural_basis()
508 [1 0] [ 0 1/2*sqrt2] [0 0]
509 [0 0], [1/2*sqrt2 0], [0 1]
514 sage: J = JordanSpinEJA(2)
516 Finite family {0: e0, 1: e1}
517 sage: J.natural_basis()
524 if self
._natural
_basis
is None:
525 M
= self
.natural_basis_space()
526 return tuple( M(b
.to_vector()) for b
in self
.basis() )
528 return self
._natural
_basis
531 def natural_basis_space(self
):
533 Return the matrix space in which this algebra's natural basis
536 if self
._natural
_basis
is None or len(self
._natural
_basis
) == 0:
537 return MatrixSpace(self
.base_ring(), self
.dimension(), 1)
539 return self
._natural
_basis
[0].matrix_space()
543 def natural_inner_product(X
,Y
):
545 Compute the inner product of two naturally-represented elements.
547 For example in the real symmetric matrix EJA, this will compute
548 the trace inner-product of two n-by-n symmetric matrices. The
549 default should work for the real cartesian product EJA, the
550 Jordan spin EJA, and the real symmetric matrices. The others
551 will have to be overridden.
553 return (X
.conjugate_transpose()*Y
).trace()
559 Return the unit element of this algebra.
563 sage: from mjo.eja.eja_algebra import (RealCartesianProductEJA,
568 sage: J = RealCartesianProductEJA(5)
570 e0 + e1 + e2 + e3 + e4
574 The identity element acts like the identity::
576 sage: set_random_seed()
577 sage: J = random_eja()
578 sage: x = J.random_element()
579 sage: J.one()*x == x and x*J.one() == x
582 The matrix of the unit element's operator is the identity::
584 sage: set_random_seed()
585 sage: J = random_eja()
586 sage: actual = J.one().operator().matrix()
587 sage: expected = matrix.identity(J.base_ring(), J.dimension())
588 sage: actual == expected
592 # We can brute-force compute the matrices of the operators
593 # that correspond to the basis elements of this algebra.
594 # If some linear combination of those basis elements is the
595 # algebra identity, then the same linear combination of
596 # their matrices has to be the identity matrix.
598 # Of course, matrices aren't vectors in sage, so we have to
599 # appeal to the "long vectors" isometry.
600 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
602 # Now we use basis linear algebra to find the coefficients,
603 # of the matrices-as-vectors-linear-combination, which should
604 # work for the original algebra basis too.
605 A
= matrix
.column(self
.base_ring(), oper_vecs
)
607 # We used the isometry on the left-hand side already, but we
608 # still need to do it for the right-hand side. Recall that we
609 # wanted something that summed to the identity matrix.
610 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
612 # Now if there's an identity element in the algebra, this should work.
613 coeffs
= A
.solve_right(b
)
614 return self
.linear_combination(zip(self
.gens(), coeffs
))
617 def random_elements(self
, count
):
619 Return ``count`` random elements as a tuple.
623 sage: from mjo.eja.eja_algebra import JordanSpinEJA
627 sage: J = JordanSpinEJA(3)
628 sage: x,y,z = J.random_elements(3)
629 sage: all( [ x in J, y in J, z in J ])
631 sage: len( J.random_elements(10) ) == 10
635 return tuple( self
.random_element() for idx
in xrange(count
) )
640 Return the rank of this EJA.
644 The author knows of no algorithm to compute the rank of an EJA
645 where only the multiplication table is known. In lieu of one, we
646 require the rank to be specified when the algebra is created,
647 and simply pass along that number here.
651 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
652 ....: RealSymmetricEJA,
653 ....: ComplexHermitianEJA,
654 ....: QuaternionHermitianEJA,
659 The rank of the Jordan spin algebra is always two::
661 sage: JordanSpinEJA(2).rank()
663 sage: JordanSpinEJA(3).rank()
665 sage: JordanSpinEJA(4).rank()
668 The rank of the `n`-by-`n` Hermitian real, complex, or
669 quaternion matrices is `n`::
671 sage: RealSymmetricEJA(4).rank()
673 sage: ComplexHermitianEJA(3).rank()
675 sage: QuaternionHermitianEJA(2).rank()
680 Ensure that every EJA that we know how to construct has a
681 positive integer rank, unless the algebra is trivial in
682 which case its rank will be zero::
684 sage: set_random_seed()
685 sage: J = random_eja()
689 sage: r > 0 or (r == 0 and J.is_trivial())
696 def vector_space(self
):
698 Return the vector space that underlies this algebra.
702 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
706 sage: J = RealSymmetricEJA(2)
707 sage: J.vector_space()
708 Vector space of dimension 3 over...
711 return self
.zero().to_vector().parent().ambient_vector_space()
714 Element
= FiniteDimensionalEuclideanJordanAlgebraElement
717 class KnownRankEJA(object):
719 A class for algebras that we actually know we can construct. The
720 main issue is that, for most of our methods to make sense, we need
721 to know the rank of our algebra. Thus we can't simply generate a
722 "random" algebra, or even check that a given basis and product
723 satisfy the axioms; because even if everything looks OK, we wouldn't
724 know the rank we need to actuallty build the thing.
726 Not really a subclass of FDEJA because doing that causes method
727 resolution errors, e.g.
729 TypeError: Error when calling the metaclass bases
730 Cannot create a consistent method resolution
731 order (MRO) for bases FiniteDimensionalEuclideanJordanAlgebra,
736 def _max_test_case_size():
738 Return an integer "size" that is an upper bound on the size of
739 this algebra when it is used in a random test
740 case. Unfortunately, the term "size" is quite vague -- when
741 dealing with `R^n` under either the Hadamard or Jordan spin
742 product, the "size" refers to the dimension `n`. When dealing
743 with a matrix algebra (real symmetric or complex/quaternion
744 Hermitian), it refers to the size of the matrix, which is
745 far less than the dimension of the underlying vector space.
747 We default to five in this class, which is safe in `R^n`. The
748 matrix algebra subclasses (or any class where the "size" is
749 interpreted to be far less than the dimension) should override
750 with a smaller number.
755 def random_instance(cls
, field
=QQ
, **kwargs
):
757 Return a random instance of this type of algebra.
759 Beware, this will crash for "most instances" because the
760 constructor below looks wrong.
762 n
= ZZ
.random_element(cls
._max
_test
_case
_size
()) + 1
763 return cls(n
, field
, **kwargs
)
766 class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra
,
769 Return the Euclidean Jordan Algebra corresponding to the set
770 `R^n` under the Hadamard product.
772 Note: this is nothing more than the Cartesian product of ``n``
773 copies of the spin algebra. Once Cartesian product algebras
774 are implemented, this can go.
778 sage: from mjo.eja.eja_algebra import RealCartesianProductEJA
782 This multiplication table can be verified by hand::
784 sage: J = RealCartesianProductEJA(3)
785 sage: e0,e1,e2 = J.gens()
801 We can change the generator prefix::
803 sage: RealCartesianProductEJA(3, prefix='r').gens()
807 def __init__(self
, n
, field
=QQ
, **kwargs
):
808 V
= VectorSpace(field
, n
)
809 mult_table
= [ [ V
.gen(i
)*(i
== j
) for j
in xrange(n
) ]
812 fdeja
= super(RealCartesianProductEJA
, self
)
813 return fdeja
.__init
__(field
, mult_table
, rank
=n
, **kwargs
)
815 def inner_product(self
, x
, y
):
817 Faster to reimplement than to use natural representations.
821 sage: from mjo.eja.eja_algebra import RealCartesianProductEJA
825 Ensure that this is the usual inner product for the algebras
828 sage: set_random_seed()
829 sage: J = RealCartesianProductEJA.random_instance()
830 sage: x,y = J.random_elements(2)
831 sage: X = x.natural_representation()
832 sage: Y = y.natural_representation()
833 sage: x.inner_product(y) == J.natural_inner_product(X,Y)
837 return x
.to_vector().inner_product(y
.to_vector())
840 def random_eja(field
=QQ
):
842 Return a "random" finite-dimensional Euclidean Jordan Algebra.
846 For now, we choose a random natural number ``n`` (greater than zero)
847 and then give you back one of the following:
849 * The cartesian product of the rational numbers ``n`` times; this is
850 ``QQ^n`` with the Hadamard product.
852 * The Jordan spin algebra on ``QQ^n``.
854 * The ``n``-by-``n`` rational symmetric matrices with the symmetric
857 * The ``n``-by-``n`` complex-rational Hermitian matrices embedded
858 in the space of ``2n``-by-``2n`` real symmetric matrices.
860 * The ``n``-by-``n`` quaternion-rational Hermitian matrices embedded
861 in the space of ``4n``-by-``4n`` real symmetric matrices.
863 Later this might be extended to return Cartesian products of the
868 sage: from mjo.eja.eja_algebra import random_eja
873 Euclidean Jordan algebra of dimension...
876 classname
= choice(KnownRankEJA
.__subclasses
__())
877 return classname
.random_instance(field
=field
)
884 class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra
):
886 def _max_test_case_size():
887 # Play it safe, since this will be squared and the underlying
888 # field can have dimension 4 (quaternions) too.
891 def __init__(self
, field
, basis
, rank
, normalize_basis
=True, **kwargs
):
893 Compared to the superclass constructor, we take a basis instead of
894 a multiplication table because the latter can be computed in terms
895 of the former when the product is known (like it is here).
897 # Used in this class's fast _charpoly_coeff() override.
898 self
._basis
_normalizers
= None
900 # We're going to loop through this a few times, so now's a good
901 # time to ensure that it isn't a generator expression.
904 if rank
> 1 and normalize_basis
:
905 # We'll need sqrt(2) to normalize the basis, and this
906 # winds up in the multiplication table, so the whole
907 # algebra needs to be over the field extension.
908 R
= PolynomialRing(field
, 'z')
911 if p
.is_irreducible():
912 field
= field
.extension(p
, 'sqrt2', embedding
=RLF(2).sqrt())
913 basis
= tuple( s
.change_ring(field
) for s
in basis
)
914 self
._basis
_normalizers
= tuple(
915 ~
(self
.natural_inner_product(s
,s
).sqrt()) for s
in basis
)
916 basis
= tuple(s
*c
for (s
,c
) in izip(basis
,self
._basis
_normalizers
))
918 Qs
= self
.multiplication_table_from_matrix_basis(basis
)
920 fdeja
= super(MatrixEuclideanJordanAlgebra
, self
)
921 return fdeja
.__init
__(field
,
929 def _charpoly_coeff(self
, i
):
931 Override the parent method with something that tries to compute
932 over a faster (non-extension) field.
934 if self
._basis
_normalizers
is None:
935 # We didn't normalize, so assume that the basis we started
936 # with had entries in a nice field.
937 return super(MatrixEuclideanJordanAlgebra
, self
)._charpoly
_coeff
(i
)
939 basis
= ( (b
/n
) for (b
,n
) in izip(self
.natural_basis(),
940 self
._basis
_normalizers
) )
942 # Do this over the rationals and convert back at the end.
943 J
= MatrixEuclideanJordanAlgebra(QQ
,
946 normalize_basis
=False)
947 (_
,x
,_
,_
) = J
._charpoly
_matrix
_system
()
948 p
= J
._charpoly
_coeff
(i
)
949 # p might be missing some vars, have to substitute "optionally"
950 pairs
= izip(x
.base_ring().gens(), self
._basis
_normalizers
)
951 substitutions
= { v: v*c for (v,c) in pairs }
952 result
= p
.subs(substitutions
)
954 # The result of "subs" can be either a coefficient-ring
955 # element or a polynomial. Gotta handle both cases.
957 return self
.base_ring()(result
)
959 return result
.change_ring(self
.base_ring())
963 def multiplication_table_from_matrix_basis(basis
):
965 At least three of the five simple Euclidean Jordan algebras have the
966 symmetric multiplication (A,B) |-> (AB + BA)/2, where the
967 multiplication on the right is matrix multiplication. Given a basis
968 for the underlying matrix space, this function returns a
969 multiplication table (obtained by looping through the basis
970 elements) for an algebra of those matrices.
972 # In S^2, for example, we nominally have four coordinates even
973 # though the space is of dimension three only. The vector space V
974 # is supposed to hold the entire long vector, and the subspace W
975 # of V will be spanned by the vectors that arise from symmetric
976 # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
977 field
= basis
[0].base_ring()
978 dimension
= basis
[0].nrows()
980 V
= VectorSpace(field
, dimension
**2)
981 W
= V
.span_of_basis( _mat2vec(s
) for s
in basis
)
983 mult_table
= [[W
.zero() for j
in xrange(n
)] for i
in xrange(n
)]
986 mat_entry
= (basis
[i
]*basis
[j
] + basis
[j
]*basis
[i
])/2
987 mult_table
[i
][j
] = W
.coordinate_vector(_mat2vec(mat_entry
))
995 Embed the matrix ``M`` into a space of real matrices.
997 The matrix ``M`` can have entries in any field at the moment:
998 the real numbers, complex numbers, or quaternions. And although
999 they are not a field, we can probably support octonions at some
1000 point, too. This function returns a real matrix that "acts like"
1001 the original with respect to matrix multiplication; i.e.
1003 real_embed(M*N) = real_embed(M)*real_embed(N)
1006 raise NotImplementedError
1010 def real_unembed(M
):
1012 The inverse of :meth:`real_embed`.
1014 raise NotImplementedError
1018 def natural_inner_product(cls
,X
,Y
):
1019 Xu
= cls
.real_unembed(X
)
1020 Yu
= cls
.real_unembed(Y
)
1021 tr
= (Xu
*Yu
).trace()
1024 # It's real already.
1027 # Otherwise, try the thing that works for complex numbers; and
1028 # if that doesn't work, the thing that works for quaternions.
1030 return tr
.vector()[0] # real part, imag part is index 1
1031 except AttributeError:
1032 # A quaternions doesn't have a vector() method, but does
1033 # have coefficient_tuple() method that returns the
1034 # coefficients of 1, i, j, and k -- in that order.
1035 return tr
.coefficient_tuple()[0]
1038 class RealMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1042 The identity function, for embedding real matrices into real
1048 def real_unembed(M
):
1050 The identity function, for unembedding real matrices from real
1056 class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra
, KnownRankEJA
):
1058 The rank-n simple EJA consisting of real symmetric n-by-n
1059 matrices, the usual symmetric Jordan product, and the trace inner
1060 product. It has dimension `(n^2 + n)/2` over the reals.
1064 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1068 sage: J = RealSymmetricEJA(2)
1069 sage: e0, e1, e2 = J.gens()
1077 In theory, our "field" can be any subfield of the reals::
1079 sage: RealSymmetricEJA(2, AA)
1080 Euclidean Jordan algebra of dimension 3 over Algebraic Real Field
1081 sage: RealSymmetricEJA(2, RR)
1082 Euclidean Jordan algebra of dimension 3 over Real Field with
1083 53 bits of precision
1087 The dimension of this algebra is `(n^2 + n) / 2`::
1089 sage: set_random_seed()
1090 sage: n_max = RealSymmetricEJA._max_test_case_size()
1091 sage: n = ZZ.random_element(1, n_max)
1092 sage: J = RealSymmetricEJA(n)
1093 sage: J.dimension() == (n^2 + n)/2
1096 The Jordan multiplication is what we think it is::
1098 sage: set_random_seed()
1099 sage: J = RealSymmetricEJA.random_instance()
1100 sage: x,y = J.random_elements(2)
1101 sage: actual = (x*y).natural_representation()
1102 sage: X = x.natural_representation()
1103 sage: Y = y.natural_representation()
1104 sage: expected = (X*Y + Y*X)/2
1105 sage: actual == expected
1107 sage: J(expected) == x*y
1110 We can change the generator prefix::
1112 sage: RealSymmetricEJA(3, prefix='q').gens()
1113 (q0, q1, q2, q3, q4, q5)
1115 Our natural basis is normalized with respect to the natural inner
1116 product unless we specify otherwise::
1118 sage: set_random_seed()
1119 sage: J = RealSymmetricEJA.random_instance()
1120 sage: all( b.norm() == 1 for b in J.gens() )
1123 Since our natural basis is normalized with respect to the natural
1124 inner product, and since we know that this algebra is an EJA, any
1125 left-multiplication operator's matrix will be symmetric because
1126 natural->EJA basis representation is an isometry and within the EJA
1127 the operator is self-adjoint by the Jordan axiom::
1129 sage: set_random_seed()
1130 sage: x = RealSymmetricEJA.random_instance().random_element()
1131 sage: x.operator().matrix().is_symmetric()
1136 def _denormalized_basis(cls
, n
, field
):
1138 Return a basis for the space of real symmetric n-by-n matrices.
1142 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1146 sage: set_random_seed()
1147 sage: n = ZZ.random_element(1,5)
1148 sage: B = RealSymmetricEJA._denormalized_basis(n,QQ)
1149 sage: all( M.is_symmetric() for M in B)
1153 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1157 for j
in xrange(i
+1):
1158 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
1162 Sij
= Eij
+ Eij
.transpose()
1168 def _max_test_case_size():
1169 return 4 # Dimension 10
1172 def __init__(self
, n
, field
=QQ
, **kwargs
):
1173 basis
= self
._denormalized
_basis
(n
, field
)
1174 super(RealSymmetricEJA
, self
).__init
__(field
, basis
, n
, **kwargs
)
1177 class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1181 Embed the n-by-n complex matrix ``M`` into the space of real
1182 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1183 bi` to the block matrix ``[[a,b],[-b,a]]``.
1187 sage: from mjo.eja.eja_algebra import \
1188 ....: ComplexMatrixEuclideanJordanAlgebra
1192 sage: F = QuadraticField(-1, 'i')
1193 sage: x1 = F(4 - 2*i)
1194 sage: x2 = F(1 + 2*i)
1197 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1198 sage: ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1207 Embedding is a homomorphism (isomorphism, in fact)::
1209 sage: set_random_seed()
1210 sage: n_max = ComplexMatrixEuclideanJordanAlgebra._max_test_case_size()
1211 sage: n = ZZ.random_element(n_max)
1212 sage: F = QuadraticField(-1, 'i')
1213 sage: X = random_matrix(F, n)
1214 sage: Y = random_matrix(F, n)
1215 sage: Xe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X)
1216 sage: Ye = ComplexMatrixEuclideanJordanAlgebra.real_embed(Y)
1217 sage: XYe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1224 raise ValueError("the matrix 'M' must be square")
1226 # We don't need any adjoined elements...
1227 field
= M
.base_ring().base_ring()
1231 a
= z
.list()[0] # real part, I guess
1232 b
= z
.list()[1] # imag part, I guess
1233 blocks
.append(matrix(field
, 2, [[a
,b
],[-b
,a
]]))
1235 return matrix
.block(field
, n
, blocks
)
1239 def real_unembed(M
):
1241 The inverse of _embed_complex_matrix().
1245 sage: from mjo.eja.eja_algebra import \
1246 ....: ComplexMatrixEuclideanJordanAlgebra
1250 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1251 ....: [-2, 1, -4, 3],
1252 ....: [ 9, 10, 11, 12],
1253 ....: [-10, 9, -12, 11] ])
1254 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(A)
1256 [ 10*i + 9 12*i + 11]
1260 Unembedding is the inverse of embedding::
1262 sage: set_random_seed()
1263 sage: F = QuadraticField(-1, 'i')
1264 sage: M = random_matrix(F, 3)
1265 sage: Me = ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1266 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1272 raise ValueError("the matrix 'M' must be square")
1273 if not n
.mod(2).is_zero():
1274 raise ValueError("the matrix 'M' must be a complex embedding")
1276 # If "M" was normalized, its base ring might have roots
1277 # adjoined and they can stick around after unembedding.
1278 field
= M
.base_ring()
1279 R
= PolynomialRing(field
, 'z')
1281 F
= field
.extension(z
**2 + 1, 'i', embedding
=CLF(-1).sqrt())
1284 # Go top-left to bottom-right (reading order), converting every
1285 # 2-by-2 block we see to a single complex element.
1287 for k
in xrange(n
/2):
1288 for j
in xrange(n
/2):
1289 submat
= M
[2*k
:2*k
+2,2*j
:2*j
+2]
1290 if submat
[0,0] != submat
[1,1]:
1291 raise ValueError('bad on-diagonal submatrix')
1292 if submat
[0,1] != -submat
[1,0]:
1293 raise ValueError('bad off-diagonal submatrix')
1294 z
= submat
[0,0] + submat
[0,1]*i
1297 return matrix(F
, n
/2, elements
)
1301 def natural_inner_product(cls
,X
,Y
):
1303 Compute a natural inner product in this algebra directly from
1308 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1312 This gives the same answer as the slow, default method implemented
1313 in :class:`MatrixEuclideanJordanAlgebra`::
1315 sage: set_random_seed()
1316 sage: J = ComplexHermitianEJA.random_instance()
1317 sage: x,y = J.random_elements(2)
1318 sage: Xe = x.natural_representation()
1319 sage: Ye = y.natural_representation()
1320 sage: X = ComplexHermitianEJA.real_unembed(Xe)
1321 sage: Y = ComplexHermitianEJA.real_unembed(Ye)
1322 sage: expected = (X*Y).trace().vector()[0]
1323 sage: actual = ComplexHermitianEJA.natural_inner_product(Xe,Ye)
1324 sage: actual == expected
1328 return RealMatrixEuclideanJordanAlgebra
.natural_inner_product(X
,Y
)/2
1331 class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra
, KnownRankEJA
):
1333 The rank-n simple EJA consisting of complex Hermitian n-by-n
1334 matrices over the real numbers, the usual symmetric Jordan product,
1335 and the real-part-of-trace inner product. It has dimension `n^2` over
1340 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1344 In theory, our "field" can be any subfield of the reals::
1346 sage: ComplexHermitianEJA(2, AA)
1347 Euclidean Jordan algebra of dimension 4 over Algebraic Real Field
1348 sage: ComplexHermitianEJA(2, RR)
1349 Euclidean Jordan algebra of dimension 4 over Real Field with
1350 53 bits of precision
1354 The dimension of this algebra is `n^2`::
1356 sage: set_random_seed()
1357 sage: n_max = ComplexHermitianEJA._max_test_case_size()
1358 sage: n = ZZ.random_element(1, n_max)
1359 sage: J = ComplexHermitianEJA(n)
1360 sage: J.dimension() == n^2
1363 The Jordan multiplication is what we think it is::
1365 sage: set_random_seed()
1366 sage: J = ComplexHermitianEJA.random_instance()
1367 sage: x,y = J.random_elements(2)
1368 sage: actual = (x*y).natural_representation()
1369 sage: X = x.natural_representation()
1370 sage: Y = y.natural_representation()
1371 sage: expected = (X*Y + Y*X)/2
1372 sage: actual == expected
1374 sage: J(expected) == x*y
1377 We can change the generator prefix::
1379 sage: ComplexHermitianEJA(2, prefix='z').gens()
1382 Our natural basis is normalized with respect to the natural inner
1383 product unless we specify otherwise::
1385 sage: set_random_seed()
1386 sage: J = ComplexHermitianEJA.random_instance()
1387 sage: all( b.norm() == 1 for b in J.gens() )
1390 Since our natural basis is normalized with respect to the natural
1391 inner product, and since we know that this algebra is an EJA, any
1392 left-multiplication operator's matrix will be symmetric because
1393 natural->EJA basis representation is an isometry and within the EJA
1394 the operator is self-adjoint by the Jordan axiom::
1396 sage: set_random_seed()
1397 sage: x = ComplexHermitianEJA.random_instance().random_element()
1398 sage: x.operator().matrix().is_symmetric()
1404 def _denormalized_basis(cls
, n
, field
):
1406 Returns a basis for the space of complex Hermitian n-by-n matrices.
1408 Why do we embed these? Basically, because all of numerical linear
1409 algebra assumes that you're working with vectors consisting of `n`
1410 entries from a field and scalars from the same field. There's no way
1411 to tell SageMath that (for example) the vectors contain complex
1412 numbers, while the scalar field is real.
1416 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1420 sage: set_random_seed()
1421 sage: n = ZZ.random_element(1,5)
1422 sage: field = QuadraticField(2, 'sqrt2')
1423 sage: B = ComplexHermitianEJA._denormalized_basis(n, field)
1424 sage: all( M.is_symmetric() for M in B)
1428 R
= PolynomialRing(field
, 'z')
1430 F
= field
.extension(z
**2 + 1, 'I')
1433 # This is like the symmetric case, but we need to be careful:
1435 # * We want conjugate-symmetry, not just symmetry.
1436 # * The diagonal will (as a result) be real.
1440 for j
in xrange(i
+1):
1441 Eij
= matrix(F
, n
, lambda k
,l
: k
==i
and l
==j
)
1443 Sij
= cls
.real_embed(Eij
)
1446 # The second one has a minus because it's conjugated.
1447 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
1449 Sij_imag
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
1452 # Since we embedded these, we can drop back to the "field" that we
1453 # started with instead of the complex extension "F".
1454 return ( s
.change_ring(field
) for s
in S
)
1457 def __init__(self
, n
, field
=QQ
, **kwargs
):
1458 basis
= self
._denormalized
_basis
(n
,field
)
1459 super(ComplexHermitianEJA
,self
).__init
__(field
, basis
, n
, **kwargs
)
1462 class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1466 Embed the n-by-n quaternion matrix ``M`` into the space of real
1467 matrices of size 4n-by-4n by first sending each quaternion entry `z
1468 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
1469 c+di],[-c + di, a-bi]]`, and then embedding those into a real
1474 sage: from mjo.eja.eja_algebra import \
1475 ....: QuaternionMatrixEuclideanJordanAlgebra
1479 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1480 sage: i,j,k = Q.gens()
1481 sage: x = 1 + 2*i + 3*j + 4*k
1482 sage: M = matrix(Q, 1, [[x]])
1483 sage: QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1489 Embedding is a homomorphism (isomorphism, in fact)::
1491 sage: set_random_seed()
1492 sage: n_max = QuaternionMatrixEuclideanJordanAlgebra._max_test_case_size()
1493 sage: n = ZZ.random_element(n_max)
1494 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1495 sage: X = random_matrix(Q, n)
1496 sage: Y = random_matrix(Q, n)
1497 sage: Xe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X)
1498 sage: Ye = QuaternionMatrixEuclideanJordanAlgebra.real_embed(Y)
1499 sage: XYe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1504 quaternions
= M
.base_ring()
1507 raise ValueError("the matrix 'M' must be square")
1509 F
= QuadraticField(-1, 'i')
1514 t
= z
.coefficient_tuple()
1519 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
1520 [-c
+ d
*i
, a
- b
*i
]])
1521 realM
= ComplexMatrixEuclideanJordanAlgebra
.real_embed(cplxM
)
1522 blocks
.append(realM
)
1524 # We should have real entries by now, so use the realest field
1525 # we've got for the return value.
1526 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
1531 def real_unembed(M
):
1533 The inverse of _embed_quaternion_matrix().
1537 sage: from mjo.eja.eja_algebra import \
1538 ....: QuaternionMatrixEuclideanJordanAlgebra
1542 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
1543 ....: [-2, 1, -4, 3],
1544 ....: [-3, 4, 1, -2],
1545 ....: [-4, -3, 2, 1]])
1546 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(M)
1547 [1 + 2*i + 3*j + 4*k]
1551 Unembedding is the inverse of embedding::
1553 sage: set_random_seed()
1554 sage: Q = QuaternionAlgebra(QQ, -1, -1)
1555 sage: M = random_matrix(Q, 3)
1556 sage: Me = QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1557 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1563 raise ValueError("the matrix 'M' must be square")
1564 if not n
.mod(4).is_zero():
1565 raise ValueError("the matrix 'M' must be a quaternion embedding")
1567 # Use the base ring of the matrix to ensure that its entries can be
1568 # multiplied by elements of the quaternion algebra.
1569 field
= M
.base_ring()
1570 Q
= QuaternionAlgebra(field
,-1,-1)
1573 # Go top-left to bottom-right (reading order), converting every
1574 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
1577 for l
in xrange(n
/4):
1578 for m
in xrange(n
/4):
1579 submat
= ComplexMatrixEuclideanJordanAlgebra
.real_unembed(
1580 M
[4*l
:4*l
+4,4*m
:4*m
+4] )
1581 if submat
[0,0] != submat
[1,1].conjugate():
1582 raise ValueError('bad on-diagonal submatrix')
1583 if submat
[0,1] != -submat
[1,0].conjugate():
1584 raise ValueError('bad off-diagonal submatrix')
1585 z
= submat
[0,0].vector()[0] # real part
1586 z
+= submat
[0,0].vector()[1]*i
# imag part
1587 z
+= submat
[0,1].vector()[0]*j
# real part
1588 z
+= submat
[0,1].vector()[1]*k
# imag part
1591 return matrix(Q
, n
/4, elements
)
1595 def natural_inner_product(cls
,X
,Y
):
1597 Compute a natural inner product in this algebra directly from
1602 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1606 This gives the same answer as the slow, default method implemented
1607 in :class:`MatrixEuclideanJordanAlgebra`::
1609 sage: set_random_seed()
1610 sage: J = QuaternionHermitianEJA.random_instance()
1611 sage: x,y = J.random_elements(2)
1612 sage: Xe = x.natural_representation()
1613 sage: Ye = y.natural_representation()
1614 sage: X = QuaternionHermitianEJA.real_unembed(Xe)
1615 sage: Y = QuaternionHermitianEJA.real_unembed(Ye)
1616 sage: expected = (X*Y).trace().coefficient_tuple()[0]
1617 sage: actual = QuaternionHermitianEJA.natural_inner_product(Xe,Ye)
1618 sage: actual == expected
1622 return RealMatrixEuclideanJordanAlgebra
.natural_inner_product(X
,Y
)/4
1625 class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra
,
1628 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
1629 matrices, the usual symmetric Jordan product, and the
1630 real-part-of-trace inner product. It has dimension `2n^2 - n` over
1635 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1639 In theory, our "field" can be any subfield of the reals::
1641 sage: QuaternionHermitianEJA(2, AA)
1642 Euclidean Jordan algebra of dimension 6 over Algebraic Real Field
1643 sage: QuaternionHermitianEJA(2, RR)
1644 Euclidean Jordan algebra of dimension 6 over Real Field with
1645 53 bits of precision
1649 The dimension of this algebra is `2*n^2 - n`::
1651 sage: set_random_seed()
1652 sage: n_max = QuaternionHermitianEJA._max_test_case_size()
1653 sage: n = ZZ.random_element(1, n_max)
1654 sage: J = QuaternionHermitianEJA(n)
1655 sage: J.dimension() == 2*(n^2) - n
1658 The Jordan multiplication is what we think it is::
1660 sage: set_random_seed()
1661 sage: J = QuaternionHermitianEJA.random_instance()
1662 sage: x,y = J.random_elements(2)
1663 sage: actual = (x*y).natural_representation()
1664 sage: X = x.natural_representation()
1665 sage: Y = y.natural_representation()
1666 sage: expected = (X*Y + Y*X)/2
1667 sage: actual == expected
1669 sage: J(expected) == x*y
1672 We can change the generator prefix::
1674 sage: QuaternionHermitianEJA(2, prefix='a').gens()
1675 (a0, a1, a2, a3, a4, a5)
1677 Our natural basis is normalized with respect to the natural inner
1678 product unless we specify otherwise::
1680 sage: set_random_seed()
1681 sage: J = QuaternionHermitianEJA.random_instance()
1682 sage: all( b.norm() == 1 for b in J.gens() )
1685 Since our natural basis is normalized with respect to the natural
1686 inner product, and since we know that this algebra is an EJA, any
1687 left-multiplication operator's matrix will be symmetric because
1688 natural->EJA basis representation is an isometry and within the EJA
1689 the operator is self-adjoint by the Jordan axiom::
1691 sage: set_random_seed()
1692 sage: x = QuaternionHermitianEJA.random_instance().random_element()
1693 sage: x.operator().matrix().is_symmetric()
1698 def _denormalized_basis(cls
, n
, field
):
1700 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
1702 Why do we embed these? Basically, because all of numerical
1703 linear algebra assumes that you're working with vectors consisting
1704 of `n` entries from a field and scalars from the same field. There's
1705 no way to tell SageMath that (for example) the vectors contain
1706 complex numbers, while the scalar field is real.
1710 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1714 sage: set_random_seed()
1715 sage: n = ZZ.random_element(1,5)
1716 sage: B = QuaternionHermitianEJA._denormalized_basis(n,QQ)
1717 sage: all( M.is_symmetric() for M in B )
1721 Q
= QuaternionAlgebra(QQ
,-1,-1)
1724 # This is like the symmetric case, but we need to be careful:
1726 # * We want conjugate-symmetry, not just symmetry.
1727 # * The diagonal will (as a result) be real.
1731 for j
in xrange(i
+1):
1732 Eij
= matrix(Q
, n
, lambda k
,l
: k
==i
and l
==j
)
1734 Sij
= cls
.real_embed(Eij
)
1737 # The second, third, and fourth ones have a minus
1738 # because they're conjugated.
1739 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
1741 Sij_I
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
1743 Sij_J
= cls
.real_embed(J
*Eij
- J
*Eij
.transpose())
1745 Sij_K
= cls
.real_embed(K
*Eij
- K
*Eij
.transpose())
1748 # Since we embedded these, we can drop back to the "field" that we
1749 # started with instead of the quaternion algebra "Q".
1750 return ( s
.change_ring(field
) for s
in S
)
1753 def __init__(self
, n
, field
=QQ
, **kwargs
):
1754 basis
= self
._denormalized
_basis
(n
,field
)
1755 super(QuaternionHermitianEJA
,self
).__init
__(field
, basis
, n
, **kwargs
)
1758 class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra
, KnownRankEJA
):
1760 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
1761 with the usual inner product and jordan product ``x*y =
1762 (<x_bar,y_bar>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
1767 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1771 This multiplication table can be verified by hand::
1773 sage: J = JordanSpinEJA(4)
1774 sage: e0,e1,e2,e3 = J.gens()
1790 We can change the generator prefix::
1792 sage: JordanSpinEJA(2, prefix='B').gens()
1796 def __init__(self
, n
, field
=QQ
, **kwargs
):
1797 V
= VectorSpace(field
, n
)
1798 mult_table
= [[V
.zero() for j
in xrange(n
)] for i
in xrange(n
)]
1808 z0
= x
.inner_product(y
)
1809 zbar
= y0
*xbar
+ x0
*ybar
1810 z
= V([z0
] + zbar
.list())
1811 mult_table
[i
][j
] = z
1813 # The rank of the spin algebra is two, unless we're in a
1814 # one-dimensional ambient space (because the rank is bounded by
1815 # the ambient dimension).
1816 fdeja
= super(JordanSpinEJA
, self
)
1817 return fdeja
.__init
__(field
, mult_table
, rank
=min(n
,2), **kwargs
)
1819 def inner_product(self
, x
, y
):
1821 Faster to reimplement than to use natural representations.
1825 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1829 Ensure that this is the usual inner product for the algebras
1832 sage: set_random_seed()
1833 sage: J = JordanSpinEJA.random_instance()
1834 sage: x,y = J.random_elements(2)
1835 sage: X = x.natural_representation()
1836 sage: Y = y.natural_representation()
1837 sage: x.inner_product(y) == J.natural_inner_product(X,Y)
1841 return x
.to_vector().inner_product(y
.to_vector())