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
, AA
, QQbar
, 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
):
30 def _coerce_map_from_base_ring(self
):
32 Disable the map from the base ring into the algebra.
34 Performing a nonsense conversion like this automatically
35 is counterpedagogical. The fallback is to try the usual
36 element constructor, which should also fail.
40 sage: from mjo.eja.eja_algebra import random_eja
44 sage: set_random_seed()
45 sage: J = random_eja()
47 Traceback (most recent call last):
49 ValueError: not a naturally-represented algebra element
64 sage: from mjo.eja.eja_algebra import (JordanSpinEJA, random_eja)
68 By definition, Jordan multiplication commutes::
70 sage: set_random_seed()
71 sage: J = random_eja()
72 sage: x,y = J.random_elements(2)
78 The ``field`` we're given must be real::
80 sage: JordanSpinEJA(2,QQbar)
81 Traceback (most recent call last):
83 ValueError: field is not real
87 if not field
.is_subring(RR
):
88 # Note: this does return true for the real algebraic
89 # field, and any quadratic field where we've specified
91 raise ValueError('field is not real')
93 self
._natural
_basis
= natural_basis
96 category
= MagmaticAlgebras(field
).FiniteDimensional()
97 category
= category
.WithBasis().Unital()
99 fda
= super(FiniteDimensionalEuclideanJordanAlgebra
, self
)
101 range(len(mult_table
)),
104 self
.print_options(bracket
='')
106 # The multiplication table we're given is necessarily in terms
107 # of vectors, because we don't have an algebra yet for
108 # anything to be an element of. However, it's faster in the
109 # long run to have the multiplication table be in terms of
110 # algebra elements. We do this after calling the superclass
111 # constructor so that from_vector() knows what to do.
112 self
._multiplication
_table
= [
113 list(map(lambda x
: self
.from_vector(x
), ls
))
118 def _element_constructor_(self
, elt
):
120 Construct an element of this algebra from its natural
123 This gets called only after the parent element _call_ method
124 fails to find a coercion for the argument.
128 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
130 ....: RealSymmetricEJA)
134 The identity in `S^n` is converted to the identity in the EJA::
136 sage: J = RealSymmetricEJA(3)
137 sage: I = matrix.identity(QQ,3)
138 sage: J(I) == J.one()
141 This skew-symmetric matrix can't be represented in the EJA::
143 sage: J = RealSymmetricEJA(3)
144 sage: A = matrix(QQ,3, lambda i,j: i-j)
146 Traceback (most recent call last):
148 ArithmeticError: vector is not in free module
152 Ensure that we can convert any element of the two non-matrix
153 simple algebras (whose natural representations are their usual
154 vector representations) back and forth faithfully::
156 sage: set_random_seed()
157 sage: J = HadamardEJA.random_instance()
158 sage: x = J.random_element()
159 sage: J(x.to_vector().column()) == x
161 sage: J = JordanSpinEJA.random_instance()
162 sage: x = J.random_element()
163 sage: J(x.to_vector().column()) == x
167 msg
= "not a naturally-represented algebra element"
169 # The superclass implementation of random_element()
170 # needs to be able to coerce "0" into the algebra.
172 elif elt
in self
.base_ring():
173 # Ensure that no base ring -> algebra coercion is performed
174 # by this method. There's some stupidity in sage that would
175 # otherwise propagate to this method; for example, sage thinks
176 # that the integer 3 belongs to the space of 2-by-2 matrices.
177 raise ValueError(msg
)
179 natural_basis
= self
.natural_basis()
180 basis_space
= natural_basis
[0].matrix_space()
181 if elt
not in basis_space
:
182 raise ValueError(msg
)
184 # Thanks for nothing! Matrix spaces aren't vector spaces in
185 # Sage, so we have to figure out its natural-basis coordinates
186 # ourselves. We use the basis space's ring instead of the
187 # element's ring because the basis space might be an algebraic
188 # closure whereas the base ring of the 3-by-3 identity matrix
189 # could be QQ instead of QQbar.
190 V
= VectorSpace(basis_space
.base_ring(), elt
.nrows()*elt
.ncols())
191 W
= V
.span_of_basis( _mat2vec(s
) for s
in natural_basis
)
192 coords
= W
.coordinate_vector(_mat2vec(elt
))
193 return self
.from_vector(coords
)
196 def _max_test_case_size():
198 Return an integer "size" that is an upper bound on the size of
199 this algebra when it is used in a random test
200 case. Unfortunately, the term "size" is quite vague -- when
201 dealing with `R^n` under either the Hadamard or Jordan spin
202 product, the "size" refers to the dimension `n`. When dealing
203 with a matrix algebra (real symmetric or complex/quaternion
204 Hermitian), it refers to the size of the matrix, which is
205 far less than the dimension of the underlying vector space.
207 We default to five in this class, which is safe in `R^n`. The
208 matrix algebra subclasses (or any class where the "size" is
209 interpreted to be far less than the dimension) should override
210 with a smaller number.
216 Return a string representation of ``self``.
220 sage: from mjo.eja.eja_algebra import JordanSpinEJA
224 Ensure that it says what we think it says::
226 sage: JordanSpinEJA(2, field=AA)
227 Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
228 sage: JordanSpinEJA(3, field=RDF)
229 Euclidean Jordan algebra of dimension 3 over Real Double Field
232 fmt
= "Euclidean Jordan algebra of dimension {} over {}"
233 return fmt
.format(self
.dimension(), self
.base_ring())
235 def product_on_basis(self
, i
, j
):
236 return self
._multiplication
_table
[i
][j
]
238 def _a_regular_element(self
):
240 Guess a regular element. Needed to compute the basis for our
241 characteristic polynomial coefficients.
245 sage: from mjo.eja.eja_algebra import random_eja
249 Ensure that this hacky method succeeds for every algebra that we
250 know how to construct::
252 sage: set_random_seed()
253 sage: J = random_eja()
254 sage: J._a_regular_element().is_regular()
259 z
= self
.sum( (i
+1)*gs
[i
] for i
in range(len(gs
)) )
260 if not z
.is_regular():
261 raise ValueError("don't know a regular element")
266 def _charpoly_basis_space(self
):
268 Return the vector space spanned by the basis used in our
269 characteristic polynomial coefficients. This is used not only to
270 compute those coefficients, but also any time we need to
271 evaluate the coefficients (like when we compute the trace or
274 z
= self
._a
_regular
_element
()
275 # Don't use the parent vector space directly here in case this
276 # happens to be a subalgebra. In that case, we would be e.g.
277 # two-dimensional but span_of_basis() would expect three
279 V
= VectorSpace(self
.base_ring(), self
.vector_space().dimension())
280 basis
= [ (z
**k
).to_vector() for k
in range(self
.rank()) ]
281 V1
= V
.span_of_basis( basis
)
282 b
= (V1
.basis() + V1
.complement().basis())
283 return V
.span_of_basis(b
)
288 def _charpoly_coeff(self
, i
):
290 Return the coefficient polynomial "a_{i}" of this algebra's
291 general characteristic polynomial.
293 Having this be a separate cached method lets us compute and
294 store the trace/determinant (a_{r-1} and a_{0} respectively)
295 separate from the entire characteristic polynomial.
297 (A_of_x
, x
, xr
, detA
) = self
._charpoly
_matrix
_system
()
298 R
= A_of_x
.base_ring()
303 # Guaranteed by theory
306 # Danger: the in-place modification is done for performance
307 # reasons (reconstructing a matrix with huge polynomial
308 # entries is slow), but I don't know how cached_method works,
309 # so it's highly possible that we're modifying some global
310 # list variable by reference, here. In other words, you
311 # probably shouldn't call this method twice on the same
312 # algebra, at the same time, in two threads
313 Ai_orig
= A_of_x
.column(i
)
314 A_of_x
.set_column(i
,xr
)
315 numerator
= A_of_x
.det()
316 A_of_x
.set_column(i
,Ai_orig
)
318 # We're relying on the theory here to ensure that each a_i is
319 # indeed back in R, and the added negative signs are to make
320 # the whole charpoly expression sum to zero.
321 return R(-numerator
/detA
)
325 def _charpoly_matrix_system(self
):
327 Compute the matrix whose entries A_ij are polynomials in
328 X1,...,XN, the vector ``x`` of variables X1,...,XN, the vector
329 corresponding to `x^r` and the determinent of the matrix A =
330 [A_ij]. In other words, all of the fixed (cachable) data needed
331 to compute the coefficients of the characteristic polynomial.
336 # Turn my vector space into a module so that "vectors" can
337 # have multivatiate polynomial entries.
338 names
= tuple('X' + str(i
) for i
in range(1,n
+1))
339 R
= PolynomialRing(self
.base_ring(), names
)
341 # Using change_ring() on the parent's vector space doesn't work
342 # here because, in a subalgebra, that vector space has a basis
343 # and change_ring() tries to bring the basis along with it. And
344 # that doesn't work unless the new ring is a PID, which it usually
348 # Now let x = (X1,X2,...,Xn) be the vector whose entries are
352 # And figure out the "left multiplication by x" matrix in
355 monomial_matrices
= [ self
.monomial(i
).operator().matrix()
356 for i
in range(n
) ] # don't recompute these!
358 ek
= self
.monomial(k
).to_vector()
360 sum( x
[i
]*(monomial_matrices
[i
]*ek
)
361 for i
in range(n
) ) )
362 Lx
= matrix
.column(R
, lmbx_cols
)
364 # Now we can compute powers of x "symbolically"
365 x_powers
= [self
.one().to_vector(), x
]
366 for d
in range(2, r
+1):
367 x_powers
.append( Lx
*(x_powers
[-1]) )
369 idmat
= matrix
.identity(R
, n
)
371 W
= self
._charpoly
_basis
_space
()
372 W
= W
.change_ring(R
.fraction_field())
374 # Starting with the standard coordinates x = (X1,X2,...,Xn)
375 # and then converting the entries to W-coordinates allows us
376 # to pass in the standard coordinates to the charpoly and get
377 # back the right answer. Specifically, with x = (X1,X2,...,Xn),
380 # W.coordinates(x^2) eval'd at (standard z-coords)
384 # W-coords of (standard coords of x^2 eval'd at std-coords of z)
386 # We want the middle equivalent thing in our matrix, but use
387 # the first equivalent thing instead so that we can pass in
388 # standard coordinates.
389 x_powers
= [ W
.coordinate_vector(xp
) for xp
in x_powers
]
390 l2
= [idmat
.column(k
-1) for k
in range(r
+1, n
+1)]
391 A_of_x
= matrix
.column(R
, n
, (x_powers
[:r
] + l2
))
392 return (A_of_x
, x
, x_powers
[r
], A_of_x
.det())
396 def characteristic_polynomial(self
):
398 Return a characteristic polynomial that works for all elements
401 The resulting polynomial has `n+1` variables, where `n` is the
402 dimension of this algebra. The first `n` variables correspond to
403 the coordinates of an algebra element: when evaluated at the
404 coordinates of an algebra element with respect to a certain
405 basis, the result is a univariate polynomial (in the one
406 remaining variable ``t``), namely the characteristic polynomial
411 sage: from mjo.eja.eja_algebra import JordanSpinEJA, TrivialEJA
415 The characteristic polynomial in the spin algebra is given in
416 Alizadeh, Example 11.11::
418 sage: J = JordanSpinEJA(3)
419 sage: p = J.characteristic_polynomial(); p
420 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
421 sage: xvec = J.one().to_vector()
425 By definition, the characteristic polynomial is a monic
426 degree-zero polynomial in a rank-zero algebra. Note that
427 Cayley-Hamilton is indeed satisfied since the polynomial
428 ``1`` evaluates to the identity element of the algebra on
431 sage: J = TrivialEJA()
432 sage: J.characteristic_polynomial()
439 # The list of coefficient polynomials a_0, a_1, a_2, ..., a_n.
440 a
= [ self
._charpoly
_coeff
(i
) for i
in range(r
+1) ]
442 # We go to a bit of trouble here to reorder the
443 # indeterminates, so that it's easier to evaluate the
444 # characteristic polynomial at x's coordinates and get back
445 # something in terms of t, which is what we want.
447 S
= PolynomialRing(self
.base_ring(),'t')
449 S
= PolynomialRing(S
, R
.variable_names())
452 return sum( a
[k
]*(t
**k
) for k
in range(len(a
)) )
455 def inner_product(self
, x
, y
):
457 The inner product associated with this Euclidean Jordan algebra.
459 Defaults to the trace inner product, but can be overridden by
460 subclasses if they are sure that the necessary properties are
465 sage: from mjo.eja.eja_algebra import random_eja
469 Our inner product is "associative," which means the following for
470 a symmetric bilinear form::
472 sage: set_random_seed()
473 sage: J = random_eja()
474 sage: x,y,z = J.random_elements(3)
475 sage: (x*y).inner_product(z) == y.inner_product(x*z)
479 X
= x
.natural_representation()
480 Y
= y
.natural_representation()
481 return self
.natural_inner_product(X
,Y
)
484 def is_trivial(self
):
486 Return whether or not this algebra is trivial.
488 A trivial algebra contains only the zero element.
492 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
497 sage: J = ComplexHermitianEJA(3)
503 sage: J = TrivialEJA()
508 return self
.dimension() == 0
511 def multiplication_table(self
):
513 Return a visual representation of this algebra's multiplication
514 table (on basis elements).
518 sage: from mjo.eja.eja_algebra import JordanSpinEJA
522 sage: J = JordanSpinEJA(4)
523 sage: J.multiplication_table()
524 +----++----+----+----+----+
525 | * || e0 | e1 | e2 | e3 |
526 +====++====+====+====+====+
527 | e0 || e0 | e1 | e2 | e3 |
528 +----++----+----+----+----+
529 | e1 || e1 | e0 | 0 | 0 |
530 +----++----+----+----+----+
531 | e2 || e2 | 0 | e0 | 0 |
532 +----++----+----+----+----+
533 | e3 || e3 | 0 | 0 | e0 |
534 +----++----+----+----+----+
537 M
= list(self
._multiplication
_table
) # copy
538 for i
in range(len(M
)):
539 # M had better be "square"
540 M
[i
] = [self
.monomial(i
)] + M
[i
]
541 M
= [["*"] + list(self
.gens())] + M
542 return table(M
, header_row
=True, header_column
=True, frame
=True)
545 def natural_basis(self
):
547 Return a more-natural representation of this algebra's basis.
549 Every finite-dimensional Euclidean Jordan Algebra is a direct
550 sum of five simple algebras, four of which comprise Hermitian
551 matrices. This method returns the original "natural" basis
552 for our underlying vector space. (Typically, the natural basis
553 is used to construct the multiplication table in the first place.)
555 Note that this will always return a matrix. The standard basis
556 in `R^n` will be returned as `n`-by-`1` column matrices.
560 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
561 ....: RealSymmetricEJA)
565 sage: J = RealSymmetricEJA(2)
567 Finite family {0: e0, 1: e1, 2: e2}
568 sage: J.natural_basis()
570 [1 0] [ 0 0.7071067811865475?] [0 0]
571 [0 0], [0.7071067811865475? 0], [0 1]
576 sage: J = JordanSpinEJA(2)
578 Finite family {0: e0, 1: e1}
579 sage: J.natural_basis()
586 if self
._natural
_basis
is None:
587 M
= self
.natural_basis_space()
588 return tuple( M(b
.to_vector()) for b
in self
.basis() )
590 return self
._natural
_basis
593 def natural_basis_space(self
):
595 Return the matrix space in which this algebra's natural basis
598 if self
._natural
_basis
is None or len(self
._natural
_basis
) == 0:
599 return MatrixSpace(self
.base_ring(), self
.dimension(), 1)
601 return self
._natural
_basis
[0].matrix_space()
605 def natural_inner_product(X
,Y
):
607 Compute the inner product of two naturally-represented elements.
609 For example in the real symmetric matrix EJA, this will compute
610 the trace inner-product of two n-by-n symmetric matrices. The
611 default should work for the real cartesian product EJA, the
612 Jordan spin EJA, and the real symmetric matrices. The others
613 will have to be overridden.
615 return (X
.conjugate_transpose()*Y
).trace()
621 Return the unit element of this algebra.
625 sage: from mjo.eja.eja_algebra import (HadamardEJA,
630 sage: J = HadamardEJA(5)
632 e0 + e1 + e2 + e3 + e4
636 The identity element acts like the identity::
638 sage: set_random_seed()
639 sage: J = random_eja()
640 sage: x = J.random_element()
641 sage: J.one()*x == x and x*J.one() == x
644 The matrix of the unit element's operator is the identity::
646 sage: set_random_seed()
647 sage: J = random_eja()
648 sage: actual = J.one().operator().matrix()
649 sage: expected = matrix.identity(J.base_ring(), J.dimension())
650 sage: actual == expected
654 # We can brute-force compute the matrices of the operators
655 # that correspond to the basis elements of this algebra.
656 # If some linear combination of those basis elements is the
657 # algebra identity, then the same linear combination of
658 # their matrices has to be the identity matrix.
660 # Of course, matrices aren't vectors in sage, so we have to
661 # appeal to the "long vectors" isometry.
662 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
664 # Now we use basis linear algebra to find the coefficients,
665 # of the matrices-as-vectors-linear-combination, which should
666 # work for the original algebra basis too.
667 A
= matrix
.column(self
.base_ring(), oper_vecs
)
669 # We used the isometry on the left-hand side already, but we
670 # still need to do it for the right-hand side. Recall that we
671 # wanted something that summed to the identity matrix.
672 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
674 # Now if there's an identity element in the algebra, this should work.
675 coeffs
= A
.solve_right(b
)
676 return self
.linear_combination(zip(self
.gens(), coeffs
))
679 def peirce_decomposition(self
, c
):
681 The Peirce decomposition of this algebra relative to the
684 In the future, this can be extended to a complete system of
685 orthogonal idempotents.
689 - ``c`` -- an idempotent of this algebra.
693 A triple (J0, J5, J1) containing two subalgebras and one subspace
696 - ``J0`` -- the algebra on the eigenspace of ``c.operator()``
697 corresponding to the eigenvalue zero.
699 - ``J5`` -- the eigenspace (NOT a subalgebra) of ``c.operator()``
700 corresponding to the eigenvalue one-half.
702 - ``J1`` -- the algebra on the eigenspace of ``c.operator()``
703 corresponding to the eigenvalue one.
705 These are the only possible eigenspaces for that operator, and this
706 algebra is a direct sum of them. The spaces ``J0`` and ``J1`` are
707 orthogonal, and are subalgebras of this algebra with the appropriate
712 sage: from mjo.eja.eja_algebra import random_eja, RealSymmetricEJA
716 The canonical example comes from the symmetric matrices, which
717 decompose into diagonal and off-diagonal parts::
719 sage: J = RealSymmetricEJA(3)
720 sage: C = matrix(QQ, [ [1,0,0],
724 sage: J0,J5,J1 = J.peirce_decomposition(c)
726 Euclidean Jordan algebra of dimension 1...
728 Vector space of degree 6 and dimension 2...
730 Euclidean Jordan algebra of dimension 3...
734 Every algebra decomposes trivially with respect to its identity
737 sage: set_random_seed()
738 sage: J = random_eja()
739 sage: J0,J5,J1 = J.peirce_decomposition(J.one())
740 sage: J0.dimension() == 0 and J5.dimension() == 0
742 sage: J1.superalgebra() == J and J1.dimension() == J.dimension()
745 The identity elements in the two subalgebras are the
746 projections onto their respective subspaces of the
747 superalgebra's identity element::
749 sage: set_random_seed()
750 sage: J = random_eja()
751 sage: x = J.random_element()
752 sage: if not J.is_trivial():
753 ....: while x.is_nilpotent():
754 ....: x = J.random_element()
755 sage: c = x.subalgebra_idempotent()
756 sage: J0,J5,J1 = J.peirce_decomposition(c)
757 sage: J1(c) == J1.one()
759 sage: J0(J.one() - c) == J0.one()
763 if not c
.is_idempotent():
764 raise ValueError("element is not idempotent: %s" % c
)
766 # Default these to what they should be if they turn out to be
767 # trivial, because eigenspaces_left() won't return eigenvalues
768 # corresponding to trivial spaces (e.g. it returns only the
769 # eigenspace corresponding to lambda=1 if you take the
770 # decomposition relative to the identity element).
771 trivial
= FiniteDimensionalEuclideanJordanSubalgebra(self
, ())
772 J0
= trivial
# eigenvalue zero
773 J5
= VectorSpace(self
.base_ring(), 0) # eigenvalue one-half
774 J1
= trivial
# eigenvalue one
776 for (eigval
, eigspace
) in c
.operator().matrix().right_eigenspaces():
777 if eigval
== ~
(self
.base_ring()(2)):
780 gens
= tuple( self
.from_vector(b
) for b
in eigspace
.basis() )
781 subalg
= FiniteDimensionalEuclideanJordanSubalgebra(self
, gens
)
787 raise ValueError("unexpected eigenvalue: %s" % eigval
)
792 def random_elements(self
, count
):
794 Return ``count`` random elements as a tuple.
798 sage: from mjo.eja.eja_algebra import JordanSpinEJA
802 sage: J = JordanSpinEJA(3)
803 sage: x,y,z = J.random_elements(3)
804 sage: all( [ x in J, y in J, z in J ])
806 sage: len( J.random_elements(10) ) == 10
810 return tuple( self
.random_element() for idx
in range(count
) )
813 def random_instance(cls
, field
=AA
, **kwargs
):
815 Return a random instance of this type of algebra.
817 Beware, this will crash for "most instances" because the
818 constructor below looks wrong.
820 if cls
is TrivialEJA
:
821 # The TrivialEJA class doesn't take an "n" argument because
825 n
= ZZ
.random_element(cls
._max
_test
_case
_size
()) + 1
826 return cls(n
, field
, **kwargs
)
829 def _charpoly_coefficients(self
):
831 The `r` polynomial coefficients of the "characteristic polynomial
835 var_names
= [ "X" + str(z
) for z
in range(1,n
+1) ]
836 R
= PolynomialRing(self
.base_ring(), var_names
)
838 F
= R
.fraction_field()
841 # From a result in my book, these are the entries of the
842 # basis representation of L_x.
843 return sum( vars[k
]*self
.monomial(k
).operator().matrix()[i
,j
]
846 L_x
= matrix(F
, n
, n
, L_x_i_j
)
847 # Compute an extra power in case the rank is equal to
848 # the dimension (otherwise, we would stop at x^(r-1)).
849 x_powers
= [ (L_x
**k
)*self
.one().to_vector()
850 for k
in range(n
+1) ]
851 A
= matrix
.column(F
, x_powers
[:n
])
852 AE
= A
.extended_echelon_form()
858 # The theory says that only the first "r" coefficients are
859 # nonzero, and they actually live in the original polynomial
860 # ring and not the fraction field. We negate them because
861 # in the actual characteristic polynomial, they get moved
862 # to the other side where x^r lives.
863 return -A_rref
.solve_right(E
*b
).change_ring(R
)[:r
]
868 Return the rank of this EJA.
872 We first compute the polynomial "column matrices" `p_{k}` that
873 evaluate to `x^k` on the coordinates of `x`. Then, we begin
874 adding them to a matrix one at a time, and trying to solve the
875 system that makes `p_{0}`,`p_{1}`,..., `p_{s-1}` add up to
876 `p_{s}`. This will succeed only when `s` is the rank of the
877 algebra, as proven in a recent draft paper of mine.
881 sage: from mjo.eja.eja_algebra import (HadamardEJA,
883 ....: RealSymmetricEJA,
884 ....: ComplexHermitianEJA,
885 ....: QuaternionHermitianEJA,
890 The rank of the Jordan spin algebra is always two::
892 sage: JordanSpinEJA(2).rank()
894 sage: JordanSpinEJA(3).rank()
896 sage: JordanSpinEJA(4).rank()
899 The rank of the `n`-by-`n` Hermitian real, complex, or
900 quaternion matrices is `n`::
902 sage: RealSymmetricEJA(4).rank()
904 sage: ComplexHermitianEJA(3).rank()
906 sage: QuaternionHermitianEJA(2).rank()
911 Ensure that every EJA that we know how to construct has a
912 positive integer rank, unless the algebra is trivial in
913 which case its rank will be zero::
915 sage: set_random_seed()
916 sage: J = random_eja()
920 sage: r > 0 or (r == 0 and J.is_trivial())
923 Ensure that computing the rank actually works, since the ranks
924 of all simple algebras are known and will be cached by default::
926 sage: J = HadamardEJA(4)
927 sage: J.rank.clear_cache()
933 sage: J = JordanSpinEJA(4)
934 sage: J.rank.clear_cache()
940 sage: J = RealSymmetricEJA(3)
941 sage: J.rank.clear_cache()
947 sage: J = ComplexHermitianEJA(2)
948 sage: J.rank.clear_cache()
954 sage: J = QuaternionHermitianEJA(2)
955 sage: J.rank.clear_cache()
960 return len(self
._charpoly
_coefficients
())
963 def vector_space(self
):
965 Return the vector space that underlies this algebra.
969 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
973 sage: J = RealSymmetricEJA(2)
974 sage: J.vector_space()
975 Vector space of dimension 3 over...
978 return self
.zero().to_vector().parent().ambient_vector_space()
981 Element
= FiniteDimensionalEuclideanJordanAlgebraElement
984 class HadamardEJA(FiniteDimensionalEuclideanJordanAlgebra
):
986 Return the Euclidean Jordan Algebra corresponding to the set
987 `R^n` under the Hadamard product.
989 Note: this is nothing more than the Cartesian product of ``n``
990 copies of the spin algebra. Once Cartesian product algebras
991 are implemented, this can go.
995 sage: from mjo.eja.eja_algebra import HadamardEJA
999 This multiplication table can be verified by hand::
1001 sage: J = HadamardEJA(3)
1002 sage: e0,e1,e2 = J.gens()
1018 We can change the generator prefix::
1020 sage: HadamardEJA(3, prefix='r').gens()
1024 def __init__(self
, n
, field
=AA
, **kwargs
):
1025 V
= VectorSpace(field
, n
)
1026 mult_table
= [ [ V
.gen(i
)*(i
== j
) for j
in range(n
) ]
1029 fdeja
= super(HadamardEJA
, self
)
1030 fdeja
.__init
__(field
, mult_table
, **kwargs
)
1031 self
.rank
.set_cache(n
)
1033 def inner_product(self
, x
, y
):
1035 Faster to reimplement than to use natural representations.
1039 sage: from mjo.eja.eja_algebra import HadamardEJA
1043 Ensure that this is the usual inner product for the algebras
1046 sage: set_random_seed()
1047 sage: J = HadamardEJA.random_instance()
1048 sage: x,y = J.random_elements(2)
1049 sage: X = x.natural_representation()
1050 sage: Y = y.natural_representation()
1051 sage: x.inner_product(y) == J.natural_inner_product(X,Y)
1055 return x
.to_vector().inner_product(y
.to_vector())
1058 def random_eja(field
=AA
, nontrivial
=False):
1060 Return a "random" finite-dimensional Euclidean Jordan Algebra.
1064 sage: from mjo.eja.eja_algebra import random_eja
1069 Euclidean Jordan algebra of dimension...
1072 eja_classes
= [HadamardEJA
,
1075 ComplexHermitianEJA
,
1076 QuaternionHermitianEJA
]
1078 eja_classes
.append(TrivialEJA
)
1079 classname
= choice(eja_classes
)
1080 return classname
.random_instance(field
=field
)
1087 class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra
):
1089 def _max_test_case_size():
1090 # Play it safe, since this will be squared and the underlying
1091 # field can have dimension 4 (quaternions) too.
1094 def __init__(self
, field
, basis
, normalize_basis
=True, **kwargs
):
1096 Compared to the superclass constructor, we take a basis instead of
1097 a multiplication table because the latter can be computed in terms
1098 of the former when the product is known (like it is here).
1100 # Used in this class's fast _charpoly_coeff() override.
1101 self
._basis
_normalizers
= None
1103 # We're going to loop through this a few times, so now's a good
1104 # time to ensure that it isn't a generator expression.
1105 basis
= tuple(basis
)
1107 if len(basis
) > 1 and normalize_basis
:
1108 # We'll need sqrt(2) to normalize the basis, and this
1109 # winds up in the multiplication table, so the whole
1110 # algebra needs to be over the field extension.
1111 R
= PolynomialRing(field
, 'z')
1114 if p
.is_irreducible():
1115 field
= field
.extension(p
, 'sqrt2', embedding
=RLF(2).sqrt())
1116 basis
= tuple( s
.change_ring(field
) for s
in basis
)
1117 self
._basis
_normalizers
= tuple(
1118 ~
(self
.natural_inner_product(s
,s
).sqrt()) for s
in basis
)
1119 basis
= tuple(s
*c
for (s
,c
) in zip(basis
,self
._basis
_normalizers
))
1121 Qs
= self
.multiplication_table_from_matrix_basis(basis
)
1123 fdeja
= super(MatrixEuclideanJordanAlgebra
, self
)
1124 fdeja
.__init
__(field
, Qs
, natural_basis
=basis
, **kwargs
)
1131 Override the parent method with something that tries to compute
1132 over a faster (non-extension) field.
1134 if self
._basis
_normalizers
is None:
1135 # We didn't normalize, so assume that the basis we started
1136 # with had entries in a nice field.
1137 return super(MatrixEuclideanJordanAlgebra
, self
).rank()
1139 basis
= ( (b
/n
) for (b
,n
) in zip(self
.natural_basis(),
1140 self
._basis
_normalizers
) )
1142 # Do this over the rationals and convert back at the end.
1143 # Only works because we know the entries of the basis are
1145 J
= MatrixEuclideanJordanAlgebra(QQ
,
1147 normalize_basis
=False)
1151 def _charpoly_coeff(self
, i
):
1153 Override the parent method with something that tries to compute
1154 over a faster (non-extension) field.
1156 if self
._basis
_normalizers
is None:
1157 # We didn't normalize, so assume that the basis we started
1158 # with had entries in a nice field.
1159 return super(MatrixEuclideanJordanAlgebra
, self
)._charpoly
_coeff
(i
)
1161 basis
= ( (b
/n
) for (b
,n
) in zip(self
.natural_basis(),
1162 self
._basis
_normalizers
) )
1164 # Do this over the rationals and convert back at the end.
1165 J
= MatrixEuclideanJordanAlgebra(QQ
,
1167 normalize_basis
=False)
1168 (_
,x
,_
,_
) = J
._charpoly
_matrix
_system
()
1169 p
= J
._charpoly
_coeff
(i
)
1170 # p might be missing some vars, have to substitute "optionally"
1171 pairs
= zip(x
.base_ring().gens(), self
._basis
_normalizers
)
1172 substitutions
= { v: v*c for (v,c) in pairs }
1173 result
= p
.subs(substitutions
)
1175 # The result of "subs" can be either a coefficient-ring
1176 # element or a polynomial. Gotta handle both cases.
1178 return self
.base_ring()(result
)
1180 return result
.change_ring(self
.base_ring())
1184 def multiplication_table_from_matrix_basis(basis
):
1186 At least three of the five simple Euclidean Jordan algebras have the
1187 symmetric multiplication (A,B) |-> (AB + BA)/2, where the
1188 multiplication on the right is matrix multiplication. Given a basis
1189 for the underlying matrix space, this function returns a
1190 multiplication table (obtained by looping through the basis
1191 elements) for an algebra of those matrices.
1193 # In S^2, for example, we nominally have four coordinates even
1194 # though the space is of dimension three only. The vector space V
1195 # is supposed to hold the entire long vector, and the subspace W
1196 # of V will be spanned by the vectors that arise from symmetric
1197 # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
1198 field
= basis
[0].base_ring()
1199 dimension
= basis
[0].nrows()
1201 V
= VectorSpace(field
, dimension
**2)
1202 W
= V
.span_of_basis( _mat2vec(s
) for s
in basis
)
1204 mult_table
= [[W
.zero() for j
in range(n
)] for i
in range(n
)]
1207 mat_entry
= (basis
[i
]*basis
[j
] + basis
[j
]*basis
[i
])/2
1208 mult_table
[i
][j
] = W
.coordinate_vector(_mat2vec(mat_entry
))
1216 Embed the matrix ``M`` into a space of real matrices.
1218 The matrix ``M`` can have entries in any field at the moment:
1219 the real numbers, complex numbers, or quaternions. And although
1220 they are not a field, we can probably support octonions at some
1221 point, too. This function returns a real matrix that "acts like"
1222 the original with respect to matrix multiplication; i.e.
1224 real_embed(M*N) = real_embed(M)*real_embed(N)
1227 raise NotImplementedError
1231 def real_unembed(M
):
1233 The inverse of :meth:`real_embed`.
1235 raise NotImplementedError
1239 def natural_inner_product(cls
,X
,Y
):
1240 Xu
= cls
.real_unembed(X
)
1241 Yu
= cls
.real_unembed(Y
)
1242 tr
= (Xu
*Yu
).trace()
1245 # It's real already.
1248 # Otherwise, try the thing that works for complex numbers; and
1249 # if that doesn't work, the thing that works for quaternions.
1251 return tr
.vector()[0] # real part, imag part is index 1
1252 except AttributeError:
1253 # A quaternions doesn't have a vector() method, but does
1254 # have coefficient_tuple() method that returns the
1255 # coefficients of 1, i, j, and k -- in that order.
1256 return tr
.coefficient_tuple()[0]
1259 class RealMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1263 The identity function, for embedding real matrices into real
1269 def real_unembed(M
):
1271 The identity function, for unembedding real matrices from real
1277 class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra
):
1279 The rank-n simple EJA consisting of real symmetric n-by-n
1280 matrices, the usual symmetric Jordan product, and the trace inner
1281 product. It has dimension `(n^2 + n)/2` over the reals.
1285 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1289 sage: J = RealSymmetricEJA(2)
1290 sage: e0, e1, e2 = J.gens()
1298 In theory, our "field" can be any subfield of the reals::
1300 sage: RealSymmetricEJA(2, RDF)
1301 Euclidean Jordan algebra of dimension 3 over Real Double Field
1302 sage: RealSymmetricEJA(2, RR)
1303 Euclidean Jordan algebra of dimension 3 over Real Field with
1304 53 bits of precision
1308 The dimension of this algebra is `(n^2 + n) / 2`::
1310 sage: set_random_seed()
1311 sage: n_max = RealSymmetricEJA._max_test_case_size()
1312 sage: n = ZZ.random_element(1, n_max)
1313 sage: J = RealSymmetricEJA(n)
1314 sage: J.dimension() == (n^2 + n)/2
1317 The Jordan multiplication is what we think it is::
1319 sage: set_random_seed()
1320 sage: J = RealSymmetricEJA.random_instance()
1321 sage: x,y = J.random_elements(2)
1322 sage: actual = (x*y).natural_representation()
1323 sage: X = x.natural_representation()
1324 sage: Y = y.natural_representation()
1325 sage: expected = (X*Y + Y*X)/2
1326 sage: actual == expected
1328 sage: J(expected) == x*y
1331 We can change the generator prefix::
1333 sage: RealSymmetricEJA(3, prefix='q').gens()
1334 (q0, q1, q2, q3, q4, q5)
1336 Our natural basis is normalized with respect to the natural inner
1337 product unless we specify otherwise::
1339 sage: set_random_seed()
1340 sage: J = RealSymmetricEJA.random_instance()
1341 sage: all( b.norm() == 1 for b in J.gens() )
1344 Since our natural basis is normalized with respect to the natural
1345 inner product, and since we know that this algebra is an EJA, any
1346 left-multiplication operator's matrix will be symmetric because
1347 natural->EJA basis representation is an isometry and within the EJA
1348 the operator is self-adjoint by the Jordan axiom::
1350 sage: set_random_seed()
1351 sage: x = RealSymmetricEJA.random_instance().random_element()
1352 sage: x.operator().matrix().is_symmetric()
1357 def _denormalized_basis(cls
, n
, field
):
1359 Return a basis for the space of real symmetric n-by-n matrices.
1363 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1367 sage: set_random_seed()
1368 sage: n = ZZ.random_element(1,5)
1369 sage: B = RealSymmetricEJA._denormalized_basis(n,QQ)
1370 sage: all( M.is_symmetric() for M in B)
1374 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1378 for j
in range(i
+1):
1379 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
1383 Sij
= Eij
+ Eij
.transpose()
1389 def _max_test_case_size():
1390 return 4 # Dimension 10
1393 def __init__(self
, n
, field
=AA
, **kwargs
):
1394 basis
= self
._denormalized
_basis
(n
, field
)
1395 super(RealSymmetricEJA
, self
).__init
__(field
, basis
, **kwargs
)
1396 self
.rank
.set_cache(n
)
1399 class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1403 Embed the n-by-n complex matrix ``M`` into the space of real
1404 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1405 bi` to the block matrix ``[[a,b],[-b,a]]``.
1409 sage: from mjo.eja.eja_algebra import \
1410 ....: ComplexMatrixEuclideanJordanAlgebra
1414 sage: F = QuadraticField(-1, 'I')
1415 sage: x1 = F(4 - 2*i)
1416 sage: x2 = F(1 + 2*i)
1419 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1420 sage: ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1429 Embedding is a homomorphism (isomorphism, in fact)::
1431 sage: set_random_seed()
1432 sage: n_max = ComplexMatrixEuclideanJordanAlgebra._max_test_case_size()
1433 sage: n = ZZ.random_element(n_max)
1434 sage: F = QuadraticField(-1, 'I')
1435 sage: X = random_matrix(F, n)
1436 sage: Y = random_matrix(F, n)
1437 sage: Xe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X)
1438 sage: Ye = ComplexMatrixEuclideanJordanAlgebra.real_embed(Y)
1439 sage: XYe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1446 raise ValueError("the matrix 'M' must be square")
1448 # We don't need any adjoined elements...
1449 field
= M
.base_ring().base_ring()
1453 a
= z
.list()[0] # real part, I guess
1454 b
= z
.list()[1] # imag part, I guess
1455 blocks
.append(matrix(field
, 2, [[a
,b
],[-b
,a
]]))
1457 return matrix
.block(field
, n
, blocks
)
1461 def real_unembed(M
):
1463 The inverse of _embed_complex_matrix().
1467 sage: from mjo.eja.eja_algebra import \
1468 ....: ComplexMatrixEuclideanJordanAlgebra
1472 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1473 ....: [-2, 1, -4, 3],
1474 ....: [ 9, 10, 11, 12],
1475 ....: [-10, 9, -12, 11] ])
1476 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(A)
1478 [ 10*I + 9 12*I + 11]
1482 Unembedding is the inverse of embedding::
1484 sage: set_random_seed()
1485 sage: F = QuadraticField(-1, 'I')
1486 sage: M = random_matrix(F, 3)
1487 sage: Me = ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1488 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1494 raise ValueError("the matrix 'M' must be square")
1495 if not n
.mod(2).is_zero():
1496 raise ValueError("the matrix 'M' must be a complex embedding")
1498 # If "M" was normalized, its base ring might have roots
1499 # adjoined and they can stick around after unembedding.
1500 field
= M
.base_ring()
1501 R
= PolynomialRing(field
, 'z')
1504 # Sage doesn't know how to embed AA into QQbar, i.e. how
1505 # to adjoin sqrt(-1) to AA.
1508 F
= field
.extension(z
**2 + 1, 'I', embedding
=CLF(-1).sqrt())
1511 # Go top-left to bottom-right (reading order), converting every
1512 # 2-by-2 block we see to a single complex element.
1514 for k
in range(n
/2):
1515 for j
in range(n
/2):
1516 submat
= M
[2*k
:2*k
+2,2*j
:2*j
+2]
1517 if submat
[0,0] != submat
[1,1]:
1518 raise ValueError('bad on-diagonal submatrix')
1519 if submat
[0,1] != -submat
[1,0]:
1520 raise ValueError('bad off-diagonal submatrix')
1521 z
= submat
[0,0] + submat
[0,1]*i
1524 return matrix(F
, n
/2, elements
)
1528 def natural_inner_product(cls
,X
,Y
):
1530 Compute a natural inner product in this algebra directly from
1535 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1539 This gives the same answer as the slow, default method implemented
1540 in :class:`MatrixEuclideanJordanAlgebra`::
1542 sage: set_random_seed()
1543 sage: J = ComplexHermitianEJA.random_instance()
1544 sage: x,y = J.random_elements(2)
1545 sage: Xe = x.natural_representation()
1546 sage: Ye = y.natural_representation()
1547 sage: X = ComplexHermitianEJA.real_unembed(Xe)
1548 sage: Y = ComplexHermitianEJA.real_unembed(Ye)
1549 sage: expected = (X*Y).trace().real()
1550 sage: actual = ComplexHermitianEJA.natural_inner_product(Xe,Ye)
1551 sage: actual == expected
1555 return RealMatrixEuclideanJordanAlgebra
.natural_inner_product(X
,Y
)/2
1558 class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra
):
1560 The rank-n simple EJA consisting of complex Hermitian n-by-n
1561 matrices over the real numbers, the usual symmetric Jordan product,
1562 and the real-part-of-trace inner product. It has dimension `n^2` over
1567 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1571 In theory, our "field" can be any subfield of the reals::
1573 sage: ComplexHermitianEJA(2, RDF)
1574 Euclidean Jordan algebra of dimension 4 over Real Double Field
1575 sage: ComplexHermitianEJA(2, RR)
1576 Euclidean Jordan algebra of dimension 4 over Real Field with
1577 53 bits of precision
1581 The dimension of this algebra is `n^2`::
1583 sage: set_random_seed()
1584 sage: n_max = ComplexHermitianEJA._max_test_case_size()
1585 sage: n = ZZ.random_element(1, n_max)
1586 sage: J = ComplexHermitianEJA(n)
1587 sage: J.dimension() == n^2
1590 The Jordan multiplication is what we think it is::
1592 sage: set_random_seed()
1593 sage: J = ComplexHermitianEJA.random_instance()
1594 sage: x,y = J.random_elements(2)
1595 sage: actual = (x*y).natural_representation()
1596 sage: X = x.natural_representation()
1597 sage: Y = y.natural_representation()
1598 sage: expected = (X*Y + Y*X)/2
1599 sage: actual == expected
1601 sage: J(expected) == x*y
1604 We can change the generator prefix::
1606 sage: ComplexHermitianEJA(2, prefix='z').gens()
1609 Our natural basis is normalized with respect to the natural inner
1610 product unless we specify otherwise::
1612 sage: set_random_seed()
1613 sage: J = ComplexHermitianEJA.random_instance()
1614 sage: all( b.norm() == 1 for b in J.gens() )
1617 Since our natural basis is normalized with respect to the natural
1618 inner product, and since we know that this algebra is an EJA, any
1619 left-multiplication operator's matrix will be symmetric because
1620 natural->EJA basis representation is an isometry and within the EJA
1621 the operator is self-adjoint by the Jordan axiom::
1623 sage: set_random_seed()
1624 sage: x = ComplexHermitianEJA.random_instance().random_element()
1625 sage: x.operator().matrix().is_symmetric()
1631 def _denormalized_basis(cls
, n
, field
):
1633 Returns a basis for the space of complex Hermitian n-by-n matrices.
1635 Why do we embed these? Basically, because all of numerical linear
1636 algebra assumes that you're working with vectors consisting of `n`
1637 entries from a field and scalars from the same field. There's no way
1638 to tell SageMath that (for example) the vectors contain complex
1639 numbers, while the scalar field is real.
1643 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1647 sage: set_random_seed()
1648 sage: n = ZZ.random_element(1,5)
1649 sage: field = QuadraticField(2, 'sqrt2')
1650 sage: B = ComplexHermitianEJA._denormalized_basis(n, field)
1651 sage: all( M.is_symmetric() for M in B)
1655 R
= PolynomialRing(field
, 'z')
1657 F
= field
.extension(z
**2 + 1, 'I')
1660 # This is like the symmetric case, but we need to be careful:
1662 # * We want conjugate-symmetry, not just symmetry.
1663 # * The diagonal will (as a result) be real.
1667 for j
in range(i
+1):
1668 Eij
= matrix(F
, n
, lambda k
,l
: k
==i
and l
==j
)
1670 Sij
= cls
.real_embed(Eij
)
1673 # The second one has a minus because it's conjugated.
1674 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
1676 Sij_imag
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
1679 # Since we embedded these, we can drop back to the "field" that we
1680 # started with instead of the complex extension "F".
1681 return ( s
.change_ring(field
) for s
in S
)
1684 def __init__(self
, n
, field
=AA
, **kwargs
):
1685 basis
= self
._denormalized
_basis
(n
,field
)
1686 super(ComplexHermitianEJA
,self
).__init
__(field
, basis
, **kwargs
)
1687 self
.rank
.set_cache(n
)
1690 class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1694 Embed the n-by-n quaternion matrix ``M`` into the space of real
1695 matrices of size 4n-by-4n by first sending each quaternion entry `z
1696 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
1697 c+di],[-c + di, a-bi]]`, and then embedding those into a real
1702 sage: from mjo.eja.eja_algebra import \
1703 ....: QuaternionMatrixEuclideanJordanAlgebra
1707 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1708 sage: i,j,k = Q.gens()
1709 sage: x = 1 + 2*i + 3*j + 4*k
1710 sage: M = matrix(Q, 1, [[x]])
1711 sage: QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1717 Embedding is a homomorphism (isomorphism, in fact)::
1719 sage: set_random_seed()
1720 sage: n_max = QuaternionMatrixEuclideanJordanAlgebra._max_test_case_size()
1721 sage: n = ZZ.random_element(n_max)
1722 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1723 sage: X = random_matrix(Q, n)
1724 sage: Y = random_matrix(Q, n)
1725 sage: Xe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X)
1726 sage: Ye = QuaternionMatrixEuclideanJordanAlgebra.real_embed(Y)
1727 sage: XYe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1732 quaternions
= M
.base_ring()
1735 raise ValueError("the matrix 'M' must be square")
1737 F
= QuadraticField(-1, 'I')
1742 t
= z
.coefficient_tuple()
1747 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
1748 [-c
+ d
*i
, a
- b
*i
]])
1749 realM
= ComplexMatrixEuclideanJordanAlgebra
.real_embed(cplxM
)
1750 blocks
.append(realM
)
1752 # We should have real entries by now, so use the realest field
1753 # we've got for the return value.
1754 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
1759 def real_unembed(M
):
1761 The inverse of _embed_quaternion_matrix().
1765 sage: from mjo.eja.eja_algebra import \
1766 ....: QuaternionMatrixEuclideanJordanAlgebra
1770 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
1771 ....: [-2, 1, -4, 3],
1772 ....: [-3, 4, 1, -2],
1773 ....: [-4, -3, 2, 1]])
1774 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(M)
1775 [1 + 2*i + 3*j + 4*k]
1779 Unembedding is the inverse of embedding::
1781 sage: set_random_seed()
1782 sage: Q = QuaternionAlgebra(QQ, -1, -1)
1783 sage: M = random_matrix(Q, 3)
1784 sage: Me = QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1785 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1791 raise ValueError("the matrix 'M' must be square")
1792 if not n
.mod(4).is_zero():
1793 raise ValueError("the matrix 'M' must be a quaternion embedding")
1795 # Use the base ring of the matrix to ensure that its entries can be
1796 # multiplied by elements of the quaternion algebra.
1797 field
= M
.base_ring()
1798 Q
= QuaternionAlgebra(field
,-1,-1)
1801 # Go top-left to bottom-right (reading order), converting every
1802 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
1805 for l
in range(n
/4):
1806 for m
in range(n
/4):
1807 submat
= ComplexMatrixEuclideanJordanAlgebra
.real_unembed(
1808 M
[4*l
:4*l
+4,4*m
:4*m
+4] )
1809 if submat
[0,0] != submat
[1,1].conjugate():
1810 raise ValueError('bad on-diagonal submatrix')
1811 if submat
[0,1] != -submat
[1,0].conjugate():
1812 raise ValueError('bad off-diagonal submatrix')
1813 z
= submat
[0,0].real()
1814 z
+= submat
[0,0].imag()*i
1815 z
+= submat
[0,1].real()*j
1816 z
+= submat
[0,1].imag()*k
1819 return matrix(Q
, n
/4, elements
)
1823 def natural_inner_product(cls
,X
,Y
):
1825 Compute a natural inner product in this algebra directly from
1830 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1834 This gives the same answer as the slow, default method implemented
1835 in :class:`MatrixEuclideanJordanAlgebra`::
1837 sage: set_random_seed()
1838 sage: J = QuaternionHermitianEJA.random_instance()
1839 sage: x,y = J.random_elements(2)
1840 sage: Xe = x.natural_representation()
1841 sage: Ye = y.natural_representation()
1842 sage: X = QuaternionHermitianEJA.real_unembed(Xe)
1843 sage: Y = QuaternionHermitianEJA.real_unembed(Ye)
1844 sage: expected = (X*Y).trace().coefficient_tuple()[0]
1845 sage: actual = QuaternionHermitianEJA.natural_inner_product(Xe,Ye)
1846 sage: actual == expected
1850 return RealMatrixEuclideanJordanAlgebra
.natural_inner_product(X
,Y
)/4
1853 class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra
):
1855 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
1856 matrices, the usual symmetric Jordan product, and the
1857 real-part-of-trace inner product. It has dimension `2n^2 - n` over
1862 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1866 In theory, our "field" can be any subfield of the reals::
1868 sage: QuaternionHermitianEJA(2, RDF)
1869 Euclidean Jordan algebra of dimension 6 over Real Double Field
1870 sage: QuaternionHermitianEJA(2, RR)
1871 Euclidean Jordan algebra of dimension 6 over Real Field with
1872 53 bits of precision
1876 The dimension of this algebra is `2*n^2 - n`::
1878 sage: set_random_seed()
1879 sage: n_max = QuaternionHermitianEJA._max_test_case_size()
1880 sage: n = ZZ.random_element(1, n_max)
1881 sage: J = QuaternionHermitianEJA(n)
1882 sage: J.dimension() == 2*(n^2) - n
1885 The Jordan multiplication is what we think it is::
1887 sage: set_random_seed()
1888 sage: J = QuaternionHermitianEJA.random_instance()
1889 sage: x,y = J.random_elements(2)
1890 sage: actual = (x*y).natural_representation()
1891 sage: X = x.natural_representation()
1892 sage: Y = y.natural_representation()
1893 sage: expected = (X*Y + Y*X)/2
1894 sage: actual == expected
1896 sage: J(expected) == x*y
1899 We can change the generator prefix::
1901 sage: QuaternionHermitianEJA(2, prefix='a').gens()
1902 (a0, a1, a2, a3, a4, a5)
1904 Our natural basis is normalized with respect to the natural inner
1905 product unless we specify otherwise::
1907 sage: set_random_seed()
1908 sage: J = QuaternionHermitianEJA.random_instance()
1909 sage: all( b.norm() == 1 for b in J.gens() )
1912 Since our natural basis is normalized with respect to the natural
1913 inner product, and since we know that this algebra is an EJA, any
1914 left-multiplication operator's matrix will be symmetric because
1915 natural->EJA basis representation is an isometry and within the EJA
1916 the operator is self-adjoint by the Jordan axiom::
1918 sage: set_random_seed()
1919 sage: x = QuaternionHermitianEJA.random_instance().random_element()
1920 sage: x.operator().matrix().is_symmetric()
1925 def _denormalized_basis(cls
, n
, field
):
1927 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
1929 Why do we embed these? Basically, because all of numerical
1930 linear algebra assumes that you're working with vectors consisting
1931 of `n` entries from a field and scalars from the same field. There's
1932 no way to tell SageMath that (for example) the vectors contain
1933 complex numbers, while the scalar field is real.
1937 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1941 sage: set_random_seed()
1942 sage: n = ZZ.random_element(1,5)
1943 sage: B = QuaternionHermitianEJA._denormalized_basis(n,QQ)
1944 sage: all( M.is_symmetric() for M in B )
1948 Q
= QuaternionAlgebra(QQ
,-1,-1)
1951 # This is like the symmetric case, but we need to be careful:
1953 # * We want conjugate-symmetry, not just symmetry.
1954 # * The diagonal will (as a result) be real.
1958 for j
in range(i
+1):
1959 Eij
= matrix(Q
, n
, lambda k
,l
: k
==i
and l
==j
)
1961 Sij
= cls
.real_embed(Eij
)
1964 # The second, third, and fourth ones have a minus
1965 # because they're conjugated.
1966 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
1968 Sij_I
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
1970 Sij_J
= cls
.real_embed(J
*Eij
- J
*Eij
.transpose())
1972 Sij_K
= cls
.real_embed(K
*Eij
- K
*Eij
.transpose())
1975 # Since we embedded these, we can drop back to the "field" that we
1976 # started with instead of the quaternion algebra "Q".
1977 return ( s
.change_ring(field
) for s
in S
)
1980 def __init__(self
, n
, field
=AA
, **kwargs
):
1981 basis
= self
._denormalized
_basis
(n
,field
)
1982 super(QuaternionHermitianEJA
,self
).__init
__(field
, basis
, **kwargs
)
1983 self
.rank
.set_cache(n
)
1986 class BilinearFormEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1988 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
1989 with the half-trace inner product and jordan product ``x*y =
1990 (x0*y0 + <B*x_bar,y_bar>, x0*y_bar + y0*x_bar)`` where ``B`` is a
1991 symmetric positive-definite "bilinear form" matrix. It has
1992 dimension `n` over the reals, and reduces to the ``JordanSpinEJA``
1993 when ``B`` is the identity matrix of order ``n-1``.
1997 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
1998 ....: JordanSpinEJA)
2002 When no bilinear form is specified, the identity matrix is used,
2003 and the resulting algebra is the Jordan spin algebra::
2005 sage: J0 = BilinearFormEJA(3)
2006 sage: J1 = JordanSpinEJA(3)
2007 sage: J0.multiplication_table() == J0.multiplication_table()
2012 We can create a zero-dimensional algebra::
2014 sage: J = BilinearFormEJA(0)
2018 We can check the multiplication condition given in the Jordan, von
2019 Neumann, and Wigner paper (and also discussed on my "On the
2020 symmetry..." paper). Note that this relies heavily on the standard
2021 choice of basis, as does anything utilizing the bilinear form matrix::
2023 sage: set_random_seed()
2024 sage: n = ZZ.random_element(5)
2025 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2026 sage: B = M.transpose()*M
2027 sage: J = BilinearFormEJA(n, B=B)
2028 sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
2029 sage: V = J.vector_space()
2030 sage: sis = [ J.from_vector(V([0] + (M.inverse()*ei).list()))
2031 ....: for ei in eis ]
2032 sage: actual = [ sis[i]*sis[j]
2033 ....: for i in range(n-1)
2034 ....: for j in range(n-1) ]
2035 sage: expected = [ J.one() if i == j else J.zero()
2036 ....: for i in range(n-1)
2037 ....: for j in range(n-1) ]
2038 sage: actual == expected
2041 def __init__(self
, n
, field
=AA
, B
=None, **kwargs
):
2043 self
._B
= matrix
.identity(field
, max(0,n
-1))
2047 V
= VectorSpace(field
, n
)
2048 mult_table
= [[V
.zero() for j
in range(n
)] for i
in range(n
)]
2057 z0
= x0
*y0
+ (self
._B
*xbar
).inner_product(ybar
)
2058 zbar
= y0
*xbar
+ x0
*ybar
2059 z
= V([z0
] + zbar
.list())
2060 mult_table
[i
][j
] = z
2062 # The rank of this algebra is two, unless we're in a
2063 # one-dimensional ambient space (because the rank is bounded
2064 # by the ambient dimension).
2065 fdeja
= super(BilinearFormEJA
, self
)
2066 fdeja
.__init
__(field
, mult_table
, **kwargs
)
2067 self
.rank
.set_cache(min(n
,2))
2069 def inner_product(self
, x
, y
):
2071 Half of the trace inner product.
2073 This is defined so that the special case of the Jordan spin
2074 algebra gets the usual inner product.
2078 sage: from mjo.eja.eja_algebra import BilinearFormEJA
2082 Ensure that this is one-half of the trace inner-product when
2083 the algebra isn't just the reals (when ``n`` isn't one). This
2084 is in Faraut and Koranyi, and also my "On the symmetry..."
2087 sage: set_random_seed()
2088 sage: n = ZZ.random_element(2,5)
2089 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2090 sage: B = M.transpose()*M
2091 sage: J = BilinearFormEJA(n, B=B)
2092 sage: x = J.random_element()
2093 sage: y = J.random_element()
2094 sage: x.inner_product(y) == (x*y).trace()/2
2098 xvec
= x
.to_vector()
2100 yvec
= y
.to_vector()
2102 return x
[0]*y
[0] + (self
._B
*xbar
).inner_product(ybar
)
2105 class JordanSpinEJA(BilinearFormEJA
):
2107 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2108 with the usual inner product and jordan product ``x*y =
2109 (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2114 sage: from mjo.eja.eja_algebra import JordanSpinEJA
2118 This multiplication table can be verified by hand::
2120 sage: J = JordanSpinEJA(4)
2121 sage: e0,e1,e2,e3 = J.gens()
2137 We can change the generator prefix::
2139 sage: JordanSpinEJA(2, prefix='B').gens()
2144 Ensure that we have the usual inner product on `R^n`::
2146 sage: set_random_seed()
2147 sage: J = JordanSpinEJA.random_instance()
2148 sage: x,y = J.random_elements(2)
2149 sage: X = x.natural_representation()
2150 sage: Y = y.natural_representation()
2151 sage: x.inner_product(y) == J.natural_inner_product(X,Y)
2155 def __init__(self
, n
, field
=AA
, **kwargs
):
2156 # This is a special case of the BilinearFormEJA with the identity
2157 # matrix as its bilinear form.
2158 return super(JordanSpinEJA
, self
).__init
__(n
, field
, **kwargs
)
2161 class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra
):
2163 The trivial Euclidean Jordan algebra consisting of only a zero element.
2167 sage: from mjo.eja.eja_algebra import TrivialEJA
2171 sage: J = TrivialEJA()
2178 sage: 7*J.one()*12*J.one()
2180 sage: J.one().inner_product(J.one())
2182 sage: J.one().norm()
2184 sage: J.one().subalgebra_generated_by()
2185 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2190 def __init__(self
, field
=AA
, **kwargs
):
2192 fdeja
= super(TrivialEJA
, self
)
2193 # The rank is zero using my definition, namely the dimension of the
2194 # largest subalgebra generated by any element.
2195 fdeja
.__init
__(field
, mult_table
, **kwargs
)
2196 self
.rank
.set_cache(0)