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
.lazy_import
import lazy_import
17 from sage
.misc
.prandom
import choice
18 from sage
.misc
.table
import table
19 from sage
.modules
.free_module
import FreeModule
, VectorSpace
20 from sage
.rings
.all
import (ZZ
, QQ
, RR
, RLF
, CLF
,
23 from mjo
.eja
.eja_element
import FiniteDimensionalEuclideanJordanAlgebraElement
24 lazy_import('mjo.eja.eja_subalgebra',
25 'FiniteDimensionalEuclideanJordanSubalgebra')
26 from mjo
.eja
.eja_utils
import _mat2vec
28 class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule
):
29 # This is an ugly hack needed to prevent the category framework
30 # from implementing a coercion from our base ring (e.g. the
31 # rationals) into the algebra. First of all -- such a coercion is
32 # nonsense to begin with. But more importantly, it tries to do so
33 # in the category of rings, and since our algebras aren't
34 # associative they generally won't be rings.
35 _no_generic_basering_coercion
= True
48 sage: from mjo.eja.eja_algebra import (JordanSpinEJA, 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 The ``field`` we're given must be real::
64 sage: JordanSpinEJA(2,QQbar)
65 Traceback (most recent call last):
67 ValueError: field is not real
71 if not field
.is_subring(RR
):
72 # Note: this does return true for the real algebraic
73 # field, and any quadratic field where we've specified
75 raise ValueError('field is not real')
78 self
._natural
_basis
= natural_basis
81 category
= MagmaticAlgebras(field
).FiniteDimensional()
82 category
= category
.WithBasis().Unital()
84 fda
= super(FiniteDimensionalEuclideanJordanAlgebra
, self
)
86 range(len(mult_table
)),
89 self
.print_options(bracket
='')
91 # The multiplication table we're given is necessarily in terms
92 # of vectors, because we don't have an algebra yet for
93 # anything to be an element of. However, it's faster in the
94 # long run to have the multiplication table be in terms of
95 # algebra elements. We do this after calling the superclass
96 # constructor so that from_vector() knows what to do.
97 self
._multiplication
_table
= [ map(lambda x
: self
.from_vector(x
), ls
)
98 for ls
in mult_table
]
101 def _element_constructor_(self
, elt
):
103 Construct an element of this algebra from its natural
106 This gets called only after the parent element _call_ method
107 fails to find a coercion for the argument.
111 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
112 ....: RealCartesianProductEJA,
113 ....: RealSymmetricEJA)
117 The identity in `S^n` is converted to the identity in the EJA::
119 sage: J = RealSymmetricEJA(3)
120 sage: I = matrix.identity(QQ,3)
121 sage: J(I) == J.one()
124 This skew-symmetric matrix can't be represented in the EJA::
126 sage: J = RealSymmetricEJA(3)
127 sage: A = matrix(QQ,3, lambda i,j: i-j)
129 Traceback (most recent call last):
131 ArithmeticError: vector is not in free module
135 Ensure that we can convert any element of the two non-matrix
136 simple algebras (whose natural representations are their usual
137 vector representations) back and forth faithfully::
139 sage: set_random_seed()
140 sage: J = RealCartesianProductEJA.random_instance()
141 sage: x = J.random_element()
142 sage: J(x.to_vector().column()) == x
144 sage: J = JordanSpinEJA.random_instance()
145 sage: x = J.random_element()
146 sage: J(x.to_vector().column()) == x
151 # The superclass implementation of random_element()
152 # needs to be able to coerce "0" into the algebra.
155 natural_basis
= self
.natural_basis()
156 basis_space
= natural_basis
[0].matrix_space()
157 if elt
not in basis_space
:
158 raise ValueError("not a naturally-represented algebra element")
160 # Thanks for nothing! Matrix spaces aren't vector spaces in
161 # Sage, so we have to figure out its natural-basis coordinates
162 # ourselves. We use the basis space's ring instead of the
163 # element's ring because the basis space might be an algebraic
164 # closure whereas the base ring of the 3-by-3 identity matrix
165 # could be QQ instead of QQbar.
166 V
= VectorSpace(basis_space
.base_ring(), elt
.nrows()*elt
.ncols())
167 W
= V
.span_of_basis( _mat2vec(s
) for s
in natural_basis
)
168 coords
= W
.coordinate_vector(_mat2vec(elt
))
169 return self
.from_vector(coords
)
174 Return a string representation of ``self``.
178 sage: from mjo.eja.eja_algebra import JordanSpinEJA
182 Ensure that it says what we think it says::
184 sage: JordanSpinEJA(2, field=QQ)
185 Euclidean Jordan algebra of dimension 2 over Rational Field
186 sage: JordanSpinEJA(3, field=RDF)
187 Euclidean Jordan algebra of dimension 3 over Real Double Field
190 fmt
= "Euclidean Jordan algebra of dimension {} over {}"
191 return fmt
.format(self
.dimension(), self
.base_ring())
193 def product_on_basis(self
, i
, j
):
194 return self
._multiplication
_table
[i
][j
]
196 def _a_regular_element(self
):
198 Guess a regular element. Needed to compute the basis for our
199 characteristic polynomial coefficients.
203 sage: from mjo.eja.eja_algebra import random_eja
207 Ensure that this hacky method succeeds for every algebra that we
208 know how to construct::
210 sage: set_random_seed()
211 sage: J = random_eja()
212 sage: J._a_regular_element().is_regular()
217 z
= self
.sum( (i
+1)*gs
[i
] for i
in range(len(gs
)) )
218 if not z
.is_regular():
219 raise ValueError("don't know a regular element")
224 def _charpoly_basis_space(self
):
226 Return the vector space spanned by the basis used in our
227 characteristic polynomial coefficients. This is used not only to
228 compute those coefficients, but also any time we need to
229 evaluate the coefficients (like when we compute the trace or
232 z
= self
._a
_regular
_element
()
233 # Don't use the parent vector space directly here in case this
234 # happens to be a subalgebra. In that case, we would be e.g.
235 # two-dimensional but span_of_basis() would expect three
237 V
= VectorSpace(self
.base_ring(), self
.vector_space().dimension())
238 basis
= [ (z
**k
).to_vector() for k
in range(self
.rank()) ]
239 V1
= V
.span_of_basis( basis
)
240 b
= (V1
.basis() + V1
.complement().basis())
241 return V
.span_of_basis(b
)
246 def _charpoly_coeff(self
, i
):
248 Return the coefficient polynomial "a_{i}" of this algebra's
249 general characteristic polynomial.
251 Having this be a separate cached method lets us compute and
252 store the trace/determinant (a_{r-1} and a_{0} respectively)
253 separate from the entire characteristic polynomial.
255 (A_of_x
, x
, xr
, detA
) = self
._charpoly
_matrix
_system
()
256 R
= A_of_x
.base_ring()
261 # Guaranteed by theory
264 # Danger: the in-place modification is done for performance
265 # reasons (reconstructing a matrix with huge polynomial
266 # entries is slow), but I don't know how cached_method works,
267 # so it's highly possible that we're modifying some global
268 # list variable by reference, here. In other words, you
269 # probably shouldn't call this method twice on the same
270 # algebra, at the same time, in two threads
271 Ai_orig
= A_of_x
.column(i
)
272 A_of_x
.set_column(i
,xr
)
273 numerator
= A_of_x
.det()
274 A_of_x
.set_column(i
,Ai_orig
)
276 # We're relying on the theory here to ensure that each a_i is
277 # indeed back in R, and the added negative signs are to make
278 # the whole charpoly expression sum to zero.
279 return R(-numerator
/detA
)
283 def _charpoly_matrix_system(self
):
285 Compute the matrix whose entries A_ij are polynomials in
286 X1,...,XN, the vector ``x`` of variables X1,...,XN, the vector
287 corresponding to `x^r` and the determinent of the matrix A =
288 [A_ij]. In other words, all of the fixed (cachable) data needed
289 to compute the coefficients of the characteristic polynomial.
294 # Turn my vector space into a module so that "vectors" can
295 # have multivatiate polynomial entries.
296 names
= tuple('X' + str(i
) for i
in range(1,n
+1))
297 R
= PolynomialRing(self
.base_ring(), names
)
299 # Using change_ring() on the parent's vector space doesn't work
300 # here because, in a subalgebra, that vector space has a basis
301 # and change_ring() tries to bring the basis along with it. And
302 # that doesn't work unless the new ring is a PID, which it usually
306 # Now let x = (X1,X2,...,Xn) be the vector whose entries are
310 # And figure out the "left multiplication by x" matrix in
313 monomial_matrices
= [ self
.monomial(i
).operator().matrix()
314 for i
in range(n
) ] # don't recompute these!
316 ek
= self
.monomial(k
).to_vector()
318 sum( x
[i
]*(monomial_matrices
[i
]*ek
)
319 for i
in range(n
) ) )
320 Lx
= matrix
.column(R
, lmbx_cols
)
322 # Now we can compute powers of x "symbolically"
323 x_powers
= [self
.one().to_vector(), x
]
324 for d
in range(2, r
+1):
325 x_powers
.append( Lx
*(x_powers
[-1]) )
327 idmat
= matrix
.identity(R
, n
)
329 W
= self
._charpoly
_basis
_space
()
330 W
= W
.change_ring(R
.fraction_field())
332 # Starting with the standard coordinates x = (X1,X2,...,Xn)
333 # and then converting the entries to W-coordinates allows us
334 # to pass in the standard coordinates to the charpoly and get
335 # back the right answer. Specifically, with x = (X1,X2,...,Xn),
338 # W.coordinates(x^2) eval'd at (standard z-coords)
342 # W-coords of (standard coords of x^2 eval'd at std-coords of z)
344 # We want the middle equivalent thing in our matrix, but use
345 # the first equivalent thing instead so that we can pass in
346 # standard coordinates.
347 x_powers
= [ W
.coordinate_vector(xp
) for xp
in x_powers
]
348 l2
= [idmat
.column(k
-1) for k
in range(r
+1, n
+1)]
349 A_of_x
= matrix
.column(R
, n
, (x_powers
[:r
] + l2
))
350 return (A_of_x
, x
, x_powers
[r
], A_of_x
.det())
354 def characteristic_polynomial(self
):
356 Return a characteristic polynomial that works for all elements
359 The resulting polynomial has `n+1` variables, where `n` is the
360 dimension of this algebra. The first `n` variables correspond to
361 the coordinates of an algebra element: when evaluated at the
362 coordinates of an algebra element with respect to a certain
363 basis, the result is a univariate polynomial (in the one
364 remaining variable ``t``), namely the characteristic polynomial
369 sage: from mjo.eja.eja_algebra import JordanSpinEJA, TrivialEJA
373 The characteristic polynomial in the spin algebra is given in
374 Alizadeh, Example 11.11::
376 sage: J = JordanSpinEJA(3)
377 sage: p = J.characteristic_polynomial(); p
378 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
379 sage: xvec = J.one().to_vector()
383 By definition, the characteristic polynomial is a monic
384 degree-zero polynomial in a rank-zero algebra. Note that
385 Cayley-Hamilton is indeed satisfied since the polynomial
386 ``1`` evaluates to the identity element of the algebra on
389 sage: J = TrivialEJA()
390 sage: J.characteristic_polynomial()
397 # The list of coefficient polynomials a_0, a_1, a_2, ..., a_n.
398 a
= [ self
._charpoly
_coeff
(i
) for i
in range(r
+1) ]
400 # We go to a bit of trouble here to reorder the
401 # indeterminates, so that it's easier to evaluate the
402 # characteristic polynomial at x's coordinates and get back
403 # something in terms of t, which is what we want.
405 S
= PolynomialRing(self
.base_ring(),'t')
407 S
= PolynomialRing(S
, R
.variable_names())
410 return sum( a
[k
]*(t
**k
) for k
in xrange(len(a
)) )
413 def inner_product(self
, x
, y
):
415 The inner product associated with this Euclidean Jordan algebra.
417 Defaults to the trace inner product, but can be overridden by
418 subclasses if they are sure that the necessary properties are
423 sage: from mjo.eja.eja_algebra import random_eja
427 Our inner product is "associative," which means the following for
428 a symmetric bilinear form::
430 sage: set_random_seed()
431 sage: J = random_eja()
432 sage: x,y,z = J.random_elements(3)
433 sage: (x*y).inner_product(z) == y.inner_product(x*z)
437 X
= x
.natural_representation()
438 Y
= y
.natural_representation()
439 return self
.natural_inner_product(X
,Y
)
442 def is_trivial(self
):
444 Return whether or not this algebra is trivial.
446 A trivial algebra contains only the zero element.
450 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
455 sage: J = ComplexHermitianEJA(3)
461 sage: J = TrivialEJA()
466 return self
.dimension() == 0
469 def multiplication_table(self
):
471 Return a visual representation of this algebra's multiplication
472 table (on basis elements).
476 sage: from mjo.eja.eja_algebra import JordanSpinEJA
480 sage: J = JordanSpinEJA(4)
481 sage: J.multiplication_table()
482 +----++----+----+----+----+
483 | * || e0 | e1 | e2 | e3 |
484 +====++====+====+====+====+
485 | e0 || e0 | e1 | e2 | e3 |
486 +----++----+----+----+----+
487 | e1 || e1 | e0 | 0 | 0 |
488 +----++----+----+----+----+
489 | e2 || e2 | 0 | e0 | 0 |
490 +----++----+----+----+----+
491 | e3 || e3 | 0 | 0 | e0 |
492 +----++----+----+----+----+
495 M
= list(self
._multiplication
_table
) # copy
496 for i
in xrange(len(M
)):
497 # M had better be "square"
498 M
[i
] = [self
.monomial(i
)] + M
[i
]
499 M
= [["*"] + list(self
.gens())] + M
500 return table(M
, header_row
=True, header_column
=True, frame
=True)
503 def natural_basis(self
):
505 Return a more-natural representation of this algebra's basis.
507 Every finite-dimensional Euclidean Jordan Algebra is a direct
508 sum of five simple algebras, four of which comprise Hermitian
509 matrices. This method returns the original "natural" basis
510 for our underlying vector space. (Typically, the natural basis
511 is used to construct the multiplication table in the first place.)
513 Note that this will always return a matrix. The standard basis
514 in `R^n` will be returned as `n`-by-`1` column matrices.
518 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
519 ....: RealSymmetricEJA)
523 sage: J = RealSymmetricEJA(2)
525 Finite family {0: e0, 1: e1, 2: e2}
526 sage: J.natural_basis()
528 [1 0] [ 0 1/2*sqrt2] [0 0]
529 [0 0], [1/2*sqrt2 0], [0 1]
534 sage: J = JordanSpinEJA(2)
536 Finite family {0: e0, 1: e1}
537 sage: J.natural_basis()
544 if self
._natural
_basis
is None:
545 M
= self
.natural_basis_space()
546 return tuple( M(b
.to_vector()) for b
in self
.basis() )
548 return self
._natural
_basis
551 def natural_basis_space(self
):
553 Return the matrix space in which this algebra's natural basis
556 if self
._natural
_basis
is None or len(self
._natural
_basis
) == 0:
557 return MatrixSpace(self
.base_ring(), self
.dimension(), 1)
559 return self
._natural
_basis
[0].matrix_space()
563 def natural_inner_product(X
,Y
):
565 Compute the inner product of two naturally-represented elements.
567 For example in the real symmetric matrix EJA, this will compute
568 the trace inner-product of two n-by-n symmetric matrices. The
569 default should work for the real cartesian product EJA, the
570 Jordan spin EJA, and the real symmetric matrices. The others
571 will have to be overridden.
573 return (X
.conjugate_transpose()*Y
).trace()
579 Return the unit element of this algebra.
583 sage: from mjo.eja.eja_algebra import (RealCartesianProductEJA,
588 sage: J = RealCartesianProductEJA(5)
590 e0 + e1 + e2 + e3 + e4
594 The identity element acts like the identity::
596 sage: set_random_seed()
597 sage: J = random_eja()
598 sage: x = J.random_element()
599 sage: J.one()*x == x and x*J.one() == x
602 The matrix of the unit element's operator is the identity::
604 sage: set_random_seed()
605 sage: J = random_eja()
606 sage: actual = J.one().operator().matrix()
607 sage: expected = matrix.identity(J.base_ring(), J.dimension())
608 sage: actual == expected
612 # We can brute-force compute the matrices of the operators
613 # that correspond to the basis elements of this algebra.
614 # If some linear combination of those basis elements is the
615 # algebra identity, then the same linear combination of
616 # their matrices has to be the identity matrix.
618 # Of course, matrices aren't vectors in sage, so we have to
619 # appeal to the "long vectors" isometry.
620 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
622 # Now we use basis linear algebra to find the coefficients,
623 # of the matrices-as-vectors-linear-combination, which should
624 # work for the original algebra basis too.
625 A
= matrix
.column(self
.base_ring(), oper_vecs
)
627 # We used the isometry on the left-hand side already, but we
628 # still need to do it for the right-hand side. Recall that we
629 # wanted something that summed to the identity matrix.
630 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
632 # Now if there's an identity element in the algebra, this should work.
633 coeffs
= A
.solve_right(b
)
634 return self
.linear_combination(zip(self
.gens(), coeffs
))
637 def peirce_decomposition(self
, c
):
639 The Peirce decomposition of this algebra relative to the
642 In the future, this can be extended to a complete system of
643 orthogonal idempotents.
645 if not c
.is_idempotent():
646 raise ValueError("element is not idempotent: %s" % c
)
648 # Default these to what they should be if they turn out to be
649 # trivial, because eigenspaces_left() won't return eigenvalues
650 # corresponding to trivial spaces (e.g. it returns only the
651 # eigenspace corresponding to lambda=1 if you take the
652 # decomposition relative to the identity element).
653 trivial
= FiniteDimensionalEuclideanJordanSubalgebra(self
, ())
654 J0
= trivial
# eigenvalue zero
655 J2
= trivial
# eigenvalue one-half
656 J1
= trivial
# eigenvalue one
658 for (eigval
, eigspace
) in c
.operator().matrix().left_eigenspaces():
659 gens
= tuple( self
.from_vector(b
) for b
in eigspace
.basis() )
660 subalg
= FiniteDimensionalEuclideanJordanSubalgebra(self
, gens
)
663 elif eigval
== ~
(self
.base_ring()(2)):
668 raise ValueError("unexpected eigenvalue: %s" % eigval
)
673 def random_elements(self
, count
):
675 Return ``count`` random elements as a tuple.
679 sage: from mjo.eja.eja_algebra import JordanSpinEJA
683 sage: J = JordanSpinEJA(3)
684 sage: x,y,z = J.random_elements(3)
685 sage: all( [ x in J, y in J, z in J ])
687 sage: len( J.random_elements(10) ) == 10
691 return tuple( self
.random_element() for idx
in xrange(count
) )
696 Return the rank of this EJA.
700 The author knows of no algorithm to compute the rank of an EJA
701 where only the multiplication table is known. In lieu of one, we
702 require the rank to be specified when the algebra is created,
703 and simply pass along that number here.
707 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
708 ....: RealSymmetricEJA,
709 ....: ComplexHermitianEJA,
710 ....: QuaternionHermitianEJA,
715 The rank of the Jordan spin algebra is always two::
717 sage: JordanSpinEJA(2).rank()
719 sage: JordanSpinEJA(3).rank()
721 sage: JordanSpinEJA(4).rank()
724 The rank of the `n`-by-`n` Hermitian real, complex, or
725 quaternion matrices is `n`::
727 sage: RealSymmetricEJA(4).rank()
729 sage: ComplexHermitianEJA(3).rank()
731 sage: QuaternionHermitianEJA(2).rank()
736 Ensure that every EJA that we know how to construct has a
737 positive integer rank, unless the algebra is trivial in
738 which case its rank will be zero::
740 sage: set_random_seed()
741 sage: J = random_eja()
745 sage: r > 0 or (r == 0 and J.is_trivial())
752 def vector_space(self
):
754 Return the vector space that underlies this algebra.
758 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
762 sage: J = RealSymmetricEJA(2)
763 sage: J.vector_space()
764 Vector space of dimension 3 over...
767 return self
.zero().to_vector().parent().ambient_vector_space()
770 Element
= FiniteDimensionalEuclideanJordanAlgebraElement
773 class KnownRankEJA(object):
775 A class for algebras that we actually know we can construct. The
776 main issue is that, for most of our methods to make sense, we need
777 to know the rank of our algebra. Thus we can't simply generate a
778 "random" algebra, or even check that a given basis and product
779 satisfy the axioms; because even if everything looks OK, we wouldn't
780 know the rank we need to actuallty build the thing.
782 Not really a subclass of FDEJA because doing that causes method
783 resolution errors, e.g.
785 TypeError: Error when calling the metaclass bases
786 Cannot create a consistent method resolution
787 order (MRO) for bases FiniteDimensionalEuclideanJordanAlgebra,
792 def _max_test_case_size():
794 Return an integer "size" that is an upper bound on the size of
795 this algebra when it is used in a random test
796 case. Unfortunately, the term "size" is quite vague -- when
797 dealing with `R^n` under either the Hadamard or Jordan spin
798 product, the "size" refers to the dimension `n`. When dealing
799 with a matrix algebra (real symmetric or complex/quaternion
800 Hermitian), it refers to the size of the matrix, which is
801 far less than the dimension of the underlying vector space.
803 We default to five in this class, which is safe in `R^n`. The
804 matrix algebra subclasses (or any class where the "size" is
805 interpreted to be far less than the dimension) should override
806 with a smaller number.
811 def random_instance(cls
, field
=QQ
, **kwargs
):
813 Return a random instance of this type of algebra.
815 Beware, this will crash for "most instances" because the
816 constructor below looks wrong.
818 if cls
is TrivialEJA
:
819 # The TrivialEJA class doesn't take an "n" argument because
823 n
= ZZ
.random_element(cls
._max
_test
_case
_size
()) + 1
824 return cls(n
, field
, **kwargs
)
827 class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra
,
830 Return the Euclidean Jordan Algebra corresponding to the set
831 `R^n` under the Hadamard product.
833 Note: this is nothing more than the Cartesian product of ``n``
834 copies of the spin algebra. Once Cartesian product algebras
835 are implemented, this can go.
839 sage: from mjo.eja.eja_algebra import RealCartesianProductEJA
843 This multiplication table can be verified by hand::
845 sage: J = RealCartesianProductEJA(3)
846 sage: e0,e1,e2 = J.gens()
862 We can change the generator prefix::
864 sage: RealCartesianProductEJA(3, prefix='r').gens()
868 def __init__(self
, n
, field
=QQ
, **kwargs
):
869 V
= VectorSpace(field
, n
)
870 mult_table
= [ [ V
.gen(i
)*(i
== j
) for j
in xrange(n
) ]
873 fdeja
= super(RealCartesianProductEJA
, self
)
874 return fdeja
.__init
__(field
, mult_table
, rank
=n
, **kwargs
)
876 def inner_product(self
, x
, y
):
878 Faster to reimplement than to use natural representations.
882 sage: from mjo.eja.eja_algebra import RealCartesianProductEJA
886 Ensure that this is the usual inner product for the algebras
889 sage: set_random_seed()
890 sage: J = RealCartesianProductEJA.random_instance()
891 sage: x,y = J.random_elements(2)
892 sage: X = x.natural_representation()
893 sage: Y = y.natural_representation()
894 sage: x.inner_product(y) == J.natural_inner_product(X,Y)
898 return x
.to_vector().inner_product(y
.to_vector())
901 def random_eja(field
=QQ
, nontrivial
=False):
903 Return a "random" finite-dimensional Euclidean Jordan Algebra.
907 sage: from mjo.eja.eja_algebra import random_eja
912 Euclidean Jordan algebra of dimension...
915 eja_classes
= KnownRankEJA
.__subclasses
__()
917 eja_classes
.remove(TrivialEJA
)
918 classname
= choice(eja_classes
)
919 return classname
.random_instance(field
=field
)
926 class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra
):
928 def _max_test_case_size():
929 # Play it safe, since this will be squared and the underlying
930 # field can have dimension 4 (quaternions) too.
933 def __init__(self
, field
, basis
, rank
, normalize_basis
=True, **kwargs
):
935 Compared to the superclass constructor, we take a basis instead of
936 a multiplication table because the latter can be computed in terms
937 of the former when the product is known (like it is here).
939 # Used in this class's fast _charpoly_coeff() override.
940 self
._basis
_normalizers
= None
942 # We're going to loop through this a few times, so now's a good
943 # time to ensure that it isn't a generator expression.
946 if rank
> 1 and normalize_basis
:
947 # We'll need sqrt(2) to normalize the basis, and this
948 # winds up in the multiplication table, so the whole
949 # algebra needs to be over the field extension.
950 R
= PolynomialRing(field
, 'z')
953 if p
.is_irreducible():
954 field
= field
.extension(p
, 'sqrt2', embedding
=RLF(2).sqrt())
955 basis
= tuple( s
.change_ring(field
) for s
in basis
)
956 self
._basis
_normalizers
= tuple(
957 ~
(self
.natural_inner_product(s
,s
).sqrt()) for s
in basis
)
958 basis
= tuple(s
*c
for (s
,c
) in izip(basis
,self
._basis
_normalizers
))
960 Qs
= self
.multiplication_table_from_matrix_basis(basis
)
962 fdeja
= super(MatrixEuclideanJordanAlgebra
, self
)
963 return fdeja
.__init
__(field
,
971 def _charpoly_coeff(self
, i
):
973 Override the parent method with something that tries to compute
974 over a faster (non-extension) field.
976 if self
._basis
_normalizers
is None:
977 # We didn't normalize, so assume that the basis we started
978 # with had entries in a nice field.
979 return super(MatrixEuclideanJordanAlgebra
, self
)._charpoly
_coeff
(i
)
981 basis
= ( (b
/n
) for (b
,n
) in izip(self
.natural_basis(),
982 self
._basis
_normalizers
) )
984 # Do this over the rationals and convert back at the end.
985 J
= MatrixEuclideanJordanAlgebra(QQ
,
988 normalize_basis
=False)
989 (_
,x
,_
,_
) = J
._charpoly
_matrix
_system
()
990 p
= J
._charpoly
_coeff
(i
)
991 # p might be missing some vars, have to substitute "optionally"
992 pairs
= izip(x
.base_ring().gens(), self
._basis
_normalizers
)
993 substitutions
= { v: v*c for (v,c) in pairs }
994 result
= p
.subs(substitutions
)
996 # The result of "subs" can be either a coefficient-ring
997 # element or a polynomial. Gotta handle both cases.
999 return self
.base_ring()(result
)
1001 return result
.change_ring(self
.base_ring())
1005 def multiplication_table_from_matrix_basis(basis
):
1007 At least three of the five simple Euclidean Jordan algebras have the
1008 symmetric multiplication (A,B) |-> (AB + BA)/2, where the
1009 multiplication on the right is matrix multiplication. Given a basis
1010 for the underlying matrix space, this function returns a
1011 multiplication table (obtained by looping through the basis
1012 elements) for an algebra of those matrices.
1014 # In S^2, for example, we nominally have four coordinates even
1015 # though the space is of dimension three only. The vector space V
1016 # is supposed to hold the entire long vector, and the subspace W
1017 # of V will be spanned by the vectors that arise from symmetric
1018 # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
1019 field
= basis
[0].base_ring()
1020 dimension
= basis
[0].nrows()
1022 V
= VectorSpace(field
, dimension
**2)
1023 W
= V
.span_of_basis( _mat2vec(s
) for s
in basis
)
1025 mult_table
= [[W
.zero() for j
in xrange(n
)] for i
in xrange(n
)]
1028 mat_entry
= (basis
[i
]*basis
[j
] + basis
[j
]*basis
[i
])/2
1029 mult_table
[i
][j
] = W
.coordinate_vector(_mat2vec(mat_entry
))
1037 Embed the matrix ``M`` into a space of real matrices.
1039 The matrix ``M`` can have entries in any field at the moment:
1040 the real numbers, complex numbers, or quaternions. And although
1041 they are not a field, we can probably support octonions at some
1042 point, too. This function returns a real matrix that "acts like"
1043 the original with respect to matrix multiplication; i.e.
1045 real_embed(M*N) = real_embed(M)*real_embed(N)
1048 raise NotImplementedError
1052 def real_unembed(M
):
1054 The inverse of :meth:`real_embed`.
1056 raise NotImplementedError
1060 def natural_inner_product(cls
,X
,Y
):
1061 Xu
= cls
.real_unembed(X
)
1062 Yu
= cls
.real_unembed(Y
)
1063 tr
= (Xu
*Yu
).trace()
1066 # It's real already.
1069 # Otherwise, try the thing that works for complex numbers; and
1070 # if that doesn't work, the thing that works for quaternions.
1072 return tr
.vector()[0] # real part, imag part is index 1
1073 except AttributeError:
1074 # A quaternions doesn't have a vector() method, but does
1075 # have coefficient_tuple() method that returns the
1076 # coefficients of 1, i, j, and k -- in that order.
1077 return tr
.coefficient_tuple()[0]
1080 class RealMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1084 The identity function, for embedding real matrices into real
1090 def real_unembed(M
):
1092 The identity function, for unembedding real matrices from real
1098 class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra
, KnownRankEJA
):
1100 The rank-n simple EJA consisting of real symmetric n-by-n
1101 matrices, the usual symmetric Jordan product, and the trace inner
1102 product. It has dimension `(n^2 + n)/2` over the reals.
1106 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1110 sage: J = RealSymmetricEJA(2)
1111 sage: e0, e1, e2 = J.gens()
1119 In theory, our "field" can be any subfield of the reals::
1121 sage: RealSymmetricEJA(2, AA)
1122 Euclidean Jordan algebra of dimension 3 over Algebraic Real Field
1123 sage: RealSymmetricEJA(2, RR)
1124 Euclidean Jordan algebra of dimension 3 over Real Field with
1125 53 bits of precision
1129 The dimension of this algebra is `(n^2 + n) / 2`::
1131 sage: set_random_seed()
1132 sage: n_max = RealSymmetricEJA._max_test_case_size()
1133 sage: n = ZZ.random_element(1, n_max)
1134 sage: J = RealSymmetricEJA(n)
1135 sage: J.dimension() == (n^2 + n)/2
1138 The Jordan multiplication is what we think it is::
1140 sage: set_random_seed()
1141 sage: J = RealSymmetricEJA.random_instance()
1142 sage: x,y = J.random_elements(2)
1143 sage: actual = (x*y).natural_representation()
1144 sage: X = x.natural_representation()
1145 sage: Y = y.natural_representation()
1146 sage: expected = (X*Y + Y*X)/2
1147 sage: actual == expected
1149 sage: J(expected) == x*y
1152 We can change the generator prefix::
1154 sage: RealSymmetricEJA(3, prefix='q').gens()
1155 (q0, q1, q2, q3, q4, q5)
1157 Our natural basis is normalized with respect to the natural inner
1158 product unless we specify otherwise::
1160 sage: set_random_seed()
1161 sage: J = RealSymmetricEJA.random_instance()
1162 sage: all( b.norm() == 1 for b in J.gens() )
1165 Since our natural basis is normalized with respect to the natural
1166 inner product, and since we know that this algebra is an EJA, any
1167 left-multiplication operator's matrix will be symmetric because
1168 natural->EJA basis representation is an isometry and within the EJA
1169 the operator is self-adjoint by the Jordan axiom::
1171 sage: set_random_seed()
1172 sage: x = RealSymmetricEJA.random_instance().random_element()
1173 sage: x.operator().matrix().is_symmetric()
1178 def _denormalized_basis(cls
, n
, field
):
1180 Return a basis for the space of real symmetric n-by-n matrices.
1184 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1188 sage: set_random_seed()
1189 sage: n = ZZ.random_element(1,5)
1190 sage: B = RealSymmetricEJA._denormalized_basis(n,QQ)
1191 sage: all( M.is_symmetric() for M in B)
1195 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1199 for j
in xrange(i
+1):
1200 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
1204 Sij
= Eij
+ Eij
.transpose()
1210 def _max_test_case_size():
1211 return 4 # Dimension 10
1214 def __init__(self
, n
, field
=QQ
, **kwargs
):
1215 basis
= self
._denormalized
_basis
(n
, field
)
1216 super(RealSymmetricEJA
, self
).__init
__(field
, basis
, n
, **kwargs
)
1219 class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1223 Embed the n-by-n complex matrix ``M`` into the space of real
1224 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1225 bi` to the block matrix ``[[a,b],[-b,a]]``.
1229 sage: from mjo.eja.eja_algebra import \
1230 ....: ComplexMatrixEuclideanJordanAlgebra
1234 sage: F = QuadraticField(-1, 'i')
1235 sage: x1 = F(4 - 2*i)
1236 sage: x2 = F(1 + 2*i)
1239 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1240 sage: ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1249 Embedding is a homomorphism (isomorphism, in fact)::
1251 sage: set_random_seed()
1252 sage: n_max = ComplexMatrixEuclideanJordanAlgebra._max_test_case_size()
1253 sage: n = ZZ.random_element(n_max)
1254 sage: F = QuadraticField(-1, 'i')
1255 sage: X = random_matrix(F, n)
1256 sage: Y = random_matrix(F, n)
1257 sage: Xe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X)
1258 sage: Ye = ComplexMatrixEuclideanJordanAlgebra.real_embed(Y)
1259 sage: XYe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1266 raise ValueError("the matrix 'M' must be square")
1268 # We don't need any adjoined elements...
1269 field
= M
.base_ring().base_ring()
1273 a
= z
.list()[0] # real part, I guess
1274 b
= z
.list()[1] # imag part, I guess
1275 blocks
.append(matrix(field
, 2, [[a
,b
],[-b
,a
]]))
1277 return matrix
.block(field
, n
, blocks
)
1281 def real_unembed(M
):
1283 The inverse of _embed_complex_matrix().
1287 sage: from mjo.eja.eja_algebra import \
1288 ....: ComplexMatrixEuclideanJordanAlgebra
1292 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1293 ....: [-2, 1, -4, 3],
1294 ....: [ 9, 10, 11, 12],
1295 ....: [-10, 9, -12, 11] ])
1296 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(A)
1298 [ 10*i + 9 12*i + 11]
1302 Unembedding is the inverse of embedding::
1304 sage: set_random_seed()
1305 sage: F = QuadraticField(-1, 'i')
1306 sage: M = random_matrix(F, 3)
1307 sage: Me = ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1308 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1314 raise ValueError("the matrix 'M' must be square")
1315 if not n
.mod(2).is_zero():
1316 raise ValueError("the matrix 'M' must be a complex embedding")
1318 # If "M" was normalized, its base ring might have roots
1319 # adjoined and they can stick around after unembedding.
1320 field
= M
.base_ring()
1321 R
= PolynomialRing(field
, 'z')
1323 F
= field
.extension(z
**2 + 1, 'i', embedding
=CLF(-1).sqrt())
1326 # Go top-left to bottom-right (reading order), converting every
1327 # 2-by-2 block we see to a single complex element.
1329 for k
in xrange(n
/2):
1330 for j
in xrange(n
/2):
1331 submat
= M
[2*k
:2*k
+2,2*j
:2*j
+2]
1332 if submat
[0,0] != submat
[1,1]:
1333 raise ValueError('bad on-diagonal submatrix')
1334 if submat
[0,1] != -submat
[1,0]:
1335 raise ValueError('bad off-diagonal submatrix')
1336 z
= submat
[0,0] + submat
[0,1]*i
1339 return matrix(F
, n
/2, elements
)
1343 def natural_inner_product(cls
,X
,Y
):
1345 Compute a natural inner product in this algebra directly from
1350 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1354 This gives the same answer as the slow, default method implemented
1355 in :class:`MatrixEuclideanJordanAlgebra`::
1357 sage: set_random_seed()
1358 sage: J = ComplexHermitianEJA.random_instance()
1359 sage: x,y = J.random_elements(2)
1360 sage: Xe = x.natural_representation()
1361 sage: Ye = y.natural_representation()
1362 sage: X = ComplexHermitianEJA.real_unembed(Xe)
1363 sage: Y = ComplexHermitianEJA.real_unembed(Ye)
1364 sage: expected = (X*Y).trace().vector()[0]
1365 sage: actual = ComplexHermitianEJA.natural_inner_product(Xe,Ye)
1366 sage: actual == expected
1370 return RealMatrixEuclideanJordanAlgebra
.natural_inner_product(X
,Y
)/2
1373 class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra
, KnownRankEJA
):
1375 The rank-n simple EJA consisting of complex Hermitian n-by-n
1376 matrices over the real numbers, the usual symmetric Jordan product,
1377 and the real-part-of-trace inner product. It has dimension `n^2` over
1382 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1386 In theory, our "field" can be any subfield of the reals::
1388 sage: ComplexHermitianEJA(2, AA)
1389 Euclidean Jordan algebra of dimension 4 over Algebraic Real Field
1390 sage: ComplexHermitianEJA(2, RR)
1391 Euclidean Jordan algebra of dimension 4 over Real Field with
1392 53 bits of precision
1396 The dimension of this algebra is `n^2`::
1398 sage: set_random_seed()
1399 sage: n_max = ComplexHermitianEJA._max_test_case_size()
1400 sage: n = ZZ.random_element(1, n_max)
1401 sage: J = ComplexHermitianEJA(n)
1402 sage: J.dimension() == n^2
1405 The Jordan multiplication is what we think it is::
1407 sage: set_random_seed()
1408 sage: J = ComplexHermitianEJA.random_instance()
1409 sage: x,y = J.random_elements(2)
1410 sage: actual = (x*y).natural_representation()
1411 sage: X = x.natural_representation()
1412 sage: Y = y.natural_representation()
1413 sage: expected = (X*Y + Y*X)/2
1414 sage: actual == expected
1416 sage: J(expected) == x*y
1419 We can change the generator prefix::
1421 sage: ComplexHermitianEJA(2, prefix='z').gens()
1424 Our natural basis is normalized with respect to the natural inner
1425 product unless we specify otherwise::
1427 sage: set_random_seed()
1428 sage: J = ComplexHermitianEJA.random_instance()
1429 sage: all( b.norm() == 1 for b in J.gens() )
1432 Since our natural basis is normalized with respect to the natural
1433 inner product, and since we know that this algebra is an EJA, any
1434 left-multiplication operator's matrix will be symmetric because
1435 natural->EJA basis representation is an isometry and within the EJA
1436 the operator is self-adjoint by the Jordan axiom::
1438 sage: set_random_seed()
1439 sage: x = ComplexHermitianEJA.random_instance().random_element()
1440 sage: x.operator().matrix().is_symmetric()
1446 def _denormalized_basis(cls
, n
, field
):
1448 Returns a basis for the space of complex Hermitian n-by-n matrices.
1450 Why do we embed these? Basically, because all of numerical linear
1451 algebra assumes that you're working with vectors consisting of `n`
1452 entries from a field and scalars from the same field. There's no way
1453 to tell SageMath that (for example) the vectors contain complex
1454 numbers, while the scalar field is real.
1458 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1462 sage: set_random_seed()
1463 sage: n = ZZ.random_element(1,5)
1464 sage: field = QuadraticField(2, 'sqrt2')
1465 sage: B = ComplexHermitianEJA._denormalized_basis(n, field)
1466 sage: all( M.is_symmetric() for M in B)
1470 R
= PolynomialRing(field
, 'z')
1472 F
= field
.extension(z
**2 + 1, 'I')
1475 # This is like the symmetric case, but we need to be careful:
1477 # * We want conjugate-symmetry, not just symmetry.
1478 # * The diagonal will (as a result) be real.
1482 for j
in xrange(i
+1):
1483 Eij
= matrix(F
, n
, lambda k
,l
: k
==i
and l
==j
)
1485 Sij
= cls
.real_embed(Eij
)
1488 # The second one has a minus because it's conjugated.
1489 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
1491 Sij_imag
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
1494 # Since we embedded these, we can drop back to the "field" that we
1495 # started with instead of the complex extension "F".
1496 return ( s
.change_ring(field
) for s
in S
)
1499 def __init__(self
, n
, field
=QQ
, **kwargs
):
1500 basis
= self
._denormalized
_basis
(n
,field
)
1501 super(ComplexHermitianEJA
,self
).__init
__(field
, basis
, n
, **kwargs
)
1504 class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1508 Embed the n-by-n quaternion matrix ``M`` into the space of real
1509 matrices of size 4n-by-4n by first sending each quaternion entry `z
1510 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
1511 c+di],[-c + di, a-bi]]`, and then embedding those into a real
1516 sage: from mjo.eja.eja_algebra import \
1517 ....: QuaternionMatrixEuclideanJordanAlgebra
1521 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1522 sage: i,j,k = Q.gens()
1523 sage: x = 1 + 2*i + 3*j + 4*k
1524 sage: M = matrix(Q, 1, [[x]])
1525 sage: QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1531 Embedding is a homomorphism (isomorphism, in fact)::
1533 sage: set_random_seed()
1534 sage: n_max = QuaternionMatrixEuclideanJordanAlgebra._max_test_case_size()
1535 sage: n = ZZ.random_element(n_max)
1536 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1537 sage: X = random_matrix(Q, n)
1538 sage: Y = random_matrix(Q, n)
1539 sage: Xe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X)
1540 sage: Ye = QuaternionMatrixEuclideanJordanAlgebra.real_embed(Y)
1541 sage: XYe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1546 quaternions
= M
.base_ring()
1549 raise ValueError("the matrix 'M' must be square")
1551 F
= QuadraticField(-1, 'i')
1556 t
= z
.coefficient_tuple()
1561 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
1562 [-c
+ d
*i
, a
- b
*i
]])
1563 realM
= ComplexMatrixEuclideanJordanAlgebra
.real_embed(cplxM
)
1564 blocks
.append(realM
)
1566 # We should have real entries by now, so use the realest field
1567 # we've got for the return value.
1568 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
1573 def real_unembed(M
):
1575 The inverse of _embed_quaternion_matrix().
1579 sage: from mjo.eja.eja_algebra import \
1580 ....: QuaternionMatrixEuclideanJordanAlgebra
1584 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
1585 ....: [-2, 1, -4, 3],
1586 ....: [-3, 4, 1, -2],
1587 ....: [-4, -3, 2, 1]])
1588 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(M)
1589 [1 + 2*i + 3*j + 4*k]
1593 Unembedding is the inverse of embedding::
1595 sage: set_random_seed()
1596 sage: Q = QuaternionAlgebra(QQ, -1, -1)
1597 sage: M = random_matrix(Q, 3)
1598 sage: Me = QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1599 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1605 raise ValueError("the matrix 'M' must be square")
1606 if not n
.mod(4).is_zero():
1607 raise ValueError("the matrix 'M' must be a quaternion embedding")
1609 # Use the base ring of the matrix to ensure that its entries can be
1610 # multiplied by elements of the quaternion algebra.
1611 field
= M
.base_ring()
1612 Q
= QuaternionAlgebra(field
,-1,-1)
1615 # Go top-left to bottom-right (reading order), converting every
1616 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
1619 for l
in xrange(n
/4):
1620 for m
in xrange(n
/4):
1621 submat
= ComplexMatrixEuclideanJordanAlgebra
.real_unembed(
1622 M
[4*l
:4*l
+4,4*m
:4*m
+4] )
1623 if submat
[0,0] != submat
[1,1].conjugate():
1624 raise ValueError('bad on-diagonal submatrix')
1625 if submat
[0,1] != -submat
[1,0].conjugate():
1626 raise ValueError('bad off-diagonal submatrix')
1627 z
= submat
[0,0].vector()[0] # real part
1628 z
+= submat
[0,0].vector()[1]*i
# imag part
1629 z
+= submat
[0,1].vector()[0]*j
# real part
1630 z
+= submat
[0,1].vector()[1]*k
# imag part
1633 return matrix(Q
, n
/4, elements
)
1637 def natural_inner_product(cls
,X
,Y
):
1639 Compute a natural inner product in this algebra directly from
1644 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1648 This gives the same answer as the slow, default method implemented
1649 in :class:`MatrixEuclideanJordanAlgebra`::
1651 sage: set_random_seed()
1652 sage: J = QuaternionHermitianEJA.random_instance()
1653 sage: x,y = J.random_elements(2)
1654 sage: Xe = x.natural_representation()
1655 sage: Ye = y.natural_representation()
1656 sage: X = QuaternionHermitianEJA.real_unembed(Xe)
1657 sage: Y = QuaternionHermitianEJA.real_unembed(Ye)
1658 sage: expected = (X*Y).trace().coefficient_tuple()[0]
1659 sage: actual = QuaternionHermitianEJA.natural_inner_product(Xe,Ye)
1660 sage: actual == expected
1664 return RealMatrixEuclideanJordanAlgebra
.natural_inner_product(X
,Y
)/4
1667 class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra
,
1670 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
1671 matrices, the usual symmetric Jordan product, and the
1672 real-part-of-trace inner product. It has dimension `2n^2 - n` over
1677 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1681 In theory, our "field" can be any subfield of the reals::
1683 sage: QuaternionHermitianEJA(2, AA)
1684 Euclidean Jordan algebra of dimension 6 over Algebraic Real Field
1685 sage: QuaternionHermitianEJA(2, RR)
1686 Euclidean Jordan algebra of dimension 6 over Real Field with
1687 53 bits of precision
1691 The dimension of this algebra is `2*n^2 - n`::
1693 sage: set_random_seed()
1694 sage: n_max = QuaternionHermitianEJA._max_test_case_size()
1695 sage: n = ZZ.random_element(1, n_max)
1696 sage: J = QuaternionHermitianEJA(n)
1697 sage: J.dimension() == 2*(n^2) - n
1700 The Jordan multiplication is what we think it is::
1702 sage: set_random_seed()
1703 sage: J = QuaternionHermitianEJA.random_instance()
1704 sage: x,y = J.random_elements(2)
1705 sage: actual = (x*y).natural_representation()
1706 sage: X = x.natural_representation()
1707 sage: Y = y.natural_representation()
1708 sage: expected = (X*Y + Y*X)/2
1709 sage: actual == expected
1711 sage: J(expected) == x*y
1714 We can change the generator prefix::
1716 sage: QuaternionHermitianEJA(2, prefix='a').gens()
1717 (a0, a1, a2, a3, a4, a5)
1719 Our natural basis is normalized with respect to the natural inner
1720 product unless we specify otherwise::
1722 sage: set_random_seed()
1723 sage: J = QuaternionHermitianEJA.random_instance()
1724 sage: all( b.norm() == 1 for b in J.gens() )
1727 Since our natural basis is normalized with respect to the natural
1728 inner product, and since we know that this algebra is an EJA, any
1729 left-multiplication operator's matrix will be symmetric because
1730 natural->EJA basis representation is an isometry and within the EJA
1731 the operator is self-adjoint by the Jordan axiom::
1733 sage: set_random_seed()
1734 sage: x = QuaternionHermitianEJA.random_instance().random_element()
1735 sage: x.operator().matrix().is_symmetric()
1740 def _denormalized_basis(cls
, n
, field
):
1742 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
1744 Why do we embed these? Basically, because all of numerical
1745 linear algebra assumes that you're working with vectors consisting
1746 of `n` entries from a field and scalars from the same field. There's
1747 no way to tell SageMath that (for example) the vectors contain
1748 complex numbers, while the scalar field is real.
1752 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1756 sage: set_random_seed()
1757 sage: n = ZZ.random_element(1,5)
1758 sage: B = QuaternionHermitianEJA._denormalized_basis(n,QQ)
1759 sage: all( M.is_symmetric() for M in B )
1763 Q
= QuaternionAlgebra(QQ
,-1,-1)
1766 # This is like the symmetric case, but we need to be careful:
1768 # * We want conjugate-symmetry, not just symmetry.
1769 # * The diagonal will (as a result) be real.
1773 for j
in xrange(i
+1):
1774 Eij
= matrix(Q
, n
, lambda k
,l
: k
==i
and l
==j
)
1776 Sij
= cls
.real_embed(Eij
)
1779 # The second, third, and fourth ones have a minus
1780 # because they're conjugated.
1781 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
1783 Sij_I
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
1785 Sij_J
= cls
.real_embed(J
*Eij
- J
*Eij
.transpose())
1787 Sij_K
= cls
.real_embed(K
*Eij
- K
*Eij
.transpose())
1790 # Since we embedded these, we can drop back to the "field" that we
1791 # started with instead of the quaternion algebra "Q".
1792 return ( s
.change_ring(field
) for s
in S
)
1795 def __init__(self
, n
, field
=QQ
, **kwargs
):
1796 basis
= self
._denormalized
_basis
(n
,field
)
1797 super(QuaternionHermitianEJA
,self
).__init
__(field
, basis
, n
, **kwargs
)
1800 class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra
, KnownRankEJA
):
1802 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
1803 with the usual inner product and jordan product ``x*y =
1804 (<x_bar,y_bar>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
1809 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1813 This multiplication table can be verified by hand::
1815 sage: J = JordanSpinEJA(4)
1816 sage: e0,e1,e2,e3 = J.gens()
1832 We can change the generator prefix::
1834 sage: JordanSpinEJA(2, prefix='B').gens()
1838 def __init__(self
, n
, field
=QQ
, **kwargs
):
1839 V
= VectorSpace(field
, n
)
1840 mult_table
= [[V
.zero() for j
in xrange(n
)] for i
in xrange(n
)]
1850 z0
= x
.inner_product(y
)
1851 zbar
= y0
*xbar
+ x0
*ybar
1852 z
= V([z0
] + zbar
.list())
1853 mult_table
[i
][j
] = z
1855 # The rank of the spin algebra is two, unless we're in a
1856 # one-dimensional ambient space (because the rank is bounded by
1857 # the ambient dimension).
1858 fdeja
= super(JordanSpinEJA
, self
)
1859 return fdeja
.__init
__(field
, mult_table
, rank
=min(n
,2), **kwargs
)
1861 def inner_product(self
, x
, y
):
1863 Faster to reimplement than to use natural representations.
1867 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1871 Ensure that this is the usual inner product for the algebras
1874 sage: set_random_seed()
1875 sage: J = JordanSpinEJA.random_instance()
1876 sage: x,y = J.random_elements(2)
1877 sage: X = x.natural_representation()
1878 sage: Y = y.natural_representation()
1879 sage: x.inner_product(y) == J.natural_inner_product(X,Y)
1883 return x
.to_vector().inner_product(y
.to_vector())
1886 class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra
, KnownRankEJA
):
1888 The trivial Euclidean Jordan algebra consisting of only a zero element.
1892 sage: from mjo.eja.eja_algebra import TrivialEJA
1896 sage: J = TrivialEJA()
1903 sage: 7*J.one()*12*J.one()
1905 sage: J.one().inner_product(J.one())
1907 sage: J.one().norm()
1909 sage: J.one().subalgebra_generated_by()
1910 Euclidean Jordan algebra of dimension 0 over Rational Field
1915 def __init__(self
, field
=QQ
, **kwargs
):
1917 fdeja
= super(TrivialEJA
, self
)
1918 # The rank is zero using my definition, namely the dimension of the
1919 # largest subalgebra generated by any element.
1920 return fdeja
.__init
__(field
, mult_table
, rank
=0, **kwargs
)