2 Euclidean Jordan Algebras. These are formally-real Jordan Algebras;
3 specifically those where u^2 + v^2 = 0 implies that u = v = 0. They
4 are used in optimization, and have some additional nice methods beyond
5 what can be supported in a general Jordan Algebra.
8 from itertools
import repeat
10 from sage
.algebras
.quatalg
.quaternion_algebra
import QuaternionAlgebra
11 from sage
.categories
.magmatic_algebras
import MagmaticAlgebras
12 from sage
.combinat
.free_module
import CombinatorialFreeModule
13 from sage
.matrix
.constructor
import matrix
14 from sage
.matrix
.matrix_space
import MatrixSpace
15 from sage
.misc
.cachefunc
import cached_method
16 from sage
.misc
.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 range(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 range(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.
647 - ``c`` -- an idempotent of this algebra.
651 A triple (J0, J5, J1) containing two subalgebras and one subspace
654 - ``J0`` -- the algebra on the eigenspace of ``c.operator()``
655 corresponding to the eigenvalue zero.
657 - ``J5`` -- the eigenspace (NOT a subalgebra) of ``c.operator()``
658 corresponding to the eigenvalue one-half.
660 - ``J1`` -- the algebra on the eigenspace of ``c.operator()``
661 corresponding to the eigenvalue one.
663 These are the only possible eigenspaces for that operator, and this
664 algebra is a direct sum of them. The spaces ``J0`` and ``J1`` are
665 orthogonal, and are subalgebras of this algebra with the appropriate
670 sage: from mjo.eja.eja_algebra import random_eja, RealSymmetricEJA
674 The canonical example comes from the symmetric matrices, which
675 decompose into diagonal and off-diagonal parts::
677 sage: J = RealSymmetricEJA(3)
678 sage: C = matrix(QQ, [ [1,0,0],
682 sage: J0,J5,J1 = J.peirce_decomposition(c)
684 Euclidean Jordan algebra of dimension 1...
686 Vector space of degree 6 and dimension 2...
688 Euclidean Jordan algebra of dimension 3...
692 Every algebra decomposes trivially with respect to its identity
695 sage: set_random_seed()
696 sage: J = random_eja()
697 sage: J0,J5,J1 = J.peirce_decomposition(J.one())
698 sage: J0.dimension() == 0 and J5.dimension() == 0
700 sage: J1.superalgebra() == J and J1.dimension() == J.dimension()
703 The identity elements in the two subalgebras are the
704 projections onto their respective subspaces of the
705 superalgebra's identity element::
707 sage: set_random_seed()
708 sage: J = random_eja()
709 sage: x = J.random_element()
710 sage: if not J.is_trivial():
711 ....: while x.is_nilpotent():
712 ....: x = J.random_element()
713 sage: c = x.subalgebra_idempotent()
714 sage: J0,J5,J1 = J.peirce_decomposition(c)
715 sage: J1(c) == J1.one()
717 sage: J0(J.one() - c) == J0.one()
721 if not c
.is_idempotent():
722 raise ValueError("element is not idempotent: %s" % c
)
724 # Default these to what they should be if they turn out to be
725 # trivial, because eigenspaces_left() won't return eigenvalues
726 # corresponding to trivial spaces (e.g. it returns only the
727 # eigenspace corresponding to lambda=1 if you take the
728 # decomposition relative to the identity element).
729 trivial
= FiniteDimensionalEuclideanJordanSubalgebra(self
, ())
730 J0
= trivial
# eigenvalue zero
731 J5
= VectorSpace(self
.base_ring(), 0) # eigenvalue one-half
732 J1
= trivial
# eigenvalue one
734 for (eigval
, eigspace
) in c
.operator().matrix().left_eigenspaces():
735 if eigval
== ~
(self
.base_ring()(2)):
738 gens
= tuple( self
.from_vector(b
) for b
in eigspace
.basis() )
739 subalg
= FiniteDimensionalEuclideanJordanSubalgebra(self
, gens
)
745 raise ValueError("unexpected eigenvalue: %s" % eigval
)
750 def random_elements(self
, count
):
752 Return ``count`` random elements as a tuple.
756 sage: from mjo.eja.eja_algebra import JordanSpinEJA
760 sage: J = JordanSpinEJA(3)
761 sage: x,y,z = J.random_elements(3)
762 sage: all( [ x in J, y in J, z in J ])
764 sage: len( J.random_elements(10) ) == 10
768 return tuple( self
.random_element() for idx
in range(count
) )
773 Return the rank of this EJA.
777 The author knows of no algorithm to compute the rank of an EJA
778 where only the multiplication table is known. In lieu of one, we
779 require the rank to be specified when the algebra is created,
780 and simply pass along that number here.
784 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
785 ....: RealSymmetricEJA,
786 ....: ComplexHermitianEJA,
787 ....: QuaternionHermitianEJA,
792 The rank of the Jordan spin algebra is always two::
794 sage: JordanSpinEJA(2).rank()
796 sage: JordanSpinEJA(3).rank()
798 sage: JordanSpinEJA(4).rank()
801 The rank of the `n`-by-`n` Hermitian real, complex, or
802 quaternion matrices is `n`::
804 sage: RealSymmetricEJA(4).rank()
806 sage: ComplexHermitianEJA(3).rank()
808 sage: QuaternionHermitianEJA(2).rank()
813 Ensure that every EJA that we know how to construct has a
814 positive integer rank, unless the algebra is trivial in
815 which case its rank will be zero::
817 sage: set_random_seed()
818 sage: J = random_eja()
822 sage: r > 0 or (r == 0 and J.is_trivial())
829 def vector_space(self
):
831 Return the vector space that underlies this algebra.
835 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
839 sage: J = RealSymmetricEJA(2)
840 sage: J.vector_space()
841 Vector space of dimension 3 over...
844 return self
.zero().to_vector().parent().ambient_vector_space()
847 Element
= FiniteDimensionalEuclideanJordanAlgebraElement
850 class KnownRankEJA(object):
852 A class for algebras that we actually know we can construct. The
853 main issue is that, for most of our methods to make sense, we need
854 to know the rank of our algebra. Thus we can't simply generate a
855 "random" algebra, or even check that a given basis and product
856 satisfy the axioms; because even if everything looks OK, we wouldn't
857 know the rank we need to actuallty build the thing.
859 Not really a subclass of FDEJA because doing that causes method
860 resolution errors, e.g.
862 TypeError: Error when calling the metaclass bases
863 Cannot create a consistent method resolution
864 order (MRO) for bases FiniteDimensionalEuclideanJordanAlgebra,
869 def _max_test_case_size():
871 Return an integer "size" that is an upper bound on the size of
872 this algebra when it is used in a random test
873 case. Unfortunately, the term "size" is quite vague -- when
874 dealing with `R^n` under either the Hadamard or Jordan spin
875 product, the "size" refers to the dimension `n`. When dealing
876 with a matrix algebra (real symmetric or complex/quaternion
877 Hermitian), it refers to the size of the matrix, which is
878 far less than the dimension of the underlying vector space.
880 We default to five in this class, which is safe in `R^n`. The
881 matrix algebra subclasses (or any class where the "size" is
882 interpreted to be far less than the dimension) should override
883 with a smaller number.
888 def random_instance(cls
, field
=QQ
, **kwargs
):
890 Return a random instance of this type of algebra.
892 Beware, this will crash for "most instances" because the
893 constructor below looks wrong.
895 if cls
is TrivialEJA
:
896 # The TrivialEJA class doesn't take an "n" argument because
900 n
= ZZ
.random_element(cls
._max
_test
_case
_size
()) + 1
901 return cls(n
, field
, **kwargs
)
904 class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra
,
907 Return the Euclidean Jordan Algebra corresponding to the set
908 `R^n` under the Hadamard product.
910 Note: this is nothing more than the Cartesian product of ``n``
911 copies of the spin algebra. Once Cartesian product algebras
912 are implemented, this can go.
916 sage: from mjo.eja.eja_algebra import RealCartesianProductEJA
920 This multiplication table can be verified by hand::
922 sage: J = RealCartesianProductEJA(3)
923 sage: e0,e1,e2 = J.gens()
939 We can change the generator prefix::
941 sage: RealCartesianProductEJA(3, prefix='r').gens()
945 def __init__(self
, n
, field
=QQ
, **kwargs
):
946 V
= VectorSpace(field
, n
)
947 mult_table
= [ [ V
.gen(i
)*(i
== j
) for j
in range(n
) ]
950 fdeja
= super(RealCartesianProductEJA
, self
)
951 return fdeja
.__init
__(field
, mult_table
, rank
=n
, **kwargs
)
953 def inner_product(self
, x
, y
):
955 Faster to reimplement than to use natural representations.
959 sage: from mjo.eja.eja_algebra import RealCartesianProductEJA
963 Ensure that this is the usual inner product for the algebras
966 sage: set_random_seed()
967 sage: J = RealCartesianProductEJA.random_instance()
968 sage: x,y = J.random_elements(2)
969 sage: X = x.natural_representation()
970 sage: Y = y.natural_representation()
971 sage: x.inner_product(y) == J.natural_inner_product(X,Y)
975 return x
.to_vector().inner_product(y
.to_vector())
978 def random_eja(field
=QQ
, nontrivial
=False):
980 Return a "random" finite-dimensional Euclidean Jordan Algebra.
984 sage: from mjo.eja.eja_algebra import random_eja
989 Euclidean Jordan algebra of dimension...
992 eja_classes
= KnownRankEJA
.__subclasses
__()
994 eja_classes
.remove(TrivialEJA
)
995 classname
= choice(eja_classes
)
996 return classname
.random_instance(field
=field
)
1003 class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra
):
1005 def _max_test_case_size():
1006 # Play it safe, since this will be squared and the underlying
1007 # field can have dimension 4 (quaternions) too.
1010 def __init__(self
, field
, basis
, rank
, normalize_basis
=True, **kwargs
):
1012 Compared to the superclass constructor, we take a basis instead of
1013 a multiplication table because the latter can be computed in terms
1014 of the former when the product is known (like it is here).
1016 # Used in this class's fast _charpoly_coeff() override.
1017 self
._basis
_normalizers
= None
1019 # We're going to loop through this a few times, so now's a good
1020 # time to ensure that it isn't a generator expression.
1021 basis
= tuple(basis
)
1023 if rank
> 1 and normalize_basis
:
1024 # We'll need sqrt(2) to normalize the basis, and this
1025 # winds up in the multiplication table, so the whole
1026 # algebra needs to be over the field extension.
1027 R
= PolynomialRing(field
, 'z')
1030 if p
.is_irreducible():
1031 field
= field
.extension(p
, 'sqrt2', embedding
=RLF(2).sqrt())
1032 basis
= tuple( s
.change_ring(field
) for s
in basis
)
1033 self
._basis
_normalizers
= tuple(
1034 ~
(self
.natural_inner_product(s
,s
).sqrt()) for s
in basis
)
1035 basis
= tuple(s
*c
for (s
,c
) in zip(basis
,self
._basis
_normalizers
))
1037 Qs
= self
.multiplication_table_from_matrix_basis(basis
)
1039 fdeja
= super(MatrixEuclideanJordanAlgebra
, self
)
1040 return fdeja
.__init
__(field
,
1043 natural_basis
=basis
,
1048 def _charpoly_coeff(self
, i
):
1050 Override the parent method with something that tries to compute
1051 over a faster (non-extension) field.
1053 if self
._basis
_normalizers
is None:
1054 # We didn't normalize, so assume that the basis we started
1055 # with had entries in a nice field.
1056 return super(MatrixEuclideanJordanAlgebra
, self
)._charpoly
_coeff
(i
)
1058 basis
= ( (b
/n
) for (b
,n
) in zip(self
.natural_basis(),
1059 self
._basis
_normalizers
) )
1061 # Do this over the rationals and convert back at the end.
1062 J
= MatrixEuclideanJordanAlgebra(QQ
,
1065 normalize_basis
=False)
1066 (_
,x
,_
,_
) = J
._charpoly
_matrix
_system
()
1067 p
= J
._charpoly
_coeff
(i
)
1068 # p might be missing some vars, have to substitute "optionally"
1069 pairs
= zip(x
.base_ring().gens(), self
._basis
_normalizers
)
1070 substitutions
= { v: v*c for (v,c) in pairs }
1071 result
= p
.subs(substitutions
)
1073 # The result of "subs" can be either a coefficient-ring
1074 # element or a polynomial. Gotta handle both cases.
1076 return self
.base_ring()(result
)
1078 return result
.change_ring(self
.base_ring())
1082 def multiplication_table_from_matrix_basis(basis
):
1084 At least three of the five simple Euclidean Jordan algebras have the
1085 symmetric multiplication (A,B) |-> (AB + BA)/2, where the
1086 multiplication on the right is matrix multiplication. Given a basis
1087 for the underlying matrix space, this function returns a
1088 multiplication table (obtained by looping through the basis
1089 elements) for an algebra of those matrices.
1091 # In S^2, for example, we nominally have four coordinates even
1092 # though the space is of dimension three only. The vector space V
1093 # is supposed to hold the entire long vector, and the subspace W
1094 # of V will be spanned by the vectors that arise from symmetric
1095 # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
1096 field
= basis
[0].base_ring()
1097 dimension
= basis
[0].nrows()
1099 V
= VectorSpace(field
, dimension
**2)
1100 W
= V
.span_of_basis( _mat2vec(s
) for s
in basis
)
1102 mult_table
= [[W
.zero() for j
in range(n
)] for i
in range(n
)]
1105 mat_entry
= (basis
[i
]*basis
[j
] + basis
[j
]*basis
[i
])/2
1106 mult_table
[i
][j
] = W
.coordinate_vector(_mat2vec(mat_entry
))
1114 Embed the matrix ``M`` into a space of real matrices.
1116 The matrix ``M`` can have entries in any field at the moment:
1117 the real numbers, complex numbers, or quaternions. And although
1118 they are not a field, we can probably support octonions at some
1119 point, too. This function returns a real matrix that "acts like"
1120 the original with respect to matrix multiplication; i.e.
1122 real_embed(M*N) = real_embed(M)*real_embed(N)
1125 raise NotImplementedError
1129 def real_unembed(M
):
1131 The inverse of :meth:`real_embed`.
1133 raise NotImplementedError
1137 def natural_inner_product(cls
,X
,Y
):
1138 Xu
= cls
.real_unembed(X
)
1139 Yu
= cls
.real_unembed(Y
)
1140 tr
= (Xu
*Yu
).trace()
1143 # It's real already.
1146 # Otherwise, try the thing that works for complex numbers; and
1147 # if that doesn't work, the thing that works for quaternions.
1149 return tr
.vector()[0] # real part, imag part is index 1
1150 except AttributeError:
1151 # A quaternions doesn't have a vector() method, but does
1152 # have coefficient_tuple() method that returns the
1153 # coefficients of 1, i, j, and k -- in that order.
1154 return tr
.coefficient_tuple()[0]
1157 class RealMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1161 The identity function, for embedding real matrices into real
1167 def real_unembed(M
):
1169 The identity function, for unembedding real matrices from real
1175 class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra
, KnownRankEJA
):
1177 The rank-n simple EJA consisting of real symmetric n-by-n
1178 matrices, the usual symmetric Jordan product, and the trace inner
1179 product. It has dimension `(n^2 + n)/2` over the reals.
1183 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1187 sage: J = RealSymmetricEJA(2)
1188 sage: e0, e1, e2 = J.gens()
1196 In theory, our "field" can be any subfield of the reals::
1198 sage: RealSymmetricEJA(2, AA)
1199 Euclidean Jordan algebra of dimension 3 over Algebraic Real Field
1200 sage: RealSymmetricEJA(2, RR)
1201 Euclidean Jordan algebra of dimension 3 over Real Field with
1202 53 bits of precision
1206 The dimension of this algebra is `(n^2 + n) / 2`::
1208 sage: set_random_seed()
1209 sage: n_max = RealSymmetricEJA._max_test_case_size()
1210 sage: n = ZZ.random_element(1, n_max)
1211 sage: J = RealSymmetricEJA(n)
1212 sage: J.dimension() == (n^2 + n)/2
1215 The Jordan multiplication is what we think it is::
1217 sage: set_random_seed()
1218 sage: J = RealSymmetricEJA.random_instance()
1219 sage: x,y = J.random_elements(2)
1220 sage: actual = (x*y).natural_representation()
1221 sage: X = x.natural_representation()
1222 sage: Y = y.natural_representation()
1223 sage: expected = (X*Y + Y*X)/2
1224 sage: actual == expected
1226 sage: J(expected) == x*y
1229 We can change the generator prefix::
1231 sage: RealSymmetricEJA(3, prefix='q').gens()
1232 (q0, q1, q2, q3, q4, q5)
1234 Our natural basis is normalized with respect to the natural inner
1235 product unless we specify otherwise::
1237 sage: set_random_seed()
1238 sage: J = RealSymmetricEJA.random_instance()
1239 sage: all( b.norm() == 1 for b in J.gens() )
1242 Since our natural basis is normalized with respect to the natural
1243 inner product, and since we know that this algebra is an EJA, any
1244 left-multiplication operator's matrix will be symmetric because
1245 natural->EJA basis representation is an isometry and within the EJA
1246 the operator is self-adjoint by the Jordan axiom::
1248 sage: set_random_seed()
1249 sage: x = RealSymmetricEJA.random_instance().random_element()
1250 sage: x.operator().matrix().is_symmetric()
1255 def _denormalized_basis(cls
, n
, field
):
1257 Return a basis for the space of real symmetric n-by-n matrices.
1261 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1265 sage: set_random_seed()
1266 sage: n = ZZ.random_element(1,5)
1267 sage: B = RealSymmetricEJA._denormalized_basis(n,QQ)
1268 sage: all( M.is_symmetric() for M in B)
1272 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1276 for j
in range(i
+1):
1277 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
1281 Sij
= Eij
+ Eij
.transpose()
1287 def _max_test_case_size():
1288 return 4 # Dimension 10
1291 def __init__(self
, n
, field
=QQ
, **kwargs
):
1292 basis
= self
._denormalized
_basis
(n
, field
)
1293 super(RealSymmetricEJA
, self
).__init
__(field
, basis
, n
, **kwargs
)
1296 class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1300 Embed the n-by-n complex matrix ``M`` into the space of real
1301 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1302 bi` to the block matrix ``[[a,b],[-b,a]]``.
1306 sage: from mjo.eja.eja_algebra import \
1307 ....: ComplexMatrixEuclideanJordanAlgebra
1311 sage: F = QuadraticField(-1, 'i')
1312 sage: x1 = F(4 - 2*i)
1313 sage: x2 = F(1 + 2*i)
1316 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1317 sage: ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1326 Embedding is a homomorphism (isomorphism, in fact)::
1328 sage: set_random_seed()
1329 sage: n_max = ComplexMatrixEuclideanJordanAlgebra._max_test_case_size()
1330 sage: n = ZZ.random_element(n_max)
1331 sage: F = QuadraticField(-1, 'i')
1332 sage: X = random_matrix(F, n)
1333 sage: Y = random_matrix(F, n)
1334 sage: Xe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X)
1335 sage: Ye = ComplexMatrixEuclideanJordanAlgebra.real_embed(Y)
1336 sage: XYe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1343 raise ValueError("the matrix 'M' must be square")
1345 # We don't need any adjoined elements...
1346 field
= M
.base_ring().base_ring()
1350 a
= z
.list()[0] # real part, I guess
1351 b
= z
.list()[1] # imag part, I guess
1352 blocks
.append(matrix(field
, 2, [[a
,b
],[-b
,a
]]))
1354 return matrix
.block(field
, n
, blocks
)
1358 def real_unembed(M
):
1360 The inverse of _embed_complex_matrix().
1364 sage: from mjo.eja.eja_algebra import \
1365 ....: ComplexMatrixEuclideanJordanAlgebra
1369 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1370 ....: [-2, 1, -4, 3],
1371 ....: [ 9, 10, 11, 12],
1372 ....: [-10, 9, -12, 11] ])
1373 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(A)
1375 [ 10*i + 9 12*i + 11]
1379 Unembedding is the inverse of embedding::
1381 sage: set_random_seed()
1382 sage: F = QuadraticField(-1, 'i')
1383 sage: M = random_matrix(F, 3)
1384 sage: Me = ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1385 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1391 raise ValueError("the matrix 'M' must be square")
1392 if not n
.mod(2).is_zero():
1393 raise ValueError("the matrix 'M' must be a complex embedding")
1395 # If "M" was normalized, its base ring might have roots
1396 # adjoined and they can stick around after unembedding.
1397 field
= M
.base_ring()
1398 R
= PolynomialRing(field
, 'z')
1400 F
= field
.extension(z
**2 + 1, 'i', embedding
=CLF(-1).sqrt())
1403 # Go top-left to bottom-right (reading order), converting every
1404 # 2-by-2 block we see to a single complex element.
1406 for k
in range(n
/2):
1407 for j
in range(n
/2):
1408 submat
= M
[2*k
:2*k
+2,2*j
:2*j
+2]
1409 if submat
[0,0] != submat
[1,1]:
1410 raise ValueError('bad on-diagonal submatrix')
1411 if submat
[0,1] != -submat
[1,0]:
1412 raise ValueError('bad off-diagonal submatrix')
1413 z
= submat
[0,0] + submat
[0,1]*i
1416 return matrix(F
, n
/2, elements
)
1420 def natural_inner_product(cls
,X
,Y
):
1422 Compute a natural inner product in this algebra directly from
1427 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1431 This gives the same answer as the slow, default method implemented
1432 in :class:`MatrixEuclideanJordanAlgebra`::
1434 sage: set_random_seed()
1435 sage: J = ComplexHermitianEJA.random_instance()
1436 sage: x,y = J.random_elements(2)
1437 sage: Xe = x.natural_representation()
1438 sage: Ye = y.natural_representation()
1439 sage: X = ComplexHermitianEJA.real_unembed(Xe)
1440 sage: Y = ComplexHermitianEJA.real_unembed(Ye)
1441 sage: expected = (X*Y).trace().vector()[0]
1442 sage: actual = ComplexHermitianEJA.natural_inner_product(Xe,Ye)
1443 sage: actual == expected
1447 return RealMatrixEuclideanJordanAlgebra
.natural_inner_product(X
,Y
)/2
1450 class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra
, KnownRankEJA
):
1452 The rank-n simple EJA consisting of complex Hermitian n-by-n
1453 matrices over the real numbers, the usual symmetric Jordan product,
1454 and the real-part-of-trace inner product. It has dimension `n^2` over
1459 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1463 In theory, our "field" can be any subfield of the reals::
1465 sage: ComplexHermitianEJA(2, AA)
1466 Euclidean Jordan algebra of dimension 4 over Algebraic Real Field
1467 sage: ComplexHermitianEJA(2, RR)
1468 Euclidean Jordan algebra of dimension 4 over Real Field with
1469 53 bits of precision
1473 The dimension of this algebra is `n^2`::
1475 sage: set_random_seed()
1476 sage: n_max = ComplexHermitianEJA._max_test_case_size()
1477 sage: n = ZZ.random_element(1, n_max)
1478 sage: J = ComplexHermitianEJA(n)
1479 sage: J.dimension() == n^2
1482 The Jordan multiplication is what we think it is::
1484 sage: set_random_seed()
1485 sage: J = ComplexHermitianEJA.random_instance()
1486 sage: x,y = J.random_elements(2)
1487 sage: actual = (x*y).natural_representation()
1488 sage: X = x.natural_representation()
1489 sage: Y = y.natural_representation()
1490 sage: expected = (X*Y + Y*X)/2
1491 sage: actual == expected
1493 sage: J(expected) == x*y
1496 We can change the generator prefix::
1498 sage: ComplexHermitianEJA(2, prefix='z').gens()
1501 Our natural basis is normalized with respect to the natural inner
1502 product unless we specify otherwise::
1504 sage: set_random_seed()
1505 sage: J = ComplexHermitianEJA.random_instance()
1506 sage: all( b.norm() == 1 for b in J.gens() )
1509 Since our natural basis is normalized with respect to the natural
1510 inner product, and since we know that this algebra is an EJA, any
1511 left-multiplication operator's matrix will be symmetric because
1512 natural->EJA basis representation is an isometry and within the EJA
1513 the operator is self-adjoint by the Jordan axiom::
1515 sage: set_random_seed()
1516 sage: x = ComplexHermitianEJA.random_instance().random_element()
1517 sage: x.operator().matrix().is_symmetric()
1523 def _denormalized_basis(cls
, n
, field
):
1525 Returns a basis for the space of complex Hermitian n-by-n matrices.
1527 Why do we embed these? Basically, because all of numerical linear
1528 algebra assumes that you're working with vectors consisting of `n`
1529 entries from a field and scalars from the same field. There's no way
1530 to tell SageMath that (for example) the vectors contain complex
1531 numbers, while the scalar field is real.
1535 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1539 sage: set_random_seed()
1540 sage: n = ZZ.random_element(1,5)
1541 sage: field = QuadraticField(2, 'sqrt2')
1542 sage: B = ComplexHermitianEJA._denormalized_basis(n, field)
1543 sage: all( M.is_symmetric() for M in B)
1547 R
= PolynomialRing(field
, 'z')
1549 F
= field
.extension(z
**2 + 1, 'I')
1552 # This is like the symmetric case, but we need to be careful:
1554 # * We want conjugate-symmetry, not just symmetry.
1555 # * The diagonal will (as a result) be real.
1559 for j
in range(i
+1):
1560 Eij
= matrix(F
, n
, lambda k
,l
: k
==i
and l
==j
)
1562 Sij
= cls
.real_embed(Eij
)
1565 # The second one has a minus because it's conjugated.
1566 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
1568 Sij_imag
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
1571 # Since we embedded these, we can drop back to the "field" that we
1572 # started with instead of the complex extension "F".
1573 return ( s
.change_ring(field
) for s
in S
)
1576 def __init__(self
, n
, field
=QQ
, **kwargs
):
1577 basis
= self
._denormalized
_basis
(n
,field
)
1578 super(ComplexHermitianEJA
,self
).__init
__(field
, basis
, n
, **kwargs
)
1581 class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1585 Embed the n-by-n quaternion matrix ``M`` into the space of real
1586 matrices of size 4n-by-4n by first sending each quaternion entry `z
1587 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
1588 c+di],[-c + di, a-bi]]`, and then embedding those into a real
1593 sage: from mjo.eja.eja_algebra import \
1594 ....: QuaternionMatrixEuclideanJordanAlgebra
1598 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1599 sage: i,j,k = Q.gens()
1600 sage: x = 1 + 2*i + 3*j + 4*k
1601 sage: M = matrix(Q, 1, [[x]])
1602 sage: QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1608 Embedding is a homomorphism (isomorphism, in fact)::
1610 sage: set_random_seed()
1611 sage: n_max = QuaternionMatrixEuclideanJordanAlgebra._max_test_case_size()
1612 sage: n = ZZ.random_element(n_max)
1613 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1614 sage: X = random_matrix(Q, n)
1615 sage: Y = random_matrix(Q, n)
1616 sage: Xe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X)
1617 sage: Ye = QuaternionMatrixEuclideanJordanAlgebra.real_embed(Y)
1618 sage: XYe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1623 quaternions
= M
.base_ring()
1626 raise ValueError("the matrix 'M' must be square")
1628 F
= QuadraticField(-1, 'i')
1633 t
= z
.coefficient_tuple()
1638 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
1639 [-c
+ d
*i
, a
- b
*i
]])
1640 realM
= ComplexMatrixEuclideanJordanAlgebra
.real_embed(cplxM
)
1641 blocks
.append(realM
)
1643 # We should have real entries by now, so use the realest field
1644 # we've got for the return value.
1645 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
1650 def real_unembed(M
):
1652 The inverse of _embed_quaternion_matrix().
1656 sage: from mjo.eja.eja_algebra import \
1657 ....: QuaternionMatrixEuclideanJordanAlgebra
1661 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
1662 ....: [-2, 1, -4, 3],
1663 ....: [-3, 4, 1, -2],
1664 ....: [-4, -3, 2, 1]])
1665 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(M)
1666 [1 + 2*i + 3*j + 4*k]
1670 Unembedding is the inverse of embedding::
1672 sage: set_random_seed()
1673 sage: Q = QuaternionAlgebra(QQ, -1, -1)
1674 sage: M = random_matrix(Q, 3)
1675 sage: Me = QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1676 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1682 raise ValueError("the matrix 'M' must be square")
1683 if not n
.mod(4).is_zero():
1684 raise ValueError("the matrix 'M' must be a quaternion embedding")
1686 # Use the base ring of the matrix to ensure that its entries can be
1687 # multiplied by elements of the quaternion algebra.
1688 field
= M
.base_ring()
1689 Q
= QuaternionAlgebra(field
,-1,-1)
1692 # Go top-left to bottom-right (reading order), converting every
1693 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
1696 for l
in range(n
/4):
1697 for m
in range(n
/4):
1698 submat
= ComplexMatrixEuclideanJordanAlgebra
.real_unembed(
1699 M
[4*l
:4*l
+4,4*m
:4*m
+4] )
1700 if submat
[0,0] != submat
[1,1].conjugate():
1701 raise ValueError('bad on-diagonal submatrix')
1702 if submat
[0,1] != -submat
[1,0].conjugate():
1703 raise ValueError('bad off-diagonal submatrix')
1704 z
= submat
[0,0].vector()[0] # real part
1705 z
+= submat
[0,0].vector()[1]*i
# imag part
1706 z
+= submat
[0,1].vector()[0]*j
# real part
1707 z
+= submat
[0,1].vector()[1]*k
# imag part
1710 return matrix(Q
, n
/4, elements
)
1714 def natural_inner_product(cls
,X
,Y
):
1716 Compute a natural inner product in this algebra directly from
1721 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1725 This gives the same answer as the slow, default method implemented
1726 in :class:`MatrixEuclideanJordanAlgebra`::
1728 sage: set_random_seed()
1729 sage: J = QuaternionHermitianEJA.random_instance()
1730 sage: x,y = J.random_elements(2)
1731 sage: Xe = x.natural_representation()
1732 sage: Ye = y.natural_representation()
1733 sage: X = QuaternionHermitianEJA.real_unembed(Xe)
1734 sage: Y = QuaternionHermitianEJA.real_unembed(Ye)
1735 sage: expected = (X*Y).trace().coefficient_tuple()[0]
1736 sage: actual = QuaternionHermitianEJA.natural_inner_product(Xe,Ye)
1737 sage: actual == expected
1741 return RealMatrixEuclideanJordanAlgebra
.natural_inner_product(X
,Y
)/4
1744 class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra
,
1747 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
1748 matrices, the usual symmetric Jordan product, and the
1749 real-part-of-trace inner product. It has dimension `2n^2 - n` over
1754 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1758 In theory, our "field" can be any subfield of the reals::
1760 sage: QuaternionHermitianEJA(2, AA)
1761 Euclidean Jordan algebra of dimension 6 over Algebraic Real Field
1762 sage: QuaternionHermitianEJA(2, RR)
1763 Euclidean Jordan algebra of dimension 6 over Real Field with
1764 53 bits of precision
1768 The dimension of this algebra is `2*n^2 - n`::
1770 sage: set_random_seed()
1771 sage: n_max = QuaternionHermitianEJA._max_test_case_size()
1772 sage: n = ZZ.random_element(1, n_max)
1773 sage: J = QuaternionHermitianEJA(n)
1774 sage: J.dimension() == 2*(n^2) - n
1777 The Jordan multiplication is what we think it is::
1779 sage: set_random_seed()
1780 sage: J = QuaternionHermitianEJA.random_instance()
1781 sage: x,y = J.random_elements(2)
1782 sage: actual = (x*y).natural_representation()
1783 sage: X = x.natural_representation()
1784 sage: Y = y.natural_representation()
1785 sage: expected = (X*Y + Y*X)/2
1786 sage: actual == expected
1788 sage: J(expected) == x*y
1791 We can change the generator prefix::
1793 sage: QuaternionHermitianEJA(2, prefix='a').gens()
1794 (a0, a1, a2, a3, a4, a5)
1796 Our natural basis is normalized with respect to the natural inner
1797 product unless we specify otherwise::
1799 sage: set_random_seed()
1800 sage: J = QuaternionHermitianEJA.random_instance()
1801 sage: all( b.norm() == 1 for b in J.gens() )
1804 Since our natural basis is normalized with respect to the natural
1805 inner product, and since we know that this algebra is an EJA, any
1806 left-multiplication operator's matrix will be symmetric because
1807 natural->EJA basis representation is an isometry and within the EJA
1808 the operator is self-adjoint by the Jordan axiom::
1810 sage: set_random_seed()
1811 sage: x = QuaternionHermitianEJA.random_instance().random_element()
1812 sage: x.operator().matrix().is_symmetric()
1817 def _denormalized_basis(cls
, n
, field
):
1819 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
1821 Why do we embed these? Basically, because all of numerical
1822 linear algebra assumes that you're working with vectors consisting
1823 of `n` entries from a field and scalars from the same field. There's
1824 no way to tell SageMath that (for example) the vectors contain
1825 complex numbers, while the scalar field is real.
1829 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1833 sage: set_random_seed()
1834 sage: n = ZZ.random_element(1,5)
1835 sage: B = QuaternionHermitianEJA._denormalized_basis(n,QQ)
1836 sage: all( M.is_symmetric() for M in B )
1840 Q
= QuaternionAlgebra(QQ
,-1,-1)
1843 # This is like the symmetric case, but we need to be careful:
1845 # * We want conjugate-symmetry, not just symmetry.
1846 # * The diagonal will (as a result) be real.
1850 for j
in range(i
+1):
1851 Eij
= matrix(Q
, n
, lambda k
,l
: k
==i
and l
==j
)
1853 Sij
= cls
.real_embed(Eij
)
1856 # The second, third, and fourth ones have a minus
1857 # because they're conjugated.
1858 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
1860 Sij_I
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
1862 Sij_J
= cls
.real_embed(J
*Eij
- J
*Eij
.transpose())
1864 Sij_K
= cls
.real_embed(K
*Eij
- K
*Eij
.transpose())
1867 # Since we embedded these, we can drop back to the "field" that we
1868 # started with instead of the quaternion algebra "Q".
1869 return ( s
.change_ring(field
) for s
in S
)
1872 def __init__(self
, n
, field
=QQ
, **kwargs
):
1873 basis
= self
._denormalized
_basis
(n
,field
)
1874 super(QuaternionHermitianEJA
,self
).__init
__(field
, basis
, n
, **kwargs
)
1877 class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra
, KnownRankEJA
):
1879 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
1880 with the usual inner product and jordan product ``x*y =
1881 (<x_bar,y_bar>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
1886 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1890 This multiplication table can be verified by hand::
1892 sage: J = JordanSpinEJA(4)
1893 sage: e0,e1,e2,e3 = J.gens()
1909 We can change the generator prefix::
1911 sage: JordanSpinEJA(2, prefix='B').gens()
1915 def __init__(self
, n
, field
=QQ
, **kwargs
):
1916 V
= VectorSpace(field
, n
)
1917 mult_table
= [[V
.zero() for j
in range(n
)] for i
in range(n
)]
1927 z0
= x
.inner_product(y
)
1928 zbar
= y0
*xbar
+ x0
*ybar
1929 z
= V([z0
] + zbar
.list())
1930 mult_table
[i
][j
] = z
1932 # The rank of the spin algebra is two, unless we're in a
1933 # one-dimensional ambient space (because the rank is bounded by
1934 # the ambient dimension).
1935 fdeja
= super(JordanSpinEJA
, self
)
1936 return fdeja
.__init
__(field
, mult_table
, rank
=min(n
,2), **kwargs
)
1938 def inner_product(self
, x
, y
):
1940 Faster to reimplement than to use natural representations.
1944 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1948 Ensure that this is the usual inner product for the algebras
1951 sage: set_random_seed()
1952 sage: J = JordanSpinEJA.random_instance()
1953 sage: x,y = J.random_elements(2)
1954 sage: X = x.natural_representation()
1955 sage: Y = y.natural_representation()
1956 sage: x.inner_product(y) == J.natural_inner_product(X,Y)
1960 return x
.to_vector().inner_product(y
.to_vector())
1963 class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra
, KnownRankEJA
):
1965 The trivial Euclidean Jordan algebra consisting of only a zero element.
1969 sage: from mjo.eja.eja_algebra import TrivialEJA
1973 sage: J = TrivialEJA()
1980 sage: 7*J.one()*12*J.one()
1982 sage: J.one().inner_product(J.one())
1984 sage: J.one().norm()
1986 sage: J.one().subalgebra_generated_by()
1987 Euclidean Jordan algebra of dimension 0 over Rational Field
1992 def __init__(self
, field
=QQ
, **kwargs
):
1994 fdeja
= super(TrivialEJA
, self
)
1995 # The rank is zero using my definition, namely the dimension of the
1996 # largest subalgebra generated by any element.
1997 return fdeja
.__init
__(field
, mult_table
, rank
=0, **kwargs
)