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 sage
.algebras
.quatalg
.quaternion_algebra
import QuaternionAlgebra
9 from sage
.categories
.magmatic_algebras
import MagmaticAlgebras
10 from sage
.combinat
.free_module
import CombinatorialFreeModule
11 from sage
.matrix
.constructor
import matrix
12 from sage
.matrix
.matrix_space
import MatrixSpace
13 from sage
.misc
.cachefunc
import cached_method
14 from sage
.misc
.prandom
import choice
15 from sage
.misc
.table
import table
16 from sage
.modules
.free_module
import FreeModule
, VectorSpace
17 from sage
.rings
.integer_ring
import ZZ
18 from sage
.rings
.number_field
.number_field
import NumberField
, QuadraticField
19 from sage
.rings
.polynomial
.polynomial_ring_constructor
import PolynomialRing
20 from sage
.rings
.rational_field
import QQ
21 from sage
.rings
.real_lazy
import CLF
, RLF
22 from sage
.structure
.element
import is_Matrix
24 from mjo
.eja
.eja_element
import FiniteDimensionalEuclideanJordanAlgebraElement
25 from mjo
.eja
.eja_utils
import _mat2vec
27 class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule
):
28 # This is an ugly hack needed to prevent the category framework
29 # from implementing a coercion from our base ring (e.g. the
30 # rationals) into the algebra. First of all -- such a coercion is
31 # nonsense to begin with. But more importantly, it tries to do so
32 # in the category of rings, and since our algebras aren't
33 # associative they generally won't be rings.
34 _no_generic_basering_coercion
= True
46 sage: from mjo.eja.eja_algebra import random_eja
50 By definition, Jordan multiplication commutes::
52 sage: set_random_seed()
53 sage: J = random_eja()
54 sage: x = J.random_element()
55 sage: y = J.random_element()
61 self
._natural
_basis
= natural_basis
64 category
= MagmaticAlgebras(field
).FiniteDimensional()
65 category
= category
.WithBasis().Unital()
67 fda
= super(FiniteDimensionalEuclideanJordanAlgebra
, self
)
69 range(len(mult_table
)),
72 self
.print_options(bracket
='')
74 # The multiplication table we're given is necessarily in terms
75 # of vectors, because we don't have an algebra yet for
76 # anything to be an element of. However, it's faster in the
77 # long run to have the multiplication table be in terms of
78 # algebra elements. We do this after calling the superclass
79 # constructor so that from_vector() knows what to do.
80 self
._multiplication
_table
= [ map(lambda x
: self
.from_vector(x
), ls
)
81 for ls
in mult_table
]
84 def _element_constructor_(self
, elt
):
86 Construct an element of this algebra from its natural
89 This gets called only after the parent element _call_ method
90 fails to find a coercion for the argument.
94 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
95 ....: RealCartesianProductEJA,
96 ....: RealSymmetricEJA)
100 The identity in `S^n` is converted to the identity in the EJA::
102 sage: J = RealSymmetricEJA(3)
103 sage: I = matrix.identity(QQ,3)
104 sage: J(I) == J.one()
107 This skew-symmetric matrix can't be represented in the EJA::
109 sage: J = RealSymmetricEJA(3)
110 sage: A = matrix(QQ,3, lambda i,j: i-j)
112 Traceback (most recent call last):
114 ArithmeticError: vector is not in free module
118 Ensure that we can convert any element of the two non-matrix
119 simple algebras (whose natural representations are their usual
120 vector representations) back and forth faithfully::
122 sage: set_random_seed()
123 sage: J = RealCartesianProductEJA(5)
124 sage: x = J.random_element()
125 sage: J(x.to_vector().column()) == x
127 sage: J = JordanSpinEJA(5)
128 sage: x = J.random_element()
129 sage: J(x.to_vector().column()) == x
134 # The superclass implementation of random_element()
135 # needs to be able to coerce "0" into the algebra.
138 natural_basis
= self
.natural_basis()
139 basis_space
= natural_basis
[0].matrix_space()
140 if elt
not in basis_space
:
141 raise ValueError("not a naturally-represented algebra element")
143 # Thanks for nothing! Matrix spaces aren't vector spaces in
144 # Sage, so we have to figure out its natural-basis coordinates
145 # ourselves. We use the basis space's ring instead of the
146 # element's ring because the basis space might be an algebraic
147 # closure whereas the base ring of the 3-by-3 identity matrix
148 # could be QQ instead of QQbar.
149 V
= VectorSpace(basis_space
.base_ring(), elt
.nrows()*elt
.ncols())
150 W
= V
.span_of_basis( _mat2vec(s
) for s
in natural_basis
)
151 coords
= W
.coordinate_vector(_mat2vec(elt
))
152 return self
.from_vector(coords
)
157 Return a string representation of ``self``.
161 sage: from mjo.eja.eja_algebra import JordanSpinEJA
165 Ensure that it says what we think it says::
167 sage: JordanSpinEJA(2, field=QQ)
168 Euclidean Jordan algebra of dimension 2 over Rational Field
169 sage: JordanSpinEJA(3, field=RDF)
170 Euclidean Jordan algebra of dimension 3 over Real Double Field
173 fmt
= "Euclidean Jordan algebra of dimension {} over {}"
174 return fmt
.format(self
.dimension(), self
.base_ring())
176 def product_on_basis(self
, i
, j
):
177 return self
._multiplication
_table
[i
][j
]
179 def _a_regular_element(self
):
181 Guess a regular element. Needed to compute the basis for our
182 characteristic polynomial coefficients.
186 sage: from mjo.eja.eja_algebra import random_eja
190 Ensure that this hacky method succeeds for every algebra that we
191 know how to construct::
193 sage: set_random_seed()
194 sage: J = random_eja()
195 sage: J._a_regular_element().is_regular()
200 z
= self
.sum( (i
+1)*gs
[i
] for i
in range(len(gs
)) )
201 if not z
.is_regular():
202 raise ValueError("don't know a regular element")
207 def _charpoly_basis_space(self
):
209 Return the vector space spanned by the basis used in our
210 characteristic polynomial coefficients. This is used not only to
211 compute those coefficients, but also any time we need to
212 evaluate the coefficients (like when we compute the trace or
215 z
= self
._a
_regular
_element
()
216 # Don't use the parent vector space directly here in case this
217 # happens to be a subalgebra. In that case, we would be e.g.
218 # two-dimensional but span_of_basis() would expect three
220 V
= VectorSpace(self
.base_ring(), self
.vector_space().dimension())
221 basis
= [ (z
**k
).to_vector() for k
in range(self
.rank()) ]
222 V1
= V
.span_of_basis( basis
)
223 b
= (V1
.basis() + V1
.complement().basis())
224 return V
.span_of_basis(b
)
228 def _charpoly_coeff(self
, i
):
230 Return the coefficient polynomial "a_{i}" of this algebra's
231 general characteristic polynomial.
233 Having this be a separate cached method lets us compute and
234 store the trace/determinant (a_{r-1} and a_{0} respectively)
235 separate from the entire characteristic polynomial.
237 (A_of_x
, x
, xr
, detA
) = self
._charpoly
_matrix
_system
()
238 R
= A_of_x
.base_ring()
240 # Guaranteed by theory
243 # Danger: the in-place modification is done for performance
244 # reasons (reconstructing a matrix with huge polynomial
245 # entries is slow), but I don't know how cached_method works,
246 # so it's highly possible that we're modifying some global
247 # list variable by reference, here. In other words, you
248 # probably shouldn't call this method twice on the same
249 # algebra, at the same time, in two threads
250 Ai_orig
= A_of_x
.column(i
)
251 A_of_x
.set_column(i
,xr
)
252 numerator
= A_of_x
.det()
253 A_of_x
.set_column(i
,Ai_orig
)
255 # We're relying on the theory here to ensure that each a_i is
256 # indeed back in R, and the added negative signs are to make
257 # the whole charpoly expression sum to zero.
258 return R(-numerator
/detA
)
262 def _charpoly_matrix_system(self
):
264 Compute the matrix whose entries A_ij are polynomials in
265 X1,...,XN, the vector ``x`` of variables X1,...,XN, the vector
266 corresponding to `x^r` and the determinent of the matrix A =
267 [A_ij]. In other words, all of the fixed (cachable) data needed
268 to compute the coefficients of the characteristic polynomial.
273 # Turn my vector space into a module so that "vectors" can
274 # have multivatiate polynomial entries.
275 names
= tuple('X' + str(i
) for i
in range(1,n
+1))
276 R
= PolynomialRing(self
.base_ring(), names
)
278 # Using change_ring() on the parent's vector space doesn't work
279 # here because, in a subalgebra, that vector space has a basis
280 # and change_ring() tries to bring the basis along with it. And
281 # that doesn't work unless the new ring is a PID, which it usually
285 # Now let x = (X1,X2,...,Xn) be the vector whose entries are
289 # And figure out the "left multiplication by x" matrix in
292 monomial_matrices
= [ self
.monomial(i
).operator().matrix()
293 for i
in range(n
) ] # don't recompute these!
295 ek
= self
.monomial(k
).to_vector()
297 sum( x
[i
]*(monomial_matrices
[i
]*ek
)
298 for i
in range(n
) ) )
299 Lx
= matrix
.column(R
, lmbx_cols
)
301 # Now we can compute powers of x "symbolically"
302 x_powers
= [self
.one().to_vector(), x
]
303 for d
in range(2, r
+1):
304 x_powers
.append( Lx
*(x_powers
[-1]) )
306 idmat
= matrix
.identity(R
, n
)
308 W
= self
._charpoly
_basis
_space
()
309 W
= W
.change_ring(R
.fraction_field())
311 # Starting with the standard coordinates x = (X1,X2,...,Xn)
312 # and then converting the entries to W-coordinates allows us
313 # to pass in the standard coordinates to the charpoly and get
314 # back the right answer. Specifically, with x = (X1,X2,...,Xn),
317 # W.coordinates(x^2) eval'd at (standard z-coords)
321 # W-coords of (standard coords of x^2 eval'd at std-coords of z)
323 # We want the middle equivalent thing in our matrix, but use
324 # the first equivalent thing instead so that we can pass in
325 # standard coordinates.
326 x_powers
= [ W
.coordinate_vector(xp
) for xp
in x_powers
]
327 l2
= [idmat
.column(k
-1) for k
in range(r
+1, n
+1)]
328 A_of_x
= matrix
.column(R
, n
, (x_powers
[:r
] + l2
))
329 return (A_of_x
, x
, x_powers
[r
], A_of_x
.det())
333 def characteristic_polynomial(self
):
335 Return a characteristic polynomial that works for all elements
338 The resulting polynomial has `n+1` variables, where `n` is the
339 dimension of this algebra. The first `n` variables correspond to
340 the coordinates of an algebra element: when evaluated at the
341 coordinates of an algebra element with respect to a certain
342 basis, the result is a univariate polynomial (in the one
343 remaining variable ``t``), namely the characteristic polynomial
348 sage: from mjo.eja.eja_algebra import JordanSpinEJA
352 The characteristic polynomial in the spin algebra is given in
353 Alizadeh, Example 11.11::
355 sage: J = JordanSpinEJA(3)
356 sage: p = J.characteristic_polynomial(); p
357 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
358 sage: xvec = J.one().to_vector()
366 # The list of coefficient polynomials a_1, a_2, ..., a_n.
367 a
= [ self
._charpoly
_coeff
(i
) for i
in range(n
) ]
369 # We go to a bit of trouble here to reorder the
370 # indeterminates, so that it's easier to evaluate the
371 # characteristic polynomial at x's coordinates and get back
372 # something in terms of t, which is what we want.
374 S
= PolynomialRing(self
.base_ring(),'t')
376 S
= PolynomialRing(S
, R
.variable_names())
379 # Note: all entries past the rth should be zero. The
380 # coefficient of the highest power (x^r) is 1, but it doesn't
381 # appear in the solution vector which contains coefficients
382 # for the other powers (to make them sum to x^r).
384 a
[r
] = 1 # corresponds to x^r
386 # When the rank is equal to the dimension, trying to
387 # assign a[r] goes out-of-bounds.
388 a
.append(1) # corresponds to x^r
390 return sum( a
[k
]*(t
**k
) for k
in range(len(a
)) )
393 def inner_product(self
, x
, y
):
395 The inner product associated with this Euclidean Jordan algebra.
397 Defaults to the trace inner product, but can be overridden by
398 subclasses if they are sure that the necessary properties are
403 sage: from mjo.eja.eja_algebra import random_eja
407 The inner product must satisfy its axiom for this algebra to truly
408 be a Euclidean Jordan Algebra::
410 sage: set_random_seed()
411 sage: J = random_eja()
412 sage: x = J.random_element()
413 sage: y = J.random_element()
414 sage: z = J.random_element()
415 sage: (x*y).inner_product(z) == y.inner_product(x*z)
419 X
= x
.natural_representation()
420 Y
= y
.natural_representation()
421 return self
.__class
__.natural_inner_product(X
,Y
)
424 def is_trivial(self
):
426 Return whether or not this algebra is trivial.
428 A trivial algebra contains only the zero element.
432 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
436 sage: J = ComplexHermitianEJA(3)
439 sage: A = J.zero().subalgebra_generated_by()
444 return self
.dimension() == 0
447 def multiplication_table(self
):
449 Return a visual representation of this algebra's multiplication
450 table (on basis elements).
454 sage: from mjo.eja.eja_algebra import JordanSpinEJA
458 sage: J = JordanSpinEJA(4)
459 sage: J.multiplication_table()
460 +----++----+----+----+----+
461 | * || e0 | e1 | e2 | e3 |
462 +====++====+====+====+====+
463 | e0 || e0 | e1 | e2 | e3 |
464 +----++----+----+----+----+
465 | e1 || e1 | e0 | 0 | 0 |
466 +----++----+----+----+----+
467 | e2 || e2 | 0 | e0 | 0 |
468 +----++----+----+----+----+
469 | e3 || e3 | 0 | 0 | e0 |
470 +----++----+----+----+----+
473 M
= list(self
._multiplication
_table
) # copy
474 for i
in range(len(M
)):
475 # M had better be "square"
476 M
[i
] = [self
.monomial(i
)] + M
[i
]
477 M
= [["*"] + list(self
.gens())] + M
478 return table(M
, header_row
=True, header_column
=True, frame
=True)
481 def natural_basis(self
):
483 Return a more-natural representation of this algebra's basis.
485 Every finite-dimensional Euclidean Jordan Algebra is a direct
486 sum of five simple algebras, four of which comprise Hermitian
487 matrices. This method returns the original "natural" basis
488 for our underlying vector space. (Typically, the natural basis
489 is used to construct the multiplication table in the first place.)
491 Note that this will always return a matrix. The standard basis
492 in `R^n` will be returned as `n`-by-`1` column matrices.
496 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
497 ....: RealSymmetricEJA)
501 sage: J = RealSymmetricEJA(2)
503 Finite family {0: e0, 1: e1, 2: e2}
504 sage: J.natural_basis()
506 [1 0] [ 0 1/2*sqrt2] [0 0]
507 [0 0], [1/2*sqrt2 0], [0 1]
512 sage: J = JordanSpinEJA(2)
514 Finite family {0: e0, 1: e1}
515 sage: J.natural_basis()
522 if self
._natural
_basis
is None:
523 M
= self
.natural_basis_space()
524 return tuple( M(b
.to_vector()) for b
in self
.basis() )
526 return self
._natural
_basis
529 def natural_basis_space(self
):
531 Return the matrix space in which this algebra's natural basis
534 if self
._natural
_basis
is None or len(self
._natural
_basis
) == 0:
535 return MatrixSpace(self
.base_ring(), self
.dimension(), 1)
537 return self
._natural
_basis
[0].matrix_space()
541 def natural_inner_product(X
,Y
):
543 Compute the inner product of two naturally-represented elements.
545 For example in the real symmetric matrix EJA, this will compute
546 the trace inner-product of two n-by-n symmetric matrices. The
547 default should work for the real cartesian product EJA, the
548 Jordan spin EJA, and the real symmetric matrices. The others
549 will have to be overridden.
551 return (X
.conjugate_transpose()*Y
).trace()
557 Return the unit element of this algebra.
561 sage: from mjo.eja.eja_algebra import (RealCartesianProductEJA,
566 sage: J = RealCartesianProductEJA(5)
568 e0 + e1 + e2 + e3 + e4
572 The identity element acts like the identity::
574 sage: set_random_seed()
575 sage: J = random_eja()
576 sage: x = J.random_element()
577 sage: J.one()*x == x and x*J.one() == x
580 The matrix of the unit element's operator is the identity::
582 sage: set_random_seed()
583 sage: J = random_eja()
584 sage: actual = J.one().operator().matrix()
585 sage: expected = matrix.identity(J.base_ring(), J.dimension())
586 sage: actual == expected
590 # We can brute-force compute the matrices of the operators
591 # that correspond to the basis elements of this algebra.
592 # If some linear combination of those basis elements is the
593 # algebra identity, then the same linear combination of
594 # their matrices has to be the identity matrix.
596 # Of course, matrices aren't vectors in sage, so we have to
597 # appeal to the "long vectors" isometry.
598 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
600 # Now we use basis linear algebra to find the coefficients,
601 # of the matrices-as-vectors-linear-combination, which should
602 # work for the original algebra basis too.
603 A
= matrix
.column(self
.base_ring(), oper_vecs
)
605 # We used the isometry on the left-hand side already, but we
606 # still need to do it for the right-hand side. Recall that we
607 # wanted something that summed to the identity matrix.
608 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
610 # Now if there's an identity element in the algebra, this should work.
611 coeffs
= A
.solve_right(b
)
612 return self
.linear_combination(zip(self
.gens(), coeffs
))
615 def random_element(self
):
616 # Temporary workaround for https://trac.sagemath.org/ticket/28327
617 if self
.is_trivial():
620 s
= super(FiniteDimensionalEuclideanJordanAlgebra
, self
)
621 return s
.random_element()
626 Return the rank of this EJA.
630 The author knows of no algorithm to compute the rank of an EJA
631 where only the multiplication table is known. In lieu of one, we
632 require the rank to be specified when the algebra is created,
633 and simply pass along that number here.
637 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
638 ....: RealSymmetricEJA,
639 ....: ComplexHermitianEJA,
640 ....: QuaternionHermitianEJA,
645 The rank of the Jordan spin algebra is always two::
647 sage: JordanSpinEJA(2).rank()
649 sage: JordanSpinEJA(3).rank()
651 sage: JordanSpinEJA(4).rank()
654 The rank of the `n`-by-`n` Hermitian real, complex, or
655 quaternion matrices is `n`::
657 sage: RealSymmetricEJA(2).rank()
659 sage: ComplexHermitianEJA(2).rank()
661 sage: QuaternionHermitianEJA(2).rank()
663 sage: RealSymmetricEJA(5).rank()
665 sage: ComplexHermitianEJA(5).rank()
667 sage: QuaternionHermitianEJA(5).rank()
672 Ensure that every EJA that we know how to construct has a
673 positive integer rank::
675 sage: set_random_seed()
676 sage: r = random_eja().rank()
677 sage: r in ZZ and r > 0
684 def vector_space(self
):
686 Return the vector space that underlies this algebra.
690 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
694 sage: J = RealSymmetricEJA(2)
695 sage: J.vector_space()
696 Vector space of dimension 3 over...
699 return self
.zero().to_vector().parent().ambient_vector_space()
702 Element
= FiniteDimensionalEuclideanJordanAlgebraElement
705 class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra
):
707 Return the Euclidean Jordan Algebra corresponding to the set
708 `R^n` under the Hadamard product.
710 Note: this is nothing more than the Cartesian product of ``n``
711 copies of the spin algebra. Once Cartesian product algebras
712 are implemented, this can go.
716 sage: from mjo.eja.eja_algebra import RealCartesianProductEJA
720 This multiplication table can be verified by hand::
722 sage: J = RealCartesianProductEJA(3)
723 sage: e0,e1,e2 = J.gens()
739 We can change the generator prefix::
741 sage: RealCartesianProductEJA(3, prefix='r').gens()
744 Our inner product satisfies the Jordan axiom::
746 sage: set_random_seed()
747 sage: n = ZZ.random_element(1,5)
748 sage: J = RealCartesianProductEJA(n)
749 sage: x = J.random_element()
750 sage: y = J.random_element()
751 sage: z = J.random_element()
752 sage: (x*y).inner_product(z) == y.inner_product(x*z)
756 def __init__(self
, n
, field
=QQ
, **kwargs
):
757 V
= VectorSpace(field
, n
)
758 mult_table
= [ [ V
.gen(i
)*(i
== j
) for j
in range(n
) ]
761 fdeja
= super(RealCartesianProductEJA
, self
)
762 return fdeja
.__init
__(field
, mult_table
, rank
=n
, **kwargs
)
764 def inner_product(self
, x
, y
):
766 Faster to reimplement than to use natural representations.
770 sage: from mjo.eja.eja_algebra import RealCartesianProductEJA
774 Ensure that this is the usual inner product for the algebras
777 sage: set_random_seed()
778 sage: n = ZZ.random_element(1,5)
779 sage: J = RealCartesianProductEJA(n)
780 sage: x = J.random_element()
781 sage: y = J.random_element()
782 sage: X = x.natural_representation()
783 sage: Y = y.natural_representation()
784 sage: x.inner_product(y) == J.__class__.natural_inner_product(X,Y)
788 return x
.to_vector().inner_product(y
.to_vector())
793 Return a "random" finite-dimensional Euclidean Jordan Algebra.
797 For now, we choose a random natural number ``n`` (greater than zero)
798 and then give you back one of the following:
800 * The cartesian product of the rational numbers ``n`` times; this is
801 ``QQ^n`` with the Hadamard product.
803 * The Jordan spin algebra on ``QQ^n``.
805 * The ``n``-by-``n`` rational symmetric matrices with the symmetric
808 * The ``n``-by-``n`` complex-rational Hermitian matrices embedded
809 in the space of ``2n``-by-``2n`` real symmetric matrices.
811 * The ``n``-by-``n`` quaternion-rational Hermitian matrices embedded
812 in the space of ``4n``-by-``4n`` real symmetric matrices.
814 Later this might be extended to return Cartesian products of the
819 sage: from mjo.eja.eja_algebra import random_eja
824 Euclidean Jordan algebra of dimension...
828 # The max_n component lets us choose different upper bounds on the
829 # value "n" that gets passed to the constructor. This is needed
830 # because e.g. R^{10} is reasonable to test, while the Hermitian
831 # 10-by-10 quaternion matrices are not.
832 (constructor
, max_n
) = choice([(RealCartesianProductEJA
, 6),
834 (RealSymmetricEJA
, 5),
835 (ComplexHermitianEJA
, 4),
836 (QuaternionHermitianEJA
, 3)])
837 n
= ZZ
.random_element(1, max_n
)
838 return constructor(n
, field
=QQ
)
842 def _real_symmetric_basis(n
, field
):
844 Return a basis for the space of real symmetric n-by-n matrices.
848 sage: from mjo.eja.eja_algebra import _real_symmetric_basis
852 sage: set_random_seed()
853 sage: n = ZZ.random_element(1,5)
854 sage: B = _real_symmetric_basis(n, QQ)
855 sage: all( M.is_symmetric() for M in B)
859 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
863 for j
in xrange(i
+1):
864 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
868 Sij
= Eij
+ Eij
.transpose()
873 def _complex_hermitian_basis(n
, field
):
875 Returns a basis for the space of complex Hermitian n-by-n matrices.
877 Why do we embed these? Basically, because all of numerical linear
878 algebra assumes that you're working with vectors consisting of `n`
879 entries from a field and scalars from the same field. There's no way
880 to tell SageMath that (for example) the vectors contain complex
881 numbers, while the scalar field is real.
885 sage: from mjo.eja.eja_algebra import _complex_hermitian_basis
889 sage: set_random_seed()
890 sage: n = ZZ.random_element(1,5)
891 sage: field = QuadraticField(2, 'sqrt2')
892 sage: B = _complex_hermitian_basis(n, field)
893 sage: all( M.is_symmetric() for M in B)
897 R
= PolynomialRing(field
, 'z')
899 F
= NumberField(z
**2 + 1, 'I', embedding
=CLF(-1).sqrt())
902 # This is like the symmetric case, but we need to be careful:
904 # * We want conjugate-symmetry, not just symmetry.
905 # * The diagonal will (as a result) be real.
909 for j
in xrange(i
+1):
910 Eij
= matrix(F
, n
, lambda k
,l
: k
==i
and l
==j
)
912 Sij
= _embed_complex_matrix(Eij
)
915 # The second one has a minus because it's conjugated.
916 Sij_real
= _embed_complex_matrix(Eij
+ Eij
.transpose())
918 Sij_imag
= _embed_complex_matrix(I
*Eij
- I
*Eij
.transpose())
921 # Since we embedded these, we can drop back to the "field" that we
922 # started with instead of the complex extension "F".
923 return tuple( s
.change_ring(field
) for s
in S
)
927 def _quaternion_hermitian_basis(n
, field
, normalize
):
929 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
931 Why do we embed these? Basically, because all of numerical linear
932 algebra assumes that you're working with vectors consisting of `n`
933 entries from a field and scalars from the same field. There's no way
934 to tell SageMath that (for example) the vectors contain complex
935 numbers, while the scalar field is real.
939 sage: from mjo.eja.eja_algebra import _quaternion_hermitian_basis
943 sage: set_random_seed()
944 sage: n = ZZ.random_element(1,5)
945 sage: B = _quaternion_hermitian_basis(n, QQ, False)
946 sage: all( M.is_symmetric() for M in B )
950 Q
= QuaternionAlgebra(QQ
,-1,-1)
953 # This is like the symmetric case, but we need to be careful:
955 # * We want conjugate-symmetry, not just symmetry.
956 # * The diagonal will (as a result) be real.
960 for j
in xrange(i
+1):
961 Eij
= matrix(Q
, n
, lambda k
,l
: k
==i
and l
==j
)
963 Sij
= _embed_quaternion_matrix(Eij
)
966 # Beware, orthogonal but not normalized! The second,
967 # third, and fourth ones have a minus because they're
969 Sij_real
= _embed_quaternion_matrix(Eij
+ Eij
.transpose())
971 Sij_I
= _embed_quaternion_matrix(I
*Eij
- I
*Eij
.transpose())
973 Sij_J
= _embed_quaternion_matrix(J
*Eij
- J
*Eij
.transpose())
975 Sij_K
= _embed_quaternion_matrix(K
*Eij
- K
*Eij
.transpose())
981 def _multiplication_table_from_matrix_basis(basis
):
983 At least three of the five simple Euclidean Jordan algebras have the
984 symmetric multiplication (A,B) |-> (AB + BA)/2, where the
985 multiplication on the right is matrix multiplication. Given a basis
986 for the underlying matrix space, this function returns a
987 multiplication table (obtained by looping through the basis
988 elements) for an algebra of those matrices.
990 # In S^2, for example, we nominally have four coordinates even
991 # though the space is of dimension three only. The vector space V
992 # is supposed to hold the entire long vector, and the subspace W
993 # of V will be spanned by the vectors that arise from symmetric
994 # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
995 field
= basis
[0].base_ring()
996 dimension
= basis
[0].nrows()
998 V
= VectorSpace(field
, dimension
**2)
999 W
= V
.span_of_basis( _mat2vec(s
) for s
in basis
)
1001 mult_table
= [[W
.zero() for j
in range(n
)] for i
in range(n
)]
1004 mat_entry
= (basis
[i
]*basis
[j
] + basis
[j
]*basis
[i
])/2
1005 mult_table
[i
][j
] = W
.coordinate_vector(_mat2vec(mat_entry
))
1010 def _embed_complex_matrix(M
):
1012 Embed the n-by-n complex matrix ``M`` into the space of real
1013 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1014 bi` to the block matrix ``[[a,b],[-b,a]]``.
1018 sage: from mjo.eja.eja_algebra import _embed_complex_matrix
1022 sage: F = QuadraticField(-1, 'i')
1023 sage: x1 = F(4 - 2*i)
1024 sage: x2 = F(1 + 2*i)
1027 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1028 sage: _embed_complex_matrix(M)
1037 Embedding is a homomorphism (isomorphism, in fact)::
1039 sage: set_random_seed()
1040 sage: n = ZZ.random_element(5)
1041 sage: F = QuadraticField(-1, 'i')
1042 sage: X = random_matrix(F, n)
1043 sage: Y = random_matrix(F, n)
1044 sage: actual = _embed_complex_matrix(X) * _embed_complex_matrix(Y)
1045 sage: expected = _embed_complex_matrix(X*Y)
1046 sage: actual == expected
1052 raise ValueError("the matrix 'M' must be square")
1053 field
= M
.base_ring()
1056 a
= z
.vector()[0] # real part, I guess
1057 b
= z
.vector()[1] # imag part, I guess
1058 blocks
.append(matrix(field
, 2, [[a
,b
],[-b
,a
]]))
1060 # We can drop the imaginaries here.
1061 return matrix
.block(field
.base_ring(), n
, blocks
)
1064 def _unembed_complex_matrix(M
):
1066 The inverse of _embed_complex_matrix().
1070 sage: from mjo.eja.eja_algebra import (_embed_complex_matrix,
1071 ....: _unembed_complex_matrix)
1075 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1076 ....: [-2, 1, -4, 3],
1077 ....: [ 9, 10, 11, 12],
1078 ....: [-10, 9, -12, 11] ])
1079 sage: _unembed_complex_matrix(A)
1081 [ 10*i + 9 12*i + 11]
1085 Unembedding is the inverse of embedding::
1087 sage: set_random_seed()
1088 sage: F = QuadraticField(-1, 'i')
1089 sage: M = random_matrix(F, 3)
1090 sage: _unembed_complex_matrix(_embed_complex_matrix(M)) == M
1096 raise ValueError("the matrix 'M' must be square")
1097 if not n
.mod(2).is_zero():
1098 raise ValueError("the matrix 'M' must be a complex embedding")
1100 field
= M
.base_ring() # This should already have sqrt2
1101 R
= PolynomialRing(field
, 'z')
1103 F
= NumberField(z
**2 + 1,'i', embedding
=CLF(-1).sqrt())
1106 # Go top-left to bottom-right (reading order), converting every
1107 # 2-by-2 block we see to a single complex element.
1109 for k
in xrange(n
/2):
1110 for j
in xrange(n
/2):
1111 submat
= M
[2*k
:2*k
+2,2*j
:2*j
+2]
1112 if submat
[0,0] != submat
[1,1]:
1113 raise ValueError('bad on-diagonal submatrix')
1114 if submat
[0,1] != -submat
[1,0]:
1115 raise ValueError('bad off-diagonal submatrix')
1116 z
= submat
[0,0] + submat
[0,1]*i
1119 return matrix(F
, n
/2, elements
)
1122 def _embed_quaternion_matrix(M
):
1124 Embed the n-by-n quaternion matrix ``M`` into the space of real
1125 matrices of size 4n-by-4n by first sending each quaternion entry
1126 `z = a + bi + cj + dk` to the block-complex matrix
1127 ``[[a + bi, c+di],[-c + di, a-bi]]`, and then embedding those into
1132 sage: from mjo.eja.eja_algebra import _embed_quaternion_matrix
1136 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1137 sage: i,j,k = Q.gens()
1138 sage: x = 1 + 2*i + 3*j + 4*k
1139 sage: M = matrix(Q, 1, [[x]])
1140 sage: _embed_quaternion_matrix(M)
1146 Embedding is a homomorphism (isomorphism, in fact)::
1148 sage: set_random_seed()
1149 sage: n = ZZ.random_element(5)
1150 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1151 sage: X = random_matrix(Q, n)
1152 sage: Y = random_matrix(Q, n)
1153 sage: actual = _embed_quaternion_matrix(X)*_embed_quaternion_matrix(Y)
1154 sage: expected = _embed_quaternion_matrix(X*Y)
1155 sage: actual == expected
1159 quaternions
= M
.base_ring()
1162 raise ValueError("the matrix 'M' must be square")
1164 F
= QuadraticField(-1, 'i')
1169 t
= z
.coefficient_tuple()
1174 cplx_matrix
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
1175 [-c
+ d
*i
, a
- b
*i
]])
1176 blocks
.append(_embed_complex_matrix(cplx_matrix
))
1178 # We should have real entries by now, so use the realest field
1179 # we've got for the return value.
1180 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
1183 def _unembed_quaternion_matrix(M
):
1185 The inverse of _embed_quaternion_matrix().
1189 sage: from mjo.eja.eja_algebra import (_embed_quaternion_matrix,
1190 ....: _unembed_quaternion_matrix)
1194 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
1195 ....: [-2, 1, -4, 3],
1196 ....: [-3, 4, 1, -2],
1197 ....: [-4, -3, 2, 1]])
1198 sage: _unembed_quaternion_matrix(M)
1199 [1 + 2*i + 3*j + 4*k]
1203 Unembedding is the inverse of embedding::
1205 sage: set_random_seed()
1206 sage: Q = QuaternionAlgebra(QQ, -1, -1)
1207 sage: M = random_matrix(Q, 3)
1208 sage: _unembed_quaternion_matrix(_embed_quaternion_matrix(M)) == M
1214 raise ValueError("the matrix 'M' must be square")
1215 if not n
.mod(4).is_zero():
1216 raise ValueError("the matrix 'M' must be a complex embedding")
1218 Q
= QuaternionAlgebra(QQ
,-1,-1)
1221 # Go top-left to bottom-right (reading order), converting every
1222 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
1225 for l
in xrange(n
/4):
1226 for m
in xrange(n
/4):
1227 submat
= _unembed_complex_matrix(M
[4*l
:4*l
+4,4*m
:4*m
+4])
1228 if submat
[0,0] != submat
[1,1].conjugate():
1229 raise ValueError('bad on-diagonal submatrix')
1230 if submat
[0,1] != -submat
[1,0].conjugate():
1231 raise ValueError('bad off-diagonal submatrix')
1232 z
= submat
[0,0].real() + submat
[0,0].imag()*i
1233 z
+= submat
[0,1].real()*j
+ submat
[0,1].imag()*k
1236 return matrix(Q
, n
/4, elements
)
1239 # The inner product used for the real symmetric simple EJA.
1240 # We keep it as a separate function because e.g. the complex
1241 # algebra uses the same inner product, except divided by 2.
1242 def _matrix_ip(X
,Y
):
1243 X_mat
= X
.natural_representation()
1244 Y_mat
= Y
.natural_representation()
1245 return (X_mat
*Y_mat
).trace()
1248 class RealSymmetricEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1250 The rank-n simple EJA consisting of real symmetric n-by-n
1251 matrices, the usual symmetric Jordan product, and the trace inner
1252 product. It has dimension `(n^2 + n)/2` over the reals.
1256 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1260 sage: J = RealSymmetricEJA(2)
1261 sage: e0, e1, e2 = J.gens()
1271 The dimension of this algebra is `(n^2 + n) / 2`::
1273 sage: set_random_seed()
1274 sage: n = ZZ.random_element(1,5)
1275 sage: J = RealSymmetricEJA(n)
1276 sage: J.dimension() == (n^2 + n)/2
1279 The Jordan multiplication is what we think it is::
1281 sage: set_random_seed()
1282 sage: n = ZZ.random_element(1,5)
1283 sage: J = RealSymmetricEJA(n)
1284 sage: x = J.random_element()
1285 sage: y = J.random_element()
1286 sage: actual = (x*y).natural_representation()
1287 sage: X = x.natural_representation()
1288 sage: Y = y.natural_representation()
1289 sage: expected = (X*Y + Y*X)/2
1290 sage: actual == expected
1292 sage: J(expected) == x*y
1295 We can change the generator prefix::
1297 sage: RealSymmetricEJA(3, prefix='q').gens()
1298 (q0, q1, q2, q3, q4, q5)
1300 Our inner product satisfies the Jordan axiom::
1302 sage: set_random_seed()
1303 sage: n = ZZ.random_element(1,5)
1304 sage: J = RealSymmetricEJA(n)
1305 sage: x = J.random_element()
1306 sage: y = J.random_element()
1307 sage: z = J.random_element()
1308 sage: (x*y).inner_product(z) == y.inner_product(x*z)
1311 Our basis is normalized with respect to the natural inner product::
1313 sage: set_random_seed()
1314 sage: n = ZZ.random_element(1,5)
1315 sage: J = RealSymmetricEJA(n)
1316 sage: all( b.norm() == 1 for b in J.gens() )
1319 Left-multiplication operators are symmetric because they satisfy
1322 sage: set_random_seed()
1323 sage: n = ZZ.random_element(1,5)
1324 sage: x = RealSymmetricEJA(n).random_element()
1325 sage: x.operator().matrix().is_symmetric()
1329 def __init__(self
, n
, field
=QQ
, normalize_basis
=True, **kwargs
):
1330 S
= _real_symmetric_basis(n
, field
)
1332 if n
> 1 and normalize_basis
:
1333 # We'll need sqrt(2) to normalize the basis, and this
1334 # winds up in the multiplication table, so the whole
1335 # algebra needs to be over the field extension.
1336 R
= PolynomialRing(field
, 'z')
1339 if p
.is_irreducible():
1340 field
= NumberField(p
, 'sqrt2', embedding
=RLF(2).sqrt())
1341 S
= [ s
.change_ring(field
) for s
in S
]
1342 self
._basis
_denormalizers
= tuple(
1343 self
.__class
__.natural_inner_product(s
,s
).sqrt()
1345 S
= tuple( s
/c
for (s
,c
) in zip(S
,self
._basis
_denormalizers
) )
1347 Qs
= _multiplication_table_from_matrix_basis(S
)
1349 fdeja
= super(RealSymmetricEJA
, self
)
1350 return fdeja
.__init
__(field
,
1358 class ComplexHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1360 The rank-n simple EJA consisting of complex Hermitian n-by-n
1361 matrices over the real numbers, the usual symmetric Jordan product,
1362 and the real-part-of-trace inner product. It has dimension `n^2` over
1367 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1371 The dimension of this algebra is `n^2`::
1373 sage: set_random_seed()
1374 sage: n = ZZ.random_element(1,5)
1375 sage: J = ComplexHermitianEJA(n)
1376 sage: J.dimension() == n^2
1379 The Jordan multiplication is what we think it is::
1381 sage: set_random_seed()
1382 sage: n = ZZ.random_element(1,5)
1383 sage: J = ComplexHermitianEJA(n)
1384 sage: x = J.random_element()
1385 sage: y = J.random_element()
1386 sage: actual = (x*y).natural_representation()
1387 sage: X = x.natural_representation()
1388 sage: Y = y.natural_representation()
1389 sage: expected = (X*Y + Y*X)/2
1390 sage: actual == expected
1392 sage: J(expected) == x*y
1395 We can change the generator prefix::
1397 sage: ComplexHermitianEJA(2, prefix='z').gens()
1400 Our inner product satisfies the Jordan axiom::
1402 sage: set_random_seed()
1403 sage: n = ZZ.random_element(1,5)
1404 sage: J = ComplexHermitianEJA(n)
1405 sage: x = J.random_element()
1406 sage: y = J.random_element()
1407 sage: z = J.random_element()
1408 sage: (x*y).inner_product(z) == y.inner_product(x*z)
1411 Our basis is normalized with respect to the natural inner product::
1413 sage: set_random_seed()
1414 sage: n = ZZ.random_element(1,4)
1415 sage: J = ComplexHermitianEJA(n)
1416 sage: all( b.norm() == 1 for b in J.gens() )
1419 Left-multiplication operators are symmetric because they satisfy
1422 sage: set_random_seed()
1423 sage: n = ZZ.random_element(1,5)
1424 sage: x = ComplexHermitianEJA(n).random_element()
1425 sage: x.operator().matrix().is_symmetric()
1429 def __init__(self
, n
, field
=QQ
, normalize_basis
=True, **kwargs
):
1430 S
= _complex_hermitian_basis(n
, field
)
1432 if n
> 1 and normalize_basis
:
1433 # We'll need sqrt(2) to normalize the basis, and this
1434 # winds up in the multiplication table, so the whole
1435 # algebra needs to be over the field extension.
1436 R
= PolynomialRing(field
, 'z')
1439 if p
.is_irreducible():
1440 field
= NumberField(p
, 'sqrt2', embedding
=RLF(2).sqrt())
1441 S
= [ s
.change_ring(field
) for s
in S
]
1442 self
._basis
_denormalizers
= tuple(
1443 self
.__class
__.natural_inner_product(s
,s
).sqrt()
1445 S
= tuple( s
/c
for (s
,c
) in zip(S
,self
._basis
_denormalizers
) )
1447 Qs
= _multiplication_table_from_matrix_basis(S
)
1449 fdeja
= super(ComplexHermitianEJA
, self
)
1450 return fdeja
.__init
__(field
,
1458 def natural_inner_product(X
,Y
):
1459 Xu
= _unembed_complex_matrix(X
)
1460 Yu
= _unembed_complex_matrix(Y
)
1461 # The trace need not be real; consider Xu = (i*I) and Yu = I.
1462 return ((Xu
*Yu
).trace()).vector()[0] # real part, I guess
1464 class QuaternionHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1466 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
1467 matrices, the usual symmetric Jordan product, and the
1468 real-part-of-trace inner product. It has dimension `2n^2 - n` over
1473 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1477 The dimension of this algebra is `n^2`::
1479 sage: set_random_seed()
1480 sage: n = ZZ.random_element(1,5)
1481 sage: J = QuaternionHermitianEJA(n)
1482 sage: J.dimension() == 2*(n^2) - n
1485 The Jordan multiplication is what we think it is::
1487 sage: set_random_seed()
1488 sage: n = ZZ.random_element(1,5)
1489 sage: J = QuaternionHermitianEJA(n)
1490 sage: x = J.random_element()
1491 sage: y = J.random_element()
1492 sage: actual = (x*y).natural_representation()
1493 sage: X = x.natural_representation()
1494 sage: Y = y.natural_representation()
1495 sage: expected = (X*Y + Y*X)/2
1496 sage: actual == expected
1498 sage: J(expected) == x*y
1501 We can change the generator prefix::
1503 sage: QuaternionHermitianEJA(2, prefix='a').gens()
1504 (a0, a1, a2, a3, a4, a5)
1506 Our inner product satisfies the Jordan axiom::
1508 sage: set_random_seed()
1509 sage: n = ZZ.random_element(1,5)
1510 sage: J = QuaternionHermitianEJA(n)
1511 sage: x = J.random_element()
1512 sage: y = J.random_element()
1513 sage: z = J.random_element()
1514 sage: (x*y).inner_product(z) == y.inner_product(x*z)
1518 def __init__(self
, n
, field
=QQ
, normalize_basis
=True, **kwargs
):
1519 S
= _quaternion_hermitian_basis(n
, field
, normalize_basis
)
1520 Qs
= _multiplication_table_from_matrix_basis(S
)
1522 fdeja
= super(QuaternionHermitianEJA
, self
)
1523 return fdeja
.__init
__(field
,
1529 def inner_product(self
, x
, y
):
1530 # Since a+bi+cj+dk on the diagonal is represented as
1532 # a + bi +cj + dk = [ a b c d]
1537 # we'll quadruple-count the "a" entries if we take the trace of
1539 return _matrix_ip(x
,y
)/4
1542 class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1544 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
1545 with the usual inner product and jordan product ``x*y =
1546 (<x_bar,y_bar>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
1551 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1555 This multiplication table can be verified by hand::
1557 sage: J = JordanSpinEJA(4)
1558 sage: e0,e1,e2,e3 = J.gens()
1574 We can change the generator prefix::
1576 sage: JordanSpinEJA(2, prefix='B').gens()
1579 Our inner product satisfies the Jordan axiom::
1581 sage: set_random_seed()
1582 sage: n = ZZ.random_element(1,5)
1583 sage: J = JordanSpinEJA(n)
1584 sage: x = J.random_element()
1585 sage: y = J.random_element()
1586 sage: z = J.random_element()
1587 sage: (x*y).inner_product(z) == y.inner_product(x*z)
1591 def __init__(self
, n
, field
=QQ
, **kwargs
):
1592 V
= VectorSpace(field
, n
)
1593 mult_table
= [[V
.zero() for j
in range(n
)] for i
in range(n
)]
1603 z0
= x
.inner_product(y
)
1604 zbar
= y0
*xbar
+ x0
*ybar
1605 z
= V([z0
] + zbar
.list())
1606 mult_table
[i
][j
] = z
1608 # The rank of the spin algebra is two, unless we're in a
1609 # one-dimensional ambient space (because the rank is bounded by
1610 # the ambient dimension).
1611 fdeja
= super(JordanSpinEJA
, self
)
1612 return fdeja
.__init
__(field
, mult_table
, rank
=min(n
,2), **kwargs
)
1614 def inner_product(self
, x
, y
):
1616 Faster to reimplement than to use natural representations.
1620 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1624 Ensure that this is the usual inner product for the algebras
1627 sage: set_random_seed()
1628 sage: n = ZZ.random_element(1,5)
1629 sage: J = JordanSpinEJA(n)
1630 sage: x = J.random_element()
1631 sage: y = J.random_element()
1632 sage: X = x.natural_representation()
1633 sage: Y = y.natural_representation()
1634 sage: x.inner_product(y) == J.__class__.natural_inner_product(X,Y)
1638 return x
.to_vector().inner_product(y
.to_vector())