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
= [
98 list(map(lambda x
: self
.from_vector(x
), ls
))
103 def _element_constructor_(self
, elt
):
105 Construct an element of this algebra from its natural
108 This gets called only after the parent element _call_ method
109 fails to find a coercion for the argument.
113 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
115 ....: RealSymmetricEJA)
119 The identity in `S^n` is converted to the identity in the EJA::
121 sage: J = RealSymmetricEJA(3)
122 sage: I = matrix.identity(QQ,3)
123 sage: J(I) == J.one()
126 This skew-symmetric matrix can't be represented in the EJA::
128 sage: J = RealSymmetricEJA(3)
129 sage: A = matrix(QQ,3, lambda i,j: i-j)
131 Traceback (most recent call last):
133 ArithmeticError: vector is not in free module
137 Ensure that we can convert any element of the two non-matrix
138 simple algebras (whose natural representations are their usual
139 vector representations) back and forth faithfully::
141 sage: set_random_seed()
142 sage: J = HadamardEJA.random_instance()
143 sage: x = J.random_element()
144 sage: J(x.to_vector().column()) == x
146 sage: J = JordanSpinEJA.random_instance()
147 sage: x = J.random_element()
148 sage: J(x.to_vector().column()) == x
153 # The superclass implementation of random_element()
154 # needs to be able to coerce "0" into the algebra.
157 natural_basis
= self
.natural_basis()
158 basis_space
= natural_basis
[0].matrix_space()
159 if elt
not in basis_space
:
160 raise ValueError("not a naturally-represented algebra element")
162 # Thanks for nothing! Matrix spaces aren't vector spaces in
163 # Sage, so we have to figure out its natural-basis coordinates
164 # ourselves. We use the basis space's ring instead of the
165 # element's ring because the basis space might be an algebraic
166 # closure whereas the base ring of the 3-by-3 identity matrix
167 # could be QQ instead of QQbar.
168 V
= VectorSpace(basis_space
.base_ring(), elt
.nrows()*elt
.ncols())
169 W
= V
.span_of_basis( _mat2vec(s
) for s
in natural_basis
)
170 coords
= W
.coordinate_vector(_mat2vec(elt
))
171 return self
.from_vector(coords
)
176 Return a string representation of ``self``.
180 sage: from mjo.eja.eja_algebra import JordanSpinEJA
184 Ensure that it says what we think it says::
186 sage: JordanSpinEJA(2, field=QQ)
187 Euclidean Jordan algebra of dimension 2 over Rational Field
188 sage: JordanSpinEJA(3, field=RDF)
189 Euclidean Jordan algebra of dimension 3 over Real Double Field
192 fmt
= "Euclidean Jordan algebra of dimension {} over {}"
193 return fmt
.format(self
.dimension(), self
.base_ring())
195 def product_on_basis(self
, i
, j
):
196 return self
._multiplication
_table
[i
][j
]
198 def _a_regular_element(self
):
200 Guess a regular element. Needed to compute the basis for our
201 characteristic polynomial coefficients.
205 sage: from mjo.eja.eja_algebra import random_eja
209 Ensure that this hacky method succeeds for every algebra that we
210 know how to construct::
212 sage: set_random_seed()
213 sage: J = random_eja()
214 sage: J._a_regular_element().is_regular()
219 z
= self
.sum( (i
+1)*gs
[i
] for i
in range(len(gs
)) )
220 if not z
.is_regular():
221 raise ValueError("don't know a regular element")
226 def _charpoly_basis_space(self
):
228 Return the vector space spanned by the basis used in our
229 characteristic polynomial coefficients. This is used not only to
230 compute those coefficients, but also any time we need to
231 evaluate the coefficients (like when we compute the trace or
234 z
= self
._a
_regular
_element
()
235 # Don't use the parent vector space directly here in case this
236 # happens to be a subalgebra. In that case, we would be e.g.
237 # two-dimensional but span_of_basis() would expect three
239 V
= VectorSpace(self
.base_ring(), self
.vector_space().dimension())
240 basis
= [ (z
**k
).to_vector() for k
in range(self
.rank()) ]
241 V1
= V
.span_of_basis( basis
)
242 b
= (V1
.basis() + V1
.complement().basis())
243 return V
.span_of_basis(b
)
248 def _charpoly_coeff(self
, i
):
250 Return the coefficient polynomial "a_{i}" of this algebra's
251 general characteristic polynomial.
253 Having this be a separate cached method lets us compute and
254 store the trace/determinant (a_{r-1} and a_{0} respectively)
255 separate from the entire characteristic polynomial.
257 (A_of_x
, x
, xr
, detA
) = self
._charpoly
_matrix
_system
()
258 R
= A_of_x
.base_ring()
263 # Guaranteed by theory
266 # Danger: the in-place modification is done for performance
267 # reasons (reconstructing a matrix with huge polynomial
268 # entries is slow), but I don't know how cached_method works,
269 # so it's highly possible that we're modifying some global
270 # list variable by reference, here. In other words, you
271 # probably shouldn't call this method twice on the same
272 # algebra, at the same time, in two threads
273 Ai_orig
= A_of_x
.column(i
)
274 A_of_x
.set_column(i
,xr
)
275 numerator
= A_of_x
.det()
276 A_of_x
.set_column(i
,Ai_orig
)
278 # We're relying on the theory here to ensure that each a_i is
279 # indeed back in R, and the added negative signs are to make
280 # the whole charpoly expression sum to zero.
281 return R(-numerator
/detA
)
285 def _charpoly_matrix_system(self
):
287 Compute the matrix whose entries A_ij are polynomials in
288 X1,...,XN, the vector ``x`` of variables X1,...,XN, the vector
289 corresponding to `x^r` and the determinent of the matrix A =
290 [A_ij]. In other words, all of the fixed (cachable) data needed
291 to compute the coefficients of the characteristic polynomial.
296 # Turn my vector space into a module so that "vectors" can
297 # have multivatiate polynomial entries.
298 names
= tuple('X' + str(i
) for i
in range(1,n
+1))
299 R
= PolynomialRing(self
.base_ring(), names
)
301 # Using change_ring() on the parent's vector space doesn't work
302 # here because, in a subalgebra, that vector space has a basis
303 # and change_ring() tries to bring the basis along with it. And
304 # that doesn't work unless the new ring is a PID, which it usually
308 # Now let x = (X1,X2,...,Xn) be the vector whose entries are
312 # And figure out the "left multiplication by x" matrix in
315 monomial_matrices
= [ self
.monomial(i
).operator().matrix()
316 for i
in range(n
) ] # don't recompute these!
318 ek
= self
.monomial(k
).to_vector()
320 sum( x
[i
]*(monomial_matrices
[i
]*ek
)
321 for i
in range(n
) ) )
322 Lx
= matrix
.column(R
, lmbx_cols
)
324 # Now we can compute powers of x "symbolically"
325 x_powers
= [self
.one().to_vector(), x
]
326 for d
in range(2, r
+1):
327 x_powers
.append( Lx
*(x_powers
[-1]) )
329 idmat
= matrix
.identity(R
, n
)
331 W
= self
._charpoly
_basis
_space
()
332 W
= W
.change_ring(R
.fraction_field())
334 # Starting with the standard coordinates x = (X1,X2,...,Xn)
335 # and then converting the entries to W-coordinates allows us
336 # to pass in the standard coordinates to the charpoly and get
337 # back the right answer. Specifically, with x = (X1,X2,...,Xn),
340 # W.coordinates(x^2) eval'd at (standard z-coords)
344 # W-coords of (standard coords of x^2 eval'd at std-coords of z)
346 # We want the middle equivalent thing in our matrix, but use
347 # the first equivalent thing instead so that we can pass in
348 # standard coordinates.
349 x_powers
= [ W
.coordinate_vector(xp
) for xp
in x_powers
]
350 l2
= [idmat
.column(k
-1) for k
in range(r
+1, n
+1)]
351 A_of_x
= matrix
.column(R
, n
, (x_powers
[:r
] + l2
))
352 return (A_of_x
, x
, x_powers
[r
], A_of_x
.det())
356 def characteristic_polynomial(self
):
358 Return a characteristic polynomial that works for all elements
361 The resulting polynomial has `n+1` variables, where `n` is the
362 dimension of this algebra. The first `n` variables correspond to
363 the coordinates of an algebra element: when evaluated at the
364 coordinates of an algebra element with respect to a certain
365 basis, the result is a univariate polynomial (in the one
366 remaining variable ``t``), namely the characteristic polynomial
371 sage: from mjo.eja.eja_algebra import JordanSpinEJA, TrivialEJA
375 The characteristic polynomial in the spin algebra is given in
376 Alizadeh, Example 11.11::
378 sage: J = JordanSpinEJA(3)
379 sage: p = J.characteristic_polynomial(); p
380 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
381 sage: xvec = J.one().to_vector()
385 By definition, the characteristic polynomial is a monic
386 degree-zero polynomial in a rank-zero algebra. Note that
387 Cayley-Hamilton is indeed satisfied since the polynomial
388 ``1`` evaluates to the identity element of the algebra on
391 sage: J = TrivialEJA()
392 sage: J.characteristic_polynomial()
399 # The list of coefficient polynomials a_0, a_1, a_2, ..., a_n.
400 a
= [ self
._charpoly
_coeff
(i
) for i
in range(r
+1) ]
402 # We go to a bit of trouble here to reorder the
403 # indeterminates, so that it's easier to evaluate the
404 # characteristic polynomial at x's coordinates and get back
405 # something in terms of t, which is what we want.
407 S
= PolynomialRing(self
.base_ring(),'t')
409 S
= PolynomialRing(S
, R
.variable_names())
412 return sum( a
[k
]*(t
**k
) for k
in range(len(a
)) )
415 def inner_product(self
, x
, y
):
417 The inner product associated with this Euclidean Jordan algebra.
419 Defaults to the trace inner product, but can be overridden by
420 subclasses if they are sure that the necessary properties are
425 sage: from mjo.eja.eja_algebra import random_eja
429 Our inner product is "associative," which means the following for
430 a symmetric bilinear form::
432 sage: set_random_seed()
433 sage: J = random_eja()
434 sage: x,y,z = J.random_elements(3)
435 sage: (x*y).inner_product(z) == y.inner_product(x*z)
439 X
= x
.natural_representation()
440 Y
= y
.natural_representation()
441 return self
.natural_inner_product(X
,Y
)
444 def is_trivial(self
):
446 Return whether or not this algebra is trivial.
448 A trivial algebra contains only the zero element.
452 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
457 sage: J = ComplexHermitianEJA(3)
463 sage: J = TrivialEJA()
468 return self
.dimension() == 0
471 def multiplication_table(self
):
473 Return a visual representation of this algebra's multiplication
474 table (on basis elements).
478 sage: from mjo.eja.eja_algebra import JordanSpinEJA
482 sage: J = JordanSpinEJA(4)
483 sage: J.multiplication_table()
484 +----++----+----+----+----+
485 | * || e0 | e1 | e2 | e3 |
486 +====++====+====+====+====+
487 | e0 || e0 | e1 | e2 | e3 |
488 +----++----+----+----+----+
489 | e1 || e1 | e0 | 0 | 0 |
490 +----++----+----+----+----+
491 | e2 || e2 | 0 | e0 | 0 |
492 +----++----+----+----+----+
493 | e3 || e3 | 0 | 0 | e0 |
494 +----++----+----+----+----+
497 M
= list(self
._multiplication
_table
) # copy
498 for i
in range(len(M
)):
499 # M had better be "square"
500 M
[i
] = [self
.monomial(i
)] + M
[i
]
501 M
= [["*"] + list(self
.gens())] + M
502 return table(M
, header_row
=True, header_column
=True, frame
=True)
505 def natural_basis(self
):
507 Return a more-natural representation of this algebra's basis.
509 Every finite-dimensional Euclidean Jordan Algebra is a direct
510 sum of five simple algebras, four of which comprise Hermitian
511 matrices. This method returns the original "natural" basis
512 for our underlying vector space. (Typically, the natural basis
513 is used to construct the multiplication table in the first place.)
515 Note that this will always return a matrix. The standard basis
516 in `R^n` will be returned as `n`-by-`1` column matrices.
520 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
521 ....: RealSymmetricEJA)
525 sage: J = RealSymmetricEJA(2)
527 Finite family {0: e0, 1: e1, 2: e2}
528 sage: J.natural_basis()
530 [1 0] [ 0 1/2*sqrt2] [0 0]
531 [0 0], [1/2*sqrt2 0], [0 1]
536 sage: J = JordanSpinEJA(2)
538 Finite family {0: e0, 1: e1}
539 sage: J.natural_basis()
546 if self
._natural
_basis
is None:
547 M
= self
.natural_basis_space()
548 return tuple( M(b
.to_vector()) for b
in self
.basis() )
550 return self
._natural
_basis
553 def natural_basis_space(self
):
555 Return the matrix space in which this algebra's natural basis
558 if self
._natural
_basis
is None or len(self
._natural
_basis
) == 0:
559 return MatrixSpace(self
.base_ring(), self
.dimension(), 1)
561 return self
._natural
_basis
[0].matrix_space()
565 def natural_inner_product(X
,Y
):
567 Compute the inner product of two naturally-represented elements.
569 For example in the real symmetric matrix EJA, this will compute
570 the trace inner-product of two n-by-n symmetric matrices. The
571 default should work for the real cartesian product EJA, the
572 Jordan spin EJA, and the real symmetric matrices. The others
573 will have to be overridden.
575 return (X
.conjugate_transpose()*Y
).trace()
581 Return the unit element of this algebra.
585 sage: from mjo.eja.eja_algebra import (HadamardEJA,
590 sage: J = HadamardEJA(5)
592 e0 + e1 + e2 + e3 + e4
596 The identity element acts like the identity::
598 sage: set_random_seed()
599 sage: J = random_eja()
600 sage: x = J.random_element()
601 sage: J.one()*x == x and x*J.one() == x
604 The matrix of the unit element's operator is the identity::
606 sage: set_random_seed()
607 sage: J = random_eja()
608 sage: actual = J.one().operator().matrix()
609 sage: expected = matrix.identity(J.base_ring(), J.dimension())
610 sage: actual == expected
614 # We can brute-force compute the matrices of the operators
615 # that correspond to the basis elements of this algebra.
616 # If some linear combination of those basis elements is the
617 # algebra identity, then the same linear combination of
618 # their matrices has to be the identity matrix.
620 # Of course, matrices aren't vectors in sage, so we have to
621 # appeal to the "long vectors" isometry.
622 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
624 # Now we use basis linear algebra to find the coefficients,
625 # of the matrices-as-vectors-linear-combination, which should
626 # work for the original algebra basis too.
627 A
= matrix
.column(self
.base_ring(), oper_vecs
)
629 # We used the isometry on the left-hand side already, but we
630 # still need to do it for the right-hand side. Recall that we
631 # wanted something that summed to the identity matrix.
632 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
634 # Now if there's an identity element in the algebra, this should work.
635 coeffs
= A
.solve_right(b
)
636 return self
.linear_combination(zip(self
.gens(), coeffs
))
639 def peirce_decomposition(self
, c
):
641 The Peirce decomposition of this algebra relative to the
644 In the future, this can be extended to a complete system of
645 orthogonal idempotents.
649 - ``c`` -- an idempotent of this algebra.
653 A triple (J0, J5, J1) containing two subalgebras and one subspace
656 - ``J0`` -- the algebra on the eigenspace of ``c.operator()``
657 corresponding to the eigenvalue zero.
659 - ``J5`` -- the eigenspace (NOT a subalgebra) of ``c.operator()``
660 corresponding to the eigenvalue one-half.
662 - ``J1`` -- the algebra on the eigenspace of ``c.operator()``
663 corresponding to the eigenvalue one.
665 These are the only possible eigenspaces for that operator, and this
666 algebra is a direct sum of them. The spaces ``J0`` and ``J1`` are
667 orthogonal, and are subalgebras of this algebra with the appropriate
672 sage: from mjo.eja.eja_algebra import random_eja, RealSymmetricEJA
676 The canonical example comes from the symmetric matrices, which
677 decompose into diagonal and off-diagonal parts::
679 sage: J = RealSymmetricEJA(3)
680 sage: C = matrix(QQ, [ [1,0,0],
684 sage: J0,J5,J1 = J.peirce_decomposition(c)
686 Euclidean Jordan algebra of dimension 1...
688 Vector space of degree 6 and dimension 2...
690 Euclidean Jordan algebra of dimension 3...
694 Every algebra decomposes trivially with respect to its identity
697 sage: set_random_seed()
698 sage: J = random_eja()
699 sage: J0,J5,J1 = J.peirce_decomposition(J.one())
700 sage: J0.dimension() == 0 and J5.dimension() == 0
702 sage: J1.superalgebra() == J and J1.dimension() == J.dimension()
705 The identity elements in the two subalgebras are the
706 projections onto their respective subspaces of the
707 superalgebra's identity element::
709 sage: set_random_seed()
710 sage: J = random_eja()
711 sage: x = J.random_element()
712 sage: if not J.is_trivial():
713 ....: while x.is_nilpotent():
714 ....: x = J.random_element()
715 sage: c = x.subalgebra_idempotent()
716 sage: J0,J5,J1 = J.peirce_decomposition(c)
717 sage: J1(c) == J1.one()
719 sage: J0(J.one() - c) == J0.one()
723 if not c
.is_idempotent():
724 raise ValueError("element is not idempotent: %s" % c
)
726 # Default these to what they should be if they turn out to be
727 # trivial, because eigenspaces_left() won't return eigenvalues
728 # corresponding to trivial spaces (e.g. it returns only the
729 # eigenspace corresponding to lambda=1 if you take the
730 # decomposition relative to the identity element).
731 trivial
= FiniteDimensionalEuclideanJordanSubalgebra(self
, ())
732 J0
= trivial
# eigenvalue zero
733 J5
= VectorSpace(self
.base_ring(), 0) # eigenvalue one-half
734 J1
= trivial
# eigenvalue one
736 for (eigval
, eigspace
) in c
.operator().matrix().left_eigenspaces():
737 if eigval
== ~
(self
.base_ring()(2)):
740 gens
= tuple( self
.from_vector(b
) for b
in eigspace
.basis() )
741 subalg
= FiniteDimensionalEuclideanJordanSubalgebra(self
, gens
)
747 raise ValueError("unexpected eigenvalue: %s" % eigval
)
752 def random_elements(self
, count
):
754 Return ``count`` random elements as a tuple.
758 sage: from mjo.eja.eja_algebra import JordanSpinEJA
762 sage: J = JordanSpinEJA(3)
763 sage: x,y,z = J.random_elements(3)
764 sage: all( [ x in J, y in J, z in J ])
766 sage: len( J.random_elements(10) ) == 10
770 return tuple( self
.random_element() for idx
in range(count
) )
775 Return the rank of this EJA.
779 The author knows of no algorithm to compute the rank of an EJA
780 where only the multiplication table is known. In lieu of one, we
781 require the rank to be specified when the algebra is created,
782 and simply pass along that number here.
786 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
787 ....: RealSymmetricEJA,
788 ....: ComplexHermitianEJA,
789 ....: QuaternionHermitianEJA,
794 The rank of the Jordan spin algebra is always two::
796 sage: JordanSpinEJA(2).rank()
798 sage: JordanSpinEJA(3).rank()
800 sage: JordanSpinEJA(4).rank()
803 The rank of the `n`-by-`n` Hermitian real, complex, or
804 quaternion matrices is `n`::
806 sage: RealSymmetricEJA(4).rank()
808 sage: ComplexHermitianEJA(3).rank()
810 sage: QuaternionHermitianEJA(2).rank()
815 Ensure that every EJA that we know how to construct has a
816 positive integer rank, unless the algebra is trivial in
817 which case its rank will be zero::
819 sage: set_random_seed()
820 sage: J = random_eja()
824 sage: r > 0 or (r == 0 and J.is_trivial())
831 def vector_space(self
):
833 Return the vector space that underlies this algebra.
837 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
841 sage: J = RealSymmetricEJA(2)
842 sage: J.vector_space()
843 Vector space of dimension 3 over...
846 return self
.zero().to_vector().parent().ambient_vector_space()
849 Element
= FiniteDimensionalEuclideanJordanAlgebraElement
852 class KnownRankEJA(object):
854 A class for algebras that we actually know we can construct. The
855 main issue is that, for most of our methods to make sense, we need
856 to know the rank of our algebra. Thus we can't simply generate a
857 "random" algebra, or even check that a given basis and product
858 satisfy the axioms; because even if everything looks OK, we wouldn't
859 know the rank we need to actuallty build the thing.
861 Not really a subclass of FDEJA because doing that causes method
862 resolution errors, e.g.
864 TypeError: Error when calling the metaclass bases
865 Cannot create a consistent method resolution
866 order (MRO) for bases FiniteDimensionalEuclideanJordanAlgebra,
871 def _max_test_case_size():
873 Return an integer "size" that is an upper bound on the size of
874 this algebra when it is used in a random test
875 case. Unfortunately, the term "size" is quite vague -- when
876 dealing with `R^n` under either the Hadamard or Jordan spin
877 product, the "size" refers to the dimension `n`. When dealing
878 with a matrix algebra (real symmetric or complex/quaternion
879 Hermitian), it refers to the size of the matrix, which is
880 far less than the dimension of the underlying vector space.
882 We default to five in this class, which is safe in `R^n`. The
883 matrix algebra subclasses (or any class where the "size" is
884 interpreted to be far less than the dimension) should override
885 with a smaller number.
890 def random_instance(cls
, field
=QQ
, **kwargs
):
892 Return a random instance of this type of algebra.
894 Beware, this will crash for "most instances" because the
895 constructor below looks wrong.
897 if cls
is TrivialEJA
:
898 # The TrivialEJA class doesn't take an "n" argument because
902 n
= ZZ
.random_element(cls
._max
_test
_case
_size
()) + 1
903 return cls(n
, field
, **kwargs
)
906 class HadamardEJA(FiniteDimensionalEuclideanJordanAlgebra
, KnownRankEJA
):
908 Return the Euclidean Jordan Algebra corresponding to the set
909 `R^n` under the Hadamard product.
911 Note: this is nothing more than the Cartesian product of ``n``
912 copies of the spin algebra. Once Cartesian product algebras
913 are implemented, this can go.
917 sage: from mjo.eja.eja_algebra import HadamardEJA
921 This multiplication table can be verified by hand::
923 sage: J = HadamardEJA(3)
924 sage: e0,e1,e2 = J.gens()
940 We can change the generator prefix::
942 sage: HadamardEJA(3, prefix='r').gens()
946 def __init__(self
, n
, field
=QQ
, **kwargs
):
947 V
= VectorSpace(field
, n
)
948 mult_table
= [ [ V
.gen(i
)*(i
== j
) for j
in range(n
) ]
951 fdeja
= super(HadamardEJA
, self
)
952 return fdeja
.__init
__(field
, mult_table
, rank
=n
, **kwargs
)
954 def inner_product(self
, x
, y
):
956 Faster to reimplement than to use natural representations.
960 sage: from mjo.eja.eja_algebra import HadamardEJA
964 Ensure that this is the usual inner product for the algebras
967 sage: set_random_seed()
968 sage: J = HadamardEJA.random_instance()
969 sage: x,y = J.random_elements(2)
970 sage: X = x.natural_representation()
971 sage: Y = y.natural_representation()
972 sage: x.inner_product(y) == J.natural_inner_product(X,Y)
976 return x
.to_vector().inner_product(y
.to_vector())
979 def random_eja(field
=QQ
, nontrivial
=False):
981 Return a "random" finite-dimensional Euclidean Jordan Algebra.
985 sage: from mjo.eja.eja_algebra import random_eja
990 Euclidean Jordan algebra of dimension...
993 eja_classes
= KnownRankEJA
.__subclasses
__()
995 eja_classes
.remove(TrivialEJA
)
996 classname
= choice(eja_classes
)
997 return classname
.random_instance(field
=field
)
1004 class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra
):
1006 def _max_test_case_size():
1007 # Play it safe, since this will be squared and the underlying
1008 # field can have dimension 4 (quaternions) too.
1011 def __init__(self
, field
, basis
, rank
, normalize_basis
=True, **kwargs
):
1013 Compared to the superclass constructor, we take a basis instead of
1014 a multiplication table because the latter can be computed in terms
1015 of the former when the product is known (like it is here).
1017 # Used in this class's fast _charpoly_coeff() override.
1018 self
._basis
_normalizers
= None
1020 # We're going to loop through this a few times, so now's a good
1021 # time to ensure that it isn't a generator expression.
1022 basis
= tuple(basis
)
1024 if rank
> 1 and normalize_basis
:
1025 # We'll need sqrt(2) to normalize the basis, and this
1026 # winds up in the multiplication table, so the whole
1027 # algebra needs to be over the field extension.
1028 R
= PolynomialRing(field
, 'z')
1031 if p
.is_irreducible():
1032 field
= field
.extension(p
, 'sqrt2', embedding
=RLF(2).sqrt())
1033 basis
= tuple( s
.change_ring(field
) for s
in basis
)
1034 self
._basis
_normalizers
= tuple(
1035 ~
(self
.natural_inner_product(s
,s
).sqrt()) for s
in basis
)
1036 basis
= tuple(s
*c
for (s
,c
) in zip(basis
,self
._basis
_normalizers
))
1038 Qs
= self
.multiplication_table_from_matrix_basis(basis
)
1040 fdeja
= super(MatrixEuclideanJordanAlgebra
, self
)
1041 return fdeja
.__init
__(field
,
1044 natural_basis
=basis
,
1049 def _charpoly_coeff(self
, i
):
1051 Override the parent method with something that tries to compute
1052 over a faster (non-extension) field.
1054 if self
._basis
_normalizers
is None:
1055 # We didn't normalize, so assume that the basis we started
1056 # with had entries in a nice field.
1057 return super(MatrixEuclideanJordanAlgebra
, self
)._charpoly
_coeff
(i
)
1059 basis
= ( (b
/n
) for (b
,n
) in zip(self
.natural_basis(),
1060 self
._basis
_normalizers
) )
1062 # Do this over the rationals and convert back at the end.
1063 J
= MatrixEuclideanJordanAlgebra(QQ
,
1066 normalize_basis
=False)
1067 (_
,x
,_
,_
) = J
._charpoly
_matrix
_system
()
1068 p
= J
._charpoly
_coeff
(i
)
1069 # p might be missing some vars, have to substitute "optionally"
1070 pairs
= zip(x
.base_ring().gens(), self
._basis
_normalizers
)
1071 substitutions
= { v: v*c for (v,c) in pairs }
1072 result
= p
.subs(substitutions
)
1074 # The result of "subs" can be either a coefficient-ring
1075 # element or a polynomial. Gotta handle both cases.
1077 return self
.base_ring()(result
)
1079 return result
.change_ring(self
.base_ring())
1083 def multiplication_table_from_matrix_basis(basis
):
1085 At least three of the five simple Euclidean Jordan algebras have the
1086 symmetric multiplication (A,B) |-> (AB + BA)/2, where the
1087 multiplication on the right is matrix multiplication. Given a basis
1088 for the underlying matrix space, this function returns a
1089 multiplication table (obtained by looping through the basis
1090 elements) for an algebra of those matrices.
1092 # In S^2, for example, we nominally have four coordinates even
1093 # though the space is of dimension three only. The vector space V
1094 # is supposed to hold the entire long vector, and the subspace W
1095 # of V will be spanned by the vectors that arise from symmetric
1096 # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
1097 field
= basis
[0].base_ring()
1098 dimension
= basis
[0].nrows()
1100 V
= VectorSpace(field
, dimension
**2)
1101 W
= V
.span_of_basis( _mat2vec(s
) for s
in basis
)
1103 mult_table
= [[W
.zero() for j
in range(n
)] for i
in range(n
)]
1106 mat_entry
= (basis
[i
]*basis
[j
] + basis
[j
]*basis
[i
])/2
1107 mult_table
[i
][j
] = W
.coordinate_vector(_mat2vec(mat_entry
))
1115 Embed the matrix ``M`` into a space of real matrices.
1117 The matrix ``M`` can have entries in any field at the moment:
1118 the real numbers, complex numbers, or quaternions. And although
1119 they are not a field, we can probably support octonions at some
1120 point, too. This function returns a real matrix that "acts like"
1121 the original with respect to matrix multiplication; i.e.
1123 real_embed(M*N) = real_embed(M)*real_embed(N)
1126 raise NotImplementedError
1130 def real_unembed(M
):
1132 The inverse of :meth:`real_embed`.
1134 raise NotImplementedError
1138 def natural_inner_product(cls
,X
,Y
):
1139 Xu
= cls
.real_unembed(X
)
1140 Yu
= cls
.real_unembed(Y
)
1141 tr
= (Xu
*Yu
).trace()
1144 # It's real already.
1147 # Otherwise, try the thing that works for complex numbers; and
1148 # if that doesn't work, the thing that works for quaternions.
1150 return tr
.vector()[0] # real part, imag part is index 1
1151 except AttributeError:
1152 # A quaternions doesn't have a vector() method, but does
1153 # have coefficient_tuple() method that returns the
1154 # coefficients of 1, i, j, and k -- in that order.
1155 return tr
.coefficient_tuple()[0]
1158 class RealMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1162 The identity function, for embedding real matrices into real
1168 def real_unembed(M
):
1170 The identity function, for unembedding real matrices from real
1176 class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra
, KnownRankEJA
):
1178 The rank-n simple EJA consisting of real symmetric n-by-n
1179 matrices, the usual symmetric Jordan product, and the trace inner
1180 product. It has dimension `(n^2 + n)/2` over the reals.
1184 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1188 sage: J = RealSymmetricEJA(2)
1189 sage: e0, e1, e2 = J.gens()
1197 In theory, our "field" can be any subfield of the reals::
1199 sage: RealSymmetricEJA(2, AA)
1200 Euclidean Jordan algebra of dimension 3 over Algebraic Real Field
1201 sage: RealSymmetricEJA(2, RR)
1202 Euclidean Jordan algebra of dimension 3 over Real Field with
1203 53 bits of precision
1207 The dimension of this algebra is `(n^2 + n) / 2`::
1209 sage: set_random_seed()
1210 sage: n_max = RealSymmetricEJA._max_test_case_size()
1211 sage: n = ZZ.random_element(1, n_max)
1212 sage: J = RealSymmetricEJA(n)
1213 sage: J.dimension() == (n^2 + n)/2
1216 The Jordan multiplication is what we think it is::
1218 sage: set_random_seed()
1219 sage: J = RealSymmetricEJA.random_instance()
1220 sage: x,y = J.random_elements(2)
1221 sage: actual = (x*y).natural_representation()
1222 sage: X = x.natural_representation()
1223 sage: Y = y.natural_representation()
1224 sage: expected = (X*Y + Y*X)/2
1225 sage: actual == expected
1227 sage: J(expected) == x*y
1230 We can change the generator prefix::
1232 sage: RealSymmetricEJA(3, prefix='q').gens()
1233 (q0, q1, q2, q3, q4, q5)
1235 Our natural basis is normalized with respect to the natural inner
1236 product unless we specify otherwise::
1238 sage: set_random_seed()
1239 sage: J = RealSymmetricEJA.random_instance()
1240 sage: all( b.norm() == 1 for b in J.gens() )
1243 Since our natural basis is normalized with respect to the natural
1244 inner product, and since we know that this algebra is an EJA, any
1245 left-multiplication operator's matrix will be symmetric because
1246 natural->EJA basis representation is an isometry and within the EJA
1247 the operator is self-adjoint by the Jordan axiom::
1249 sage: set_random_seed()
1250 sage: x = RealSymmetricEJA.random_instance().random_element()
1251 sage: x.operator().matrix().is_symmetric()
1256 def _denormalized_basis(cls
, n
, field
):
1258 Return a basis for the space of real symmetric n-by-n matrices.
1262 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1266 sage: set_random_seed()
1267 sage: n = ZZ.random_element(1,5)
1268 sage: B = RealSymmetricEJA._denormalized_basis(n,QQ)
1269 sage: all( M.is_symmetric() for M in B)
1273 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1277 for j
in range(i
+1):
1278 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
1282 Sij
= Eij
+ Eij
.transpose()
1288 def _max_test_case_size():
1289 return 4 # Dimension 10
1292 def __init__(self
, n
, field
=QQ
, **kwargs
):
1293 basis
= self
._denormalized
_basis
(n
, field
)
1294 super(RealSymmetricEJA
, self
).__init
__(field
, basis
, n
, **kwargs
)
1297 class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1301 Embed the n-by-n complex matrix ``M`` into the space of real
1302 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1303 bi` to the block matrix ``[[a,b],[-b,a]]``.
1307 sage: from mjo.eja.eja_algebra import \
1308 ....: ComplexMatrixEuclideanJordanAlgebra
1312 sage: F = QuadraticField(-1, 'i')
1313 sage: x1 = F(4 - 2*i)
1314 sage: x2 = F(1 + 2*i)
1317 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1318 sage: ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1327 Embedding is a homomorphism (isomorphism, in fact)::
1329 sage: set_random_seed()
1330 sage: n_max = ComplexMatrixEuclideanJordanAlgebra._max_test_case_size()
1331 sage: n = ZZ.random_element(n_max)
1332 sage: F = QuadraticField(-1, 'i')
1333 sage: X = random_matrix(F, n)
1334 sage: Y = random_matrix(F, n)
1335 sage: Xe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X)
1336 sage: Ye = ComplexMatrixEuclideanJordanAlgebra.real_embed(Y)
1337 sage: XYe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1344 raise ValueError("the matrix 'M' must be square")
1346 # We don't need any adjoined elements...
1347 field
= M
.base_ring().base_ring()
1351 a
= z
.list()[0] # real part, I guess
1352 b
= z
.list()[1] # imag part, I guess
1353 blocks
.append(matrix(field
, 2, [[a
,b
],[-b
,a
]]))
1355 return matrix
.block(field
, n
, blocks
)
1359 def real_unembed(M
):
1361 The inverse of _embed_complex_matrix().
1365 sage: from mjo.eja.eja_algebra import \
1366 ....: ComplexMatrixEuclideanJordanAlgebra
1370 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1371 ....: [-2, 1, -4, 3],
1372 ....: [ 9, 10, 11, 12],
1373 ....: [-10, 9, -12, 11] ])
1374 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(A)
1376 [ 10*i + 9 12*i + 11]
1380 Unembedding is the inverse of embedding::
1382 sage: set_random_seed()
1383 sage: F = QuadraticField(-1, 'i')
1384 sage: M = random_matrix(F, 3)
1385 sage: Me = ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1386 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1392 raise ValueError("the matrix 'M' must be square")
1393 if not n
.mod(2).is_zero():
1394 raise ValueError("the matrix 'M' must be a complex embedding")
1396 # If "M" was normalized, its base ring might have roots
1397 # adjoined and they can stick around after unembedding.
1398 field
= M
.base_ring()
1399 R
= PolynomialRing(field
, 'z')
1401 F
= field
.extension(z
**2 + 1, 'i', embedding
=CLF(-1).sqrt())
1404 # Go top-left to bottom-right (reading order), converting every
1405 # 2-by-2 block we see to a single complex element.
1407 for k
in range(n
/2):
1408 for j
in range(n
/2):
1409 submat
= M
[2*k
:2*k
+2,2*j
:2*j
+2]
1410 if submat
[0,0] != submat
[1,1]:
1411 raise ValueError('bad on-diagonal submatrix')
1412 if submat
[0,1] != -submat
[1,0]:
1413 raise ValueError('bad off-diagonal submatrix')
1414 z
= submat
[0,0] + submat
[0,1]*i
1417 return matrix(F
, n
/2, elements
)
1421 def natural_inner_product(cls
,X
,Y
):
1423 Compute a natural inner product in this algebra directly from
1428 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1432 This gives the same answer as the slow, default method implemented
1433 in :class:`MatrixEuclideanJordanAlgebra`::
1435 sage: set_random_seed()
1436 sage: J = ComplexHermitianEJA.random_instance()
1437 sage: x,y = J.random_elements(2)
1438 sage: Xe = x.natural_representation()
1439 sage: Ye = y.natural_representation()
1440 sage: X = ComplexHermitianEJA.real_unembed(Xe)
1441 sage: Y = ComplexHermitianEJA.real_unembed(Ye)
1442 sage: expected = (X*Y).trace().vector()[0]
1443 sage: actual = ComplexHermitianEJA.natural_inner_product(Xe,Ye)
1444 sage: actual == expected
1448 return RealMatrixEuclideanJordanAlgebra
.natural_inner_product(X
,Y
)/2
1451 class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra
, KnownRankEJA
):
1453 The rank-n simple EJA consisting of complex Hermitian n-by-n
1454 matrices over the real numbers, the usual symmetric Jordan product,
1455 and the real-part-of-trace inner product. It has dimension `n^2` over
1460 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1464 In theory, our "field" can be any subfield of the reals::
1466 sage: ComplexHermitianEJA(2, AA)
1467 Euclidean Jordan algebra of dimension 4 over Algebraic Real Field
1468 sage: ComplexHermitianEJA(2, RR)
1469 Euclidean Jordan algebra of dimension 4 over Real Field with
1470 53 bits of precision
1474 The dimension of this algebra is `n^2`::
1476 sage: set_random_seed()
1477 sage: n_max = ComplexHermitianEJA._max_test_case_size()
1478 sage: n = ZZ.random_element(1, n_max)
1479 sage: J = ComplexHermitianEJA(n)
1480 sage: J.dimension() == n^2
1483 The Jordan multiplication is what we think it is::
1485 sage: set_random_seed()
1486 sage: J = ComplexHermitianEJA.random_instance()
1487 sage: x,y = J.random_elements(2)
1488 sage: actual = (x*y).natural_representation()
1489 sage: X = x.natural_representation()
1490 sage: Y = y.natural_representation()
1491 sage: expected = (X*Y + Y*X)/2
1492 sage: actual == expected
1494 sage: J(expected) == x*y
1497 We can change the generator prefix::
1499 sage: ComplexHermitianEJA(2, prefix='z').gens()
1502 Our natural basis is normalized with respect to the natural inner
1503 product unless we specify otherwise::
1505 sage: set_random_seed()
1506 sage: J = ComplexHermitianEJA.random_instance()
1507 sage: all( b.norm() == 1 for b in J.gens() )
1510 Since our natural basis is normalized with respect to the natural
1511 inner product, and since we know that this algebra is an EJA, any
1512 left-multiplication operator's matrix will be symmetric because
1513 natural->EJA basis representation is an isometry and within the EJA
1514 the operator is self-adjoint by the Jordan axiom::
1516 sage: set_random_seed()
1517 sage: x = ComplexHermitianEJA.random_instance().random_element()
1518 sage: x.operator().matrix().is_symmetric()
1524 def _denormalized_basis(cls
, n
, field
):
1526 Returns a basis for the space of complex Hermitian n-by-n matrices.
1528 Why do we embed these? Basically, because all of numerical linear
1529 algebra assumes that you're working with vectors consisting of `n`
1530 entries from a field and scalars from the same field. There's no way
1531 to tell SageMath that (for example) the vectors contain complex
1532 numbers, while the scalar field is real.
1536 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1540 sage: set_random_seed()
1541 sage: n = ZZ.random_element(1,5)
1542 sage: field = QuadraticField(2, 'sqrt2')
1543 sage: B = ComplexHermitianEJA._denormalized_basis(n, field)
1544 sage: all( M.is_symmetric() for M in B)
1548 R
= PolynomialRing(field
, 'z')
1550 F
= field
.extension(z
**2 + 1, 'I')
1553 # This is like the symmetric case, but we need to be careful:
1555 # * We want conjugate-symmetry, not just symmetry.
1556 # * The diagonal will (as a result) be real.
1560 for j
in range(i
+1):
1561 Eij
= matrix(F
, n
, lambda k
,l
: k
==i
and l
==j
)
1563 Sij
= cls
.real_embed(Eij
)
1566 # The second one has a minus because it's conjugated.
1567 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
1569 Sij_imag
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
1572 # Since we embedded these, we can drop back to the "field" that we
1573 # started with instead of the complex extension "F".
1574 return ( s
.change_ring(field
) for s
in S
)
1577 def __init__(self
, n
, field
=QQ
, **kwargs
):
1578 basis
= self
._denormalized
_basis
(n
,field
)
1579 super(ComplexHermitianEJA
,self
).__init
__(field
, basis
, n
, **kwargs
)
1582 class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1586 Embed the n-by-n quaternion matrix ``M`` into the space of real
1587 matrices of size 4n-by-4n by first sending each quaternion entry `z
1588 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
1589 c+di],[-c + di, a-bi]]`, and then embedding those into a real
1594 sage: from mjo.eja.eja_algebra import \
1595 ....: QuaternionMatrixEuclideanJordanAlgebra
1599 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1600 sage: i,j,k = Q.gens()
1601 sage: x = 1 + 2*i + 3*j + 4*k
1602 sage: M = matrix(Q, 1, [[x]])
1603 sage: QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1609 Embedding is a homomorphism (isomorphism, in fact)::
1611 sage: set_random_seed()
1612 sage: n_max = QuaternionMatrixEuclideanJordanAlgebra._max_test_case_size()
1613 sage: n = ZZ.random_element(n_max)
1614 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1615 sage: X = random_matrix(Q, n)
1616 sage: Y = random_matrix(Q, n)
1617 sage: Xe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X)
1618 sage: Ye = QuaternionMatrixEuclideanJordanAlgebra.real_embed(Y)
1619 sage: XYe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1624 quaternions
= M
.base_ring()
1627 raise ValueError("the matrix 'M' must be square")
1629 F
= QuadraticField(-1, 'i')
1634 t
= z
.coefficient_tuple()
1639 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
1640 [-c
+ d
*i
, a
- b
*i
]])
1641 realM
= ComplexMatrixEuclideanJordanAlgebra
.real_embed(cplxM
)
1642 blocks
.append(realM
)
1644 # We should have real entries by now, so use the realest field
1645 # we've got for the return value.
1646 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
1651 def real_unembed(M
):
1653 The inverse of _embed_quaternion_matrix().
1657 sage: from mjo.eja.eja_algebra import \
1658 ....: QuaternionMatrixEuclideanJordanAlgebra
1662 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
1663 ....: [-2, 1, -4, 3],
1664 ....: [-3, 4, 1, -2],
1665 ....: [-4, -3, 2, 1]])
1666 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(M)
1667 [1 + 2*i + 3*j + 4*k]
1671 Unembedding is the inverse of embedding::
1673 sage: set_random_seed()
1674 sage: Q = QuaternionAlgebra(QQ, -1, -1)
1675 sage: M = random_matrix(Q, 3)
1676 sage: Me = QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1677 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1683 raise ValueError("the matrix 'M' must be square")
1684 if not n
.mod(4).is_zero():
1685 raise ValueError("the matrix 'M' must be a quaternion embedding")
1687 # Use the base ring of the matrix to ensure that its entries can be
1688 # multiplied by elements of the quaternion algebra.
1689 field
= M
.base_ring()
1690 Q
= QuaternionAlgebra(field
,-1,-1)
1693 # Go top-left to bottom-right (reading order), converting every
1694 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
1697 for l
in range(n
/4):
1698 for m
in range(n
/4):
1699 submat
= ComplexMatrixEuclideanJordanAlgebra
.real_unembed(
1700 M
[4*l
:4*l
+4,4*m
:4*m
+4] )
1701 if submat
[0,0] != submat
[1,1].conjugate():
1702 raise ValueError('bad on-diagonal submatrix')
1703 if submat
[0,1] != -submat
[1,0].conjugate():
1704 raise ValueError('bad off-diagonal submatrix')
1705 z
= submat
[0,0].vector()[0] # real part
1706 z
+= submat
[0,0].vector()[1]*i
# imag part
1707 z
+= submat
[0,1].vector()[0]*j
# real part
1708 z
+= submat
[0,1].vector()[1]*k
# imag part
1711 return matrix(Q
, n
/4, elements
)
1715 def natural_inner_product(cls
,X
,Y
):
1717 Compute a natural inner product in this algebra directly from
1722 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1726 This gives the same answer as the slow, default method implemented
1727 in :class:`MatrixEuclideanJordanAlgebra`::
1729 sage: set_random_seed()
1730 sage: J = QuaternionHermitianEJA.random_instance()
1731 sage: x,y = J.random_elements(2)
1732 sage: Xe = x.natural_representation()
1733 sage: Ye = y.natural_representation()
1734 sage: X = QuaternionHermitianEJA.real_unembed(Xe)
1735 sage: Y = QuaternionHermitianEJA.real_unembed(Ye)
1736 sage: expected = (X*Y).trace().coefficient_tuple()[0]
1737 sage: actual = QuaternionHermitianEJA.natural_inner_product(Xe,Ye)
1738 sage: actual == expected
1742 return RealMatrixEuclideanJordanAlgebra
.natural_inner_product(X
,Y
)/4
1745 class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra
,
1748 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
1749 matrices, the usual symmetric Jordan product, and the
1750 real-part-of-trace inner product. It has dimension `2n^2 - n` over
1755 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1759 In theory, our "field" can be any subfield of the reals::
1761 sage: QuaternionHermitianEJA(2, AA)
1762 Euclidean Jordan algebra of dimension 6 over Algebraic Real Field
1763 sage: QuaternionHermitianEJA(2, RR)
1764 Euclidean Jordan algebra of dimension 6 over Real Field with
1765 53 bits of precision
1769 The dimension of this algebra is `2*n^2 - n`::
1771 sage: set_random_seed()
1772 sage: n_max = QuaternionHermitianEJA._max_test_case_size()
1773 sage: n = ZZ.random_element(1, n_max)
1774 sage: J = QuaternionHermitianEJA(n)
1775 sage: J.dimension() == 2*(n^2) - n
1778 The Jordan multiplication is what we think it is::
1780 sage: set_random_seed()
1781 sage: J = QuaternionHermitianEJA.random_instance()
1782 sage: x,y = J.random_elements(2)
1783 sage: actual = (x*y).natural_representation()
1784 sage: X = x.natural_representation()
1785 sage: Y = y.natural_representation()
1786 sage: expected = (X*Y + Y*X)/2
1787 sage: actual == expected
1789 sage: J(expected) == x*y
1792 We can change the generator prefix::
1794 sage: QuaternionHermitianEJA(2, prefix='a').gens()
1795 (a0, a1, a2, a3, a4, a5)
1797 Our natural basis is normalized with respect to the natural inner
1798 product unless we specify otherwise::
1800 sage: set_random_seed()
1801 sage: J = QuaternionHermitianEJA.random_instance()
1802 sage: all( b.norm() == 1 for b in J.gens() )
1805 Since our natural basis is normalized with respect to the natural
1806 inner product, and since we know that this algebra is an EJA, any
1807 left-multiplication operator's matrix will be symmetric because
1808 natural->EJA basis representation is an isometry and within the EJA
1809 the operator is self-adjoint by the Jordan axiom::
1811 sage: set_random_seed()
1812 sage: x = QuaternionHermitianEJA.random_instance().random_element()
1813 sage: x.operator().matrix().is_symmetric()
1818 def _denormalized_basis(cls
, n
, field
):
1820 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
1822 Why do we embed these? Basically, because all of numerical
1823 linear algebra assumes that you're working with vectors consisting
1824 of `n` entries from a field and scalars from the same field. There's
1825 no way to tell SageMath that (for example) the vectors contain
1826 complex numbers, while the scalar field is real.
1830 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1834 sage: set_random_seed()
1835 sage: n = ZZ.random_element(1,5)
1836 sage: B = QuaternionHermitianEJA._denormalized_basis(n,QQ)
1837 sage: all( M.is_symmetric() for M in B )
1841 Q
= QuaternionAlgebra(QQ
,-1,-1)
1844 # This is like the symmetric case, but we need to be careful:
1846 # * We want conjugate-symmetry, not just symmetry.
1847 # * The diagonal will (as a result) be real.
1851 for j
in range(i
+1):
1852 Eij
= matrix(Q
, n
, lambda k
,l
: k
==i
and l
==j
)
1854 Sij
= cls
.real_embed(Eij
)
1857 # The second, third, and fourth ones have a minus
1858 # because they're conjugated.
1859 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
1861 Sij_I
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
1863 Sij_J
= cls
.real_embed(J
*Eij
- J
*Eij
.transpose())
1865 Sij_K
= cls
.real_embed(K
*Eij
- K
*Eij
.transpose())
1868 # Since we embedded these, we can drop back to the "field" that we
1869 # started with instead of the quaternion algebra "Q".
1870 return ( s
.change_ring(field
) for s
in S
)
1873 def __init__(self
, n
, field
=QQ
, **kwargs
):
1874 basis
= self
._denormalized
_basis
(n
,field
)
1875 super(QuaternionHermitianEJA
,self
).__init
__(field
, basis
, n
, **kwargs
)
1878 class BilinearFormEJA(FiniteDimensionalEuclideanJordanAlgebra
, KnownRankEJA
):
1880 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
1881 with the half-trace inner product and jordan product ``x*y =
1882 (x0*y0 + <B*x_bar,y_bar>, x0*y_bar + y0*x_bar)`` where ``B`` is a
1883 symmetric positive-definite "bilinear form" matrix. It has
1884 dimension `n` over the reals, and reduces to the ``JordanSpinEJA``
1885 when ``B`` is the identity matrix of order ``n-1``.
1889 sage: from mjo.eja.eja_algebra import BilinearFormEJA
1893 When no bilinear form is specified, the identity matrix is used,
1894 and the resulting algebra is the Jordan spin algebra::
1896 sage: J0 = BilinearFormEJA(3)
1897 sage: J1 = JordanSpinEJA(3)
1898 sage: J0.multiplication_table() == J0.multiplication_table()
1903 We can create a zero-dimensional algebra::
1905 sage: J = BilinearFormEJA(0)
1909 We can check the multiplication condition given in the Jordan, von
1910 Neumann, and Wigner paper (and also discussed on my "On the
1911 symmetry..." paper). Note that this relies heavily on the standard
1912 choice of basis, as does anything utilizing the bilinear form matrix::
1914 sage: set_random_seed()
1915 sage: n = ZZ.random_element(5)
1916 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
1917 sage: B = M.transpose()*M
1918 sage: J = BilinearFormEJA(n, B=B)
1919 sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
1920 sage: V = J.vector_space()
1921 sage: sis = [ J.from_vector(V([0] + (M.inverse()*ei).list()))
1922 ....: for ei in eis ]
1923 sage: actual = [ sis[i]*sis[j]
1924 ....: for i in range(n-1)
1925 ....: for j in range(n-1) ]
1926 sage: expected = [ J.one() if i == j else J.zero()
1927 ....: for i in range(n-1)
1928 ....: for j in range(n-1) ]
1929 sage: actual == expected
1932 def __init__(self
, n
, field
=QQ
, B
=None, **kwargs
):
1934 self
._B
= matrix
.identity(field
, max(0,n
-1))
1938 V
= VectorSpace(field
, n
)
1939 mult_table
= [[V
.zero() for j
in range(n
)] for i
in range(n
)]
1948 z0
= x0
*y0
+ (self
._B
*xbar
).inner_product(ybar
)
1949 zbar
= y0
*xbar
+ x0
*ybar
1950 z
= V([z0
] + zbar
.list())
1951 mult_table
[i
][j
] = z
1953 # The rank of this algebra is two, unless we're in a
1954 # one-dimensional ambient space (because the rank is bounded
1955 # by the ambient dimension).
1956 fdeja
= super(BilinearFormEJA
, self
)
1957 return fdeja
.__init
__(field
, mult_table
, rank
=min(n
,2), **kwargs
)
1959 def inner_product(self
, x
, y
):
1961 Half of the trace inner product.
1963 This is defined so that the special case of the Jordan spin
1964 algebra gets the usual inner product.
1968 sage: from mjo.eja.eja_algebra import BilinearFormEJA
1972 Ensure that this is one-half of the trace inner-product::
1974 sage: set_random_seed()
1975 sage: n = ZZ.random_element(5)
1976 sage: M = matrix.random(QQ, n-1, algorithm='unimodular')
1977 sage: B = M.transpose()*M
1978 sage: J = BilinearFormEJA(n, B=B)
1979 sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
1980 sage: V = J.vector_space()
1981 sage: sis = [ J.from_vector(V([0] + (M.inverse()*ei).list()))
1982 ....: for ei in eis ]
1983 sage: actual = [ sis[i]*sis[j]
1984 ....: for i in range(n-1)
1985 ....: for j in range(n-1) ]
1986 sage: expected = [ J.one() if i == j else J.zero()
1987 ....: for i in range(n-1)
1988 ....: for j in range(n-1) ]
1991 xvec
= x
.to_vector()
1993 yvec
= y
.to_vector()
1995 return x
[0]*y
[0] + (self
._B
*xbar
).inner_product(ybar
)
1997 class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra
, KnownRankEJA
):
1999 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2000 with the usual inner product and jordan product ``x*y =
2001 (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2006 sage: from mjo.eja.eja_algebra import JordanSpinEJA
2010 This multiplication table can be verified by hand::
2012 sage: J = JordanSpinEJA(4)
2013 sage: e0,e1,e2,e3 = J.gens()
2029 We can change the generator prefix::
2031 sage: JordanSpinEJA(2, prefix='B').gens()
2035 def __init__(self
, n
, field
=QQ
, **kwargs
):
2036 V
= VectorSpace(field
, n
)
2037 mult_table
= [[V
.zero() for j
in range(n
)] for i
in range(n
)]
2047 z0
= x
.inner_product(y
)
2048 zbar
= y0
*xbar
+ x0
*ybar
2049 z
= V([z0
] + zbar
.list())
2050 mult_table
[i
][j
] = z
2052 # The rank of the spin algebra is two, unless we're in a
2053 # one-dimensional ambient space (because the rank is bounded by
2054 # the ambient dimension).
2055 fdeja
= super(JordanSpinEJA
, self
)
2056 return fdeja
.__init
__(field
, mult_table
, rank
=min(n
,2), **kwargs
)
2058 def inner_product(self
, x
, y
):
2060 Faster to reimplement than to use natural representations.
2064 sage: from mjo.eja.eja_algebra import JordanSpinEJA
2068 Ensure that this is the usual inner product for the algebras
2071 sage: set_random_seed()
2072 sage: J = JordanSpinEJA.random_instance()
2073 sage: x,y = J.random_elements(2)
2074 sage: X = x.natural_representation()
2075 sage: Y = y.natural_representation()
2076 sage: x.inner_product(y) == J.natural_inner_product(X,Y)
2080 return x
.to_vector().inner_product(y
.to_vector())
2083 class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra
, KnownRankEJA
):
2085 The trivial Euclidean Jordan algebra consisting of only a zero element.
2089 sage: from mjo.eja.eja_algebra import TrivialEJA
2093 sage: J = TrivialEJA()
2100 sage: 7*J.one()*12*J.one()
2102 sage: J.one().inner_product(J.one())
2104 sage: J.one().norm()
2106 sage: J.one().subalgebra_generated_by()
2107 Euclidean Jordan algebra of dimension 0 over Rational Field
2112 def __init__(self
, field
=QQ
, **kwargs
):
2114 fdeja
= super(TrivialEJA
, self
)
2115 # The rank is zero using my definition, namely the dimension of the
2116 # largest subalgebra generated by any element.
2117 return fdeja
.__init
__(field
, mult_table
, rank
=0, **kwargs
)