]> gitweb.michael.orlitzky.com - sage.d.git/blob - eja_algebra.py
2c689d79d3c85d5d6f4d0daca39ca3f0e649540d
[sage.d.git] / eja_algebra.py
1 """
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.
6 """
7
8 from itertools import izip, repeat
9
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.prandom import choice
17 from sage.misc.table import table
18 from sage.modules.free_module import FreeModule, VectorSpace
19 from sage.rings.all import (ZZ, QQ, RR, RLF, CLF,
20 PolynomialRing,
21 QuadraticField)
22 from mjo.eja.eja_element import FiniteDimensionalEuclideanJordanAlgebraElement
23 from mjo.eja.eja_utils import _mat2vec
24
25 class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
26 # This is an ugly hack needed to prevent the category framework
27 # from implementing a coercion from our base ring (e.g. the
28 # rationals) into the algebra. First of all -- such a coercion is
29 # nonsense to begin with. But more importantly, it tries to do so
30 # in the category of rings, and since our algebras aren't
31 # associative they generally won't be rings.
32 _no_generic_basering_coercion = True
33
34 def __init__(self,
35 field,
36 mult_table,
37 rank,
38 prefix='e',
39 category=None,
40 natural_basis=None,
41 check=True):
42 """
43 SETUP::
44
45 sage: from mjo.eja.eja_algebra import (JordanSpinEJA, random_eja)
46
47 EXAMPLES:
48
49 By definition, Jordan multiplication commutes::
50
51 sage: set_random_seed()
52 sage: J = random_eja()
53 sage: x,y = J.random_elements(2)
54 sage: x*y == y*x
55 True
56
57 TESTS:
58
59 The ``field`` we're given must be real::
60
61 sage: JordanSpinEJA(2,QQbar)
62 Traceback (most recent call last):
63 ...
64 ValueError: field is not real
65
66 """
67 if check:
68 if not field.is_subring(RR):
69 # Note: this does return true for the real algebraic
70 # field, and any quadratic field where we've specified
71 # a real embedding.
72 raise ValueError('field is not real')
73
74 self._rank = rank
75 self._natural_basis = natural_basis
76
77 if category is None:
78 category = MagmaticAlgebras(field).FiniteDimensional()
79 category = category.WithBasis().Unital()
80
81 fda = super(FiniteDimensionalEuclideanJordanAlgebra, self)
82 fda.__init__(field,
83 range(len(mult_table)),
84 prefix=prefix,
85 category=category)
86 self.print_options(bracket='')
87
88 # The multiplication table we're given is necessarily in terms
89 # of vectors, because we don't have an algebra yet for
90 # anything to be an element of. However, it's faster in the
91 # long run to have the multiplication table be in terms of
92 # algebra elements. We do this after calling the superclass
93 # constructor so that from_vector() knows what to do.
94 self._multiplication_table = [ map(lambda x: self.from_vector(x), ls)
95 for ls in mult_table ]
96
97
98 def _element_constructor_(self, elt):
99 """
100 Construct an element of this algebra from its natural
101 representation.
102
103 This gets called only after the parent element _call_ method
104 fails to find a coercion for the argument.
105
106 SETUP::
107
108 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
109 ....: RealCartesianProductEJA,
110 ....: RealSymmetricEJA)
111
112 EXAMPLES:
113
114 The identity in `S^n` is converted to the identity in the EJA::
115
116 sage: J = RealSymmetricEJA(3)
117 sage: I = matrix.identity(QQ,3)
118 sage: J(I) == J.one()
119 True
120
121 This skew-symmetric matrix can't be represented in the EJA::
122
123 sage: J = RealSymmetricEJA(3)
124 sage: A = matrix(QQ,3, lambda i,j: i-j)
125 sage: J(A)
126 Traceback (most recent call last):
127 ...
128 ArithmeticError: vector is not in free module
129
130 TESTS:
131
132 Ensure that we can convert any element of the two non-matrix
133 simple algebras (whose natural representations are their usual
134 vector representations) back and forth faithfully::
135
136 sage: set_random_seed()
137 sage: J = RealCartesianProductEJA.random_instance()
138 sage: x = J.random_element()
139 sage: J(x.to_vector().column()) == x
140 True
141 sage: J = JordanSpinEJA.random_instance()
142 sage: x = J.random_element()
143 sage: J(x.to_vector().column()) == x
144 True
145
146 """
147 if elt == 0:
148 # The superclass implementation of random_element()
149 # needs to be able to coerce "0" into the algebra.
150 return self.zero()
151
152 natural_basis = self.natural_basis()
153 basis_space = natural_basis[0].matrix_space()
154 if elt not in basis_space:
155 raise ValueError("not a naturally-represented algebra element")
156
157 # Thanks for nothing! Matrix spaces aren't vector spaces in
158 # Sage, so we have to figure out its natural-basis coordinates
159 # ourselves. We use the basis space's ring instead of the
160 # element's ring because the basis space might be an algebraic
161 # closure whereas the base ring of the 3-by-3 identity matrix
162 # could be QQ instead of QQbar.
163 V = VectorSpace(basis_space.base_ring(), elt.nrows()*elt.ncols())
164 W = V.span_of_basis( _mat2vec(s) for s in natural_basis )
165 coords = W.coordinate_vector(_mat2vec(elt))
166 return self.from_vector(coords)
167
168
169 def _repr_(self):
170 """
171 Return a string representation of ``self``.
172
173 SETUP::
174
175 sage: from mjo.eja.eja_algebra import JordanSpinEJA
176
177 TESTS:
178
179 Ensure that it says what we think it says::
180
181 sage: JordanSpinEJA(2, field=QQ)
182 Euclidean Jordan algebra of dimension 2 over Rational Field
183 sage: JordanSpinEJA(3, field=RDF)
184 Euclidean Jordan algebra of dimension 3 over Real Double Field
185
186 """
187 fmt = "Euclidean Jordan algebra of dimension {} over {}"
188 return fmt.format(self.dimension(), self.base_ring())
189
190 def product_on_basis(self, i, j):
191 return self._multiplication_table[i][j]
192
193 def _a_regular_element(self):
194 """
195 Guess a regular element. Needed to compute the basis for our
196 characteristic polynomial coefficients.
197
198 SETUP::
199
200 sage: from mjo.eja.eja_algebra import random_eja
201
202 TESTS:
203
204 Ensure that this hacky method succeeds for every algebra that we
205 know how to construct::
206
207 sage: set_random_seed()
208 sage: J = random_eja()
209 sage: J._a_regular_element().is_regular()
210 True
211
212 """
213 gs = self.gens()
214 z = self.sum( (i+1)*gs[i] for i in range(len(gs)) )
215 if not z.is_regular():
216 raise ValueError("don't know a regular element")
217 return z
218
219
220 @cached_method
221 def _charpoly_basis_space(self):
222 """
223 Return the vector space spanned by the basis used in our
224 characteristic polynomial coefficients. This is used not only to
225 compute those coefficients, but also any time we need to
226 evaluate the coefficients (like when we compute the trace or
227 determinant).
228 """
229 z = self._a_regular_element()
230 # Don't use the parent vector space directly here in case this
231 # happens to be a subalgebra. In that case, we would be e.g.
232 # two-dimensional but span_of_basis() would expect three
233 # coordinates.
234 V = VectorSpace(self.base_ring(), self.vector_space().dimension())
235 basis = [ (z**k).to_vector() for k in range(self.rank()) ]
236 V1 = V.span_of_basis( basis )
237 b = (V1.basis() + V1.complement().basis())
238 return V.span_of_basis(b)
239
240
241
242 @cached_method
243 def _charpoly_coeff(self, i):
244 """
245 Return the coefficient polynomial "a_{i}" of this algebra's
246 general characteristic polynomial.
247
248 Having this be a separate cached method lets us compute and
249 store the trace/determinant (a_{r-1} and a_{0} respectively)
250 separate from the entire characteristic polynomial.
251 """
252 (A_of_x, x, xr, detA) = self._charpoly_matrix_system()
253 R = A_of_x.base_ring()
254 if i >= self.rank():
255 # Guaranteed by theory
256 return R.zero()
257
258 # Danger: the in-place modification is done for performance
259 # reasons (reconstructing a matrix with huge polynomial
260 # entries is slow), but I don't know how cached_method works,
261 # so it's highly possible that we're modifying some global
262 # list variable by reference, here. In other words, you
263 # probably shouldn't call this method twice on the same
264 # algebra, at the same time, in two threads
265 Ai_orig = A_of_x.column(i)
266 A_of_x.set_column(i,xr)
267 numerator = A_of_x.det()
268 A_of_x.set_column(i,Ai_orig)
269
270 # We're relying on the theory here to ensure that each a_i is
271 # indeed back in R, and the added negative signs are to make
272 # the whole charpoly expression sum to zero.
273 return R(-numerator/detA)
274
275
276 @cached_method
277 def _charpoly_matrix_system(self):
278 """
279 Compute the matrix whose entries A_ij are polynomials in
280 X1,...,XN, the vector ``x`` of variables X1,...,XN, the vector
281 corresponding to `x^r` and the determinent of the matrix A =
282 [A_ij]. In other words, all of the fixed (cachable) data needed
283 to compute the coefficients of the characteristic polynomial.
284 """
285 r = self.rank()
286 n = self.dimension()
287
288 # Turn my vector space into a module so that "vectors" can
289 # have multivatiate polynomial entries.
290 names = tuple('X' + str(i) for i in range(1,n+1))
291 R = PolynomialRing(self.base_ring(), names)
292
293 # Using change_ring() on the parent's vector space doesn't work
294 # here because, in a subalgebra, that vector space has a basis
295 # and change_ring() tries to bring the basis along with it. And
296 # that doesn't work unless the new ring is a PID, which it usually
297 # won't be.
298 V = FreeModule(R,n)
299
300 # Now let x = (X1,X2,...,Xn) be the vector whose entries are
301 # indeterminates...
302 x = V(names)
303
304 # And figure out the "left multiplication by x" matrix in
305 # that setting.
306 lmbx_cols = []
307 monomial_matrices = [ self.monomial(i).operator().matrix()
308 for i in range(n) ] # don't recompute these!
309 for k in range(n):
310 ek = self.monomial(k).to_vector()
311 lmbx_cols.append(
312 sum( x[i]*(monomial_matrices[i]*ek)
313 for i in range(n) ) )
314 Lx = matrix.column(R, lmbx_cols)
315
316 # Now we can compute powers of x "symbolically"
317 x_powers = [self.one().to_vector(), x]
318 for d in range(2, r+1):
319 x_powers.append( Lx*(x_powers[-1]) )
320
321 idmat = matrix.identity(R, n)
322
323 W = self._charpoly_basis_space()
324 W = W.change_ring(R.fraction_field())
325
326 # Starting with the standard coordinates x = (X1,X2,...,Xn)
327 # and then converting the entries to W-coordinates allows us
328 # to pass in the standard coordinates to the charpoly and get
329 # back the right answer. Specifically, with x = (X1,X2,...,Xn),
330 # we have
331 #
332 # W.coordinates(x^2) eval'd at (standard z-coords)
333 # =
334 # W-coords of (z^2)
335 # =
336 # W-coords of (standard coords of x^2 eval'd at std-coords of z)
337 #
338 # We want the middle equivalent thing in our matrix, but use
339 # the first equivalent thing instead so that we can pass in
340 # standard coordinates.
341 x_powers = [ W.coordinate_vector(xp) for xp in x_powers ]
342 l2 = [idmat.column(k-1) for k in range(r+1, n+1)]
343 A_of_x = matrix.column(R, n, (x_powers[:r] + l2))
344 return (A_of_x, x, x_powers[r], A_of_x.det())
345
346
347 @cached_method
348 def characteristic_polynomial(self):
349 """
350 Return a characteristic polynomial that works for all elements
351 of this algebra.
352
353 The resulting polynomial has `n+1` variables, where `n` is the
354 dimension of this algebra. The first `n` variables correspond to
355 the coordinates of an algebra element: when evaluated at the
356 coordinates of an algebra element with respect to a certain
357 basis, the result is a univariate polynomial (in the one
358 remaining variable ``t``), namely the characteristic polynomial
359 of that element.
360
361 SETUP::
362
363 sage: from mjo.eja.eja_algebra import JordanSpinEJA
364
365 EXAMPLES:
366
367 The characteristic polynomial in the spin algebra is given in
368 Alizadeh, Example 11.11::
369
370 sage: J = JordanSpinEJA(3)
371 sage: p = J.characteristic_polynomial(); p
372 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
373 sage: xvec = J.one().to_vector()
374 sage: p(*xvec)
375 t^2 - 2*t + 1
376
377 """
378 r = self.rank()
379 n = self.dimension()
380
381 # The list of coefficient polynomials a_1, a_2, ..., a_n.
382 a = [ self._charpoly_coeff(i) for i in range(n) ]
383
384 # We go to a bit of trouble here to reorder the
385 # indeterminates, so that it's easier to evaluate the
386 # characteristic polynomial at x's coordinates and get back
387 # something in terms of t, which is what we want.
388 R = a[0].parent()
389 S = PolynomialRing(self.base_ring(),'t')
390 t = S.gen(0)
391 S = PolynomialRing(S, R.variable_names())
392 t = S(t)
393
394 # Note: all entries past the rth should be zero. The
395 # coefficient of the highest power (x^r) is 1, but it doesn't
396 # appear in the solution vector which contains coefficients
397 # for the other powers (to make them sum to x^r).
398 if (r < n):
399 a[r] = 1 # corresponds to x^r
400 else:
401 # When the rank is equal to the dimension, trying to
402 # assign a[r] goes out-of-bounds.
403 a.append(1) # corresponds to x^r
404
405 return sum( a[k]*(t**k) for k in xrange(len(a)) )
406
407
408 def inner_product(self, x, y):
409 """
410 The inner product associated with this Euclidean Jordan algebra.
411
412 Defaults to the trace inner product, but can be overridden by
413 subclasses if they are sure that the necessary properties are
414 satisfied.
415
416 SETUP::
417
418 sage: from mjo.eja.eja_algebra import random_eja
419
420 EXAMPLES:
421
422 Our inner product is "associative," which means the following for
423 a symmetric bilinear form::
424
425 sage: set_random_seed()
426 sage: J = random_eja()
427 sage: x,y,z = J.random_elements(3)
428 sage: (x*y).inner_product(z) == y.inner_product(x*z)
429 True
430
431 """
432 X = x.natural_representation()
433 Y = y.natural_representation()
434 return self.natural_inner_product(X,Y)
435
436
437 def is_trivial(self):
438 """
439 Return whether or not this algebra is trivial.
440
441 A trivial algebra contains only the zero element.
442
443 SETUP::
444
445 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
446
447 EXAMPLES::
448
449 sage: J = ComplexHermitianEJA(3)
450 sage: J.is_trivial()
451 False
452 sage: A = J.zero().subalgebra_generated_by()
453 sage: A.is_trivial()
454 True
455
456 """
457 return self.dimension() == 0
458
459
460 def multiplication_table(self):
461 """
462 Return a visual representation of this algebra's multiplication
463 table (on basis elements).
464
465 SETUP::
466
467 sage: from mjo.eja.eja_algebra import JordanSpinEJA
468
469 EXAMPLES::
470
471 sage: J = JordanSpinEJA(4)
472 sage: J.multiplication_table()
473 +----++----+----+----+----+
474 | * || e0 | e1 | e2 | e3 |
475 +====++====+====+====+====+
476 | e0 || e0 | e1 | e2 | e3 |
477 +----++----+----+----+----+
478 | e1 || e1 | e0 | 0 | 0 |
479 +----++----+----+----+----+
480 | e2 || e2 | 0 | e0 | 0 |
481 +----++----+----+----+----+
482 | e3 || e3 | 0 | 0 | e0 |
483 +----++----+----+----+----+
484
485 """
486 M = list(self._multiplication_table) # copy
487 for i in xrange(len(M)):
488 # M had better be "square"
489 M[i] = [self.monomial(i)] + M[i]
490 M = [["*"] + list(self.gens())] + M
491 return table(M, header_row=True, header_column=True, frame=True)
492
493
494 def natural_basis(self):
495 """
496 Return a more-natural representation of this algebra's basis.
497
498 Every finite-dimensional Euclidean Jordan Algebra is a direct
499 sum of five simple algebras, four of which comprise Hermitian
500 matrices. This method returns the original "natural" basis
501 for our underlying vector space. (Typically, the natural basis
502 is used to construct the multiplication table in the first place.)
503
504 Note that this will always return a matrix. The standard basis
505 in `R^n` will be returned as `n`-by-`1` column matrices.
506
507 SETUP::
508
509 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
510 ....: RealSymmetricEJA)
511
512 EXAMPLES::
513
514 sage: J = RealSymmetricEJA(2)
515 sage: J.basis()
516 Finite family {0: e0, 1: e1, 2: e2}
517 sage: J.natural_basis()
518 (
519 [1 0] [ 0 1/2*sqrt2] [0 0]
520 [0 0], [1/2*sqrt2 0], [0 1]
521 )
522
523 ::
524
525 sage: J = JordanSpinEJA(2)
526 sage: J.basis()
527 Finite family {0: e0, 1: e1}
528 sage: J.natural_basis()
529 (
530 [1] [0]
531 [0], [1]
532 )
533
534 """
535 if self._natural_basis is None:
536 M = self.natural_basis_space()
537 return tuple( M(b.to_vector()) for b in self.basis() )
538 else:
539 return self._natural_basis
540
541
542 def natural_basis_space(self):
543 """
544 Return the matrix space in which this algebra's natural basis
545 elements live.
546 """
547 if self._natural_basis is None or len(self._natural_basis) == 0:
548 return MatrixSpace(self.base_ring(), self.dimension(), 1)
549 else:
550 return self._natural_basis[0].matrix_space()
551
552
553 @staticmethod
554 def natural_inner_product(X,Y):
555 """
556 Compute the inner product of two naturally-represented elements.
557
558 For example in the real symmetric matrix EJA, this will compute
559 the trace inner-product of two n-by-n symmetric matrices. The
560 default should work for the real cartesian product EJA, the
561 Jordan spin EJA, and the real symmetric matrices. The others
562 will have to be overridden.
563 """
564 return (X.conjugate_transpose()*Y).trace()
565
566
567 @cached_method
568 def one(self):
569 """
570 Return the unit element of this algebra.
571
572 SETUP::
573
574 sage: from mjo.eja.eja_algebra import (RealCartesianProductEJA,
575 ....: random_eja)
576
577 EXAMPLES::
578
579 sage: J = RealCartesianProductEJA(5)
580 sage: J.one()
581 e0 + e1 + e2 + e3 + e4
582
583 TESTS:
584
585 The identity element acts like the identity::
586
587 sage: set_random_seed()
588 sage: J = random_eja()
589 sage: x = J.random_element()
590 sage: J.one()*x == x and x*J.one() == x
591 True
592
593 The matrix of the unit element's operator is the identity::
594
595 sage: set_random_seed()
596 sage: J = random_eja()
597 sage: actual = J.one().operator().matrix()
598 sage: expected = matrix.identity(J.base_ring(), J.dimension())
599 sage: actual == expected
600 True
601
602 """
603 # We can brute-force compute the matrices of the operators
604 # that correspond to the basis elements of this algebra.
605 # If some linear combination of those basis elements is the
606 # algebra identity, then the same linear combination of
607 # their matrices has to be the identity matrix.
608 #
609 # Of course, matrices aren't vectors in sage, so we have to
610 # appeal to the "long vectors" isometry.
611 oper_vecs = [ _mat2vec(g.operator().matrix()) for g in self.gens() ]
612
613 # Now we use basis linear algebra to find the coefficients,
614 # of the matrices-as-vectors-linear-combination, which should
615 # work for the original algebra basis too.
616 A = matrix.column(self.base_ring(), oper_vecs)
617
618 # We used the isometry on the left-hand side already, but we
619 # still need to do it for the right-hand side. Recall that we
620 # wanted something that summed to the identity matrix.
621 b = _mat2vec( matrix.identity(self.base_ring(), self.dimension()) )
622
623 # Now if there's an identity element in the algebra, this should work.
624 coeffs = A.solve_right(b)
625 return self.linear_combination(zip(self.gens(), coeffs))
626
627
628 def random_element(self):
629 # Temporary workaround for https://trac.sagemath.org/ticket/28327
630 if self.is_trivial():
631 return self.zero()
632 else:
633 s = super(FiniteDimensionalEuclideanJordanAlgebra, self)
634 return s.random_element()
635
636 def random_elements(self, count):
637 """
638 Return ``count`` random elements as a tuple.
639
640 SETUP::
641
642 sage: from mjo.eja.eja_algebra import JordanSpinEJA
643
644 EXAMPLES::
645
646 sage: J = JordanSpinEJA(3)
647 sage: x,y,z = J.random_elements(3)
648 sage: all( [ x in J, y in J, z in J ])
649 True
650 sage: len( J.random_elements(10) ) == 10
651 True
652
653 """
654 return tuple( self.random_element() for idx in xrange(count) )
655
656
657 def rank(self):
658 """
659 Return the rank of this EJA.
660
661 ALGORITHM:
662
663 The author knows of no algorithm to compute the rank of an EJA
664 where only the multiplication table is known. In lieu of one, we
665 require the rank to be specified when the algebra is created,
666 and simply pass along that number here.
667
668 SETUP::
669
670 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
671 ....: RealSymmetricEJA,
672 ....: ComplexHermitianEJA,
673 ....: QuaternionHermitianEJA,
674 ....: random_eja)
675
676 EXAMPLES:
677
678 The rank of the Jordan spin algebra is always two::
679
680 sage: JordanSpinEJA(2).rank()
681 2
682 sage: JordanSpinEJA(3).rank()
683 2
684 sage: JordanSpinEJA(4).rank()
685 2
686
687 The rank of the `n`-by-`n` Hermitian real, complex, or
688 quaternion matrices is `n`::
689
690 sage: RealSymmetricEJA(4).rank()
691 4
692 sage: ComplexHermitianEJA(3).rank()
693 3
694 sage: QuaternionHermitianEJA(2).rank()
695 2
696
697 TESTS:
698
699 Ensure that every EJA that we know how to construct has a
700 positive integer rank::
701
702 sage: set_random_seed()
703 sage: r = random_eja().rank()
704 sage: r in ZZ and r > 0
705 True
706
707 """
708 return self._rank
709
710
711 def vector_space(self):
712 """
713 Return the vector space that underlies this algebra.
714
715 SETUP::
716
717 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
718
719 EXAMPLES::
720
721 sage: J = RealSymmetricEJA(2)
722 sage: J.vector_space()
723 Vector space of dimension 3 over...
724
725 """
726 return self.zero().to_vector().parent().ambient_vector_space()
727
728
729 Element = FiniteDimensionalEuclideanJordanAlgebraElement
730
731
732 class KnownRankEJA(object):
733 """
734 A class for algebras that we actually know we can construct. The
735 main issue is that, for most of our methods to make sense, we need
736 to know the rank of our algebra. Thus we can't simply generate a
737 "random" algebra, or even check that a given basis and product
738 satisfy the axioms; because even if everything looks OK, we wouldn't
739 know the rank we need to actuallty build the thing.
740
741 Not really a subclass of FDEJA because doing that causes method
742 resolution errors, e.g.
743
744 TypeError: Error when calling the metaclass bases
745 Cannot create a consistent method resolution
746 order (MRO) for bases FiniteDimensionalEuclideanJordanAlgebra,
747 KnownRankEJA
748
749 """
750 @staticmethod
751 def _max_test_case_size():
752 """
753 Return an integer "size" that is an upper bound on the size of
754 this algebra when it is used in a random test
755 case. Unfortunately, the term "size" is quite vague -- when
756 dealing with `R^n` under either the Hadamard or Jordan spin
757 product, the "size" refers to the dimension `n`. When dealing
758 with a matrix algebra (real symmetric or complex/quaternion
759 Hermitian), it refers to the size of the matrix, which is
760 far less than the dimension of the underlying vector space.
761
762 We default to five in this class, which is safe in `R^n`. The
763 matrix algebra subclasses (or any class where the "size" is
764 interpreted to be far less than the dimension) should override
765 with a smaller number.
766 """
767 return 5
768
769 @classmethod
770 def random_instance(cls, field=QQ, **kwargs):
771 """
772 Return a random instance of this type of algebra.
773
774 Beware, this will crash for "most instances" because the
775 constructor below looks wrong.
776 """
777 n = ZZ.random_element(cls._max_test_case_size()) + 1
778 return cls(n, field, **kwargs)
779
780
781 class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra,
782 KnownRankEJA):
783 """
784 Return the Euclidean Jordan Algebra corresponding to the set
785 `R^n` under the Hadamard product.
786
787 Note: this is nothing more than the Cartesian product of ``n``
788 copies of the spin algebra. Once Cartesian product algebras
789 are implemented, this can go.
790
791 SETUP::
792
793 sage: from mjo.eja.eja_algebra import RealCartesianProductEJA
794
795 EXAMPLES:
796
797 This multiplication table can be verified by hand::
798
799 sage: J = RealCartesianProductEJA(3)
800 sage: e0,e1,e2 = J.gens()
801 sage: e0*e0
802 e0
803 sage: e0*e1
804 0
805 sage: e0*e2
806 0
807 sage: e1*e1
808 e1
809 sage: e1*e2
810 0
811 sage: e2*e2
812 e2
813
814 TESTS:
815
816 We can change the generator prefix::
817
818 sage: RealCartesianProductEJA(3, prefix='r').gens()
819 (r0, r1, r2)
820
821 """
822 def __init__(self, n, field=QQ, **kwargs):
823 V = VectorSpace(field, n)
824 mult_table = [ [ V.gen(i)*(i == j) for j in xrange(n) ]
825 for i in xrange(n) ]
826
827 fdeja = super(RealCartesianProductEJA, self)
828 return fdeja.__init__(field, mult_table, rank=n, **kwargs)
829
830 def inner_product(self, x, y):
831 """
832 Faster to reimplement than to use natural representations.
833
834 SETUP::
835
836 sage: from mjo.eja.eja_algebra import RealCartesianProductEJA
837
838 TESTS:
839
840 Ensure that this is the usual inner product for the algebras
841 over `R^n`::
842
843 sage: set_random_seed()
844 sage: J = RealCartesianProductEJA.random_instance()
845 sage: x,y = J.random_elements(2)
846 sage: X = x.natural_representation()
847 sage: Y = y.natural_representation()
848 sage: x.inner_product(y) == J.natural_inner_product(X,Y)
849 True
850
851 """
852 return x.to_vector().inner_product(y.to_vector())
853
854
855 def random_eja(field=QQ):
856 """
857 Return a "random" finite-dimensional Euclidean Jordan Algebra.
858
859 ALGORITHM:
860
861 For now, we choose a random natural number ``n`` (greater than zero)
862 and then give you back one of the following:
863
864 * The cartesian product of the rational numbers ``n`` times; this is
865 ``QQ^n`` with the Hadamard product.
866
867 * The Jordan spin algebra on ``QQ^n``.
868
869 * The ``n``-by-``n`` rational symmetric matrices with the symmetric
870 product.
871
872 * The ``n``-by-``n`` complex-rational Hermitian matrices embedded
873 in the space of ``2n``-by-``2n`` real symmetric matrices.
874
875 * The ``n``-by-``n`` quaternion-rational Hermitian matrices embedded
876 in the space of ``4n``-by-``4n`` real symmetric matrices.
877
878 Later this might be extended to return Cartesian products of the
879 EJAs above.
880
881 SETUP::
882
883 sage: from mjo.eja.eja_algebra import random_eja
884
885 TESTS::
886
887 sage: random_eja()
888 Euclidean Jordan algebra of dimension...
889
890 """
891 classname = choice(KnownRankEJA.__subclasses__())
892 return classname.random_instance(field=field)
893
894
895
896
897
898
899 class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
900 @staticmethod
901 def _max_test_case_size():
902 # Play it safe, since this will be squared and the underlying
903 # field can have dimension 4 (quaternions) too.
904 return 2
905
906 def __init__(self, field, basis, rank, normalize_basis=True, **kwargs):
907 """
908 Compared to the superclass constructor, we take a basis instead of
909 a multiplication table because the latter can be computed in terms
910 of the former when the product is known (like it is here).
911 """
912 # Used in this class's fast _charpoly_coeff() override.
913 self._basis_normalizers = None
914
915 # We're going to loop through this a few times, so now's a good
916 # time to ensure that it isn't a generator expression.
917 basis = tuple(basis)
918
919 if rank > 1 and normalize_basis:
920 # We'll need sqrt(2) to normalize the basis, and this
921 # winds up in the multiplication table, so the whole
922 # algebra needs to be over the field extension.
923 R = PolynomialRing(field, 'z')
924 z = R.gen()
925 p = z**2 - 2
926 if p.is_irreducible():
927 field = field.extension(p, 'sqrt2', embedding=RLF(2).sqrt())
928 basis = tuple( s.change_ring(field) for s in basis )
929 self._basis_normalizers = tuple(
930 ~(self.natural_inner_product(s,s).sqrt()) for s in basis )
931 basis = tuple(s*c for (s,c) in izip(basis,self._basis_normalizers))
932
933 Qs = self.multiplication_table_from_matrix_basis(basis)
934
935 fdeja = super(MatrixEuclideanJordanAlgebra, self)
936 return fdeja.__init__(field,
937 Qs,
938 rank=rank,
939 natural_basis=basis,
940 **kwargs)
941
942
943 @cached_method
944 def _charpoly_coeff(self, i):
945 """
946 Override the parent method with something that tries to compute
947 over a faster (non-extension) field.
948 """
949 if self._basis_normalizers is None:
950 # We didn't normalize, so assume that the basis we started
951 # with had entries in a nice field.
952 return super(MatrixEuclideanJordanAlgebra, self)._charpoly_coeff(i)
953 else:
954 basis = ( (b/n) for (b,n) in izip(self.natural_basis(),
955 self._basis_normalizers) )
956
957 # Do this over the rationals and convert back at the end.
958 J = MatrixEuclideanJordanAlgebra(QQ,
959 basis,
960 self.rank(),
961 normalize_basis=False)
962 (_,x,_,_) = J._charpoly_matrix_system()
963 p = J._charpoly_coeff(i)
964 # p might be missing some vars, have to substitute "optionally"
965 pairs = izip(x.base_ring().gens(), self._basis_normalizers)
966 substitutions = { v: v*c for (v,c) in pairs }
967 result = p.subs(substitutions)
968
969 # The result of "subs" can be either a coefficient-ring
970 # element or a polynomial. Gotta handle both cases.
971 if result in QQ:
972 return self.base_ring()(result)
973 else:
974 return result.change_ring(self.base_ring())
975
976
977 @staticmethod
978 def multiplication_table_from_matrix_basis(basis):
979 """
980 At least three of the five simple Euclidean Jordan algebras have the
981 symmetric multiplication (A,B) |-> (AB + BA)/2, where the
982 multiplication on the right is matrix multiplication. Given a basis
983 for the underlying matrix space, this function returns a
984 multiplication table (obtained by looping through the basis
985 elements) for an algebra of those matrices.
986 """
987 # In S^2, for example, we nominally have four coordinates even
988 # though the space is of dimension three only. The vector space V
989 # is supposed to hold the entire long vector, and the subspace W
990 # of V will be spanned by the vectors that arise from symmetric
991 # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
992 field = basis[0].base_ring()
993 dimension = basis[0].nrows()
994
995 V = VectorSpace(field, dimension**2)
996 W = V.span_of_basis( _mat2vec(s) for s in basis )
997 n = len(basis)
998 mult_table = [[W.zero() for j in xrange(n)] for i in xrange(n)]
999 for i in xrange(n):
1000 for j in xrange(n):
1001 mat_entry = (basis[i]*basis[j] + basis[j]*basis[i])/2
1002 mult_table[i][j] = W.coordinate_vector(_mat2vec(mat_entry))
1003
1004 return mult_table
1005
1006
1007 @staticmethod
1008 def real_embed(M):
1009 """
1010 Embed the matrix ``M`` into a space of real matrices.
1011
1012 The matrix ``M`` can have entries in any field at the moment:
1013 the real numbers, complex numbers, or quaternions. And although
1014 they are not a field, we can probably support octonions at some
1015 point, too. This function returns a real matrix that "acts like"
1016 the original with respect to matrix multiplication; i.e.
1017
1018 real_embed(M*N) = real_embed(M)*real_embed(N)
1019
1020 """
1021 raise NotImplementedError
1022
1023
1024 @staticmethod
1025 def real_unembed(M):
1026 """
1027 The inverse of :meth:`real_embed`.
1028 """
1029 raise NotImplementedError
1030
1031
1032 @classmethod
1033 def natural_inner_product(cls,X,Y):
1034 Xu = cls.real_unembed(X)
1035 Yu = cls.real_unembed(Y)
1036 tr = (Xu*Yu).trace()
1037
1038 if tr in RLF:
1039 # It's real already.
1040 return tr
1041
1042 # Otherwise, try the thing that works for complex numbers; and
1043 # if that doesn't work, the thing that works for quaternions.
1044 try:
1045 return tr.vector()[0] # real part, imag part is index 1
1046 except AttributeError:
1047 # A quaternions doesn't have a vector() method, but does
1048 # have coefficient_tuple() method that returns the
1049 # coefficients of 1, i, j, and k -- in that order.
1050 return tr.coefficient_tuple()[0]
1051
1052
1053 class RealMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
1054 @staticmethod
1055 def real_embed(M):
1056 """
1057 The identity function, for embedding real matrices into real
1058 matrices.
1059 """
1060 return M
1061
1062 @staticmethod
1063 def real_unembed(M):
1064 """
1065 The identity function, for unembedding real matrices from real
1066 matrices.
1067 """
1068 return M
1069
1070
1071 class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra, KnownRankEJA):
1072 """
1073 The rank-n simple EJA consisting of real symmetric n-by-n
1074 matrices, the usual symmetric Jordan product, and the trace inner
1075 product. It has dimension `(n^2 + n)/2` over the reals.
1076
1077 SETUP::
1078
1079 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1080
1081 EXAMPLES::
1082
1083 sage: J = RealSymmetricEJA(2)
1084 sage: e0, e1, e2 = J.gens()
1085 sage: e0*e0
1086 e0
1087 sage: e1*e1
1088 1/2*e0 + 1/2*e2
1089 sage: e2*e2
1090 e2
1091
1092 In theory, our "field" can be any subfield of the reals::
1093
1094 sage: RealSymmetricEJA(2, AA)
1095 Euclidean Jordan algebra of dimension 3 over Algebraic Real Field
1096 sage: RealSymmetricEJA(2, RR)
1097 Euclidean Jordan algebra of dimension 3 over Real Field with
1098 53 bits of precision
1099
1100 TESTS:
1101
1102 The dimension of this algebra is `(n^2 + n) / 2`::
1103
1104 sage: set_random_seed()
1105 sage: n_max = RealSymmetricEJA._max_test_case_size()
1106 sage: n = ZZ.random_element(1, n_max)
1107 sage: J = RealSymmetricEJA(n)
1108 sage: J.dimension() == (n^2 + n)/2
1109 True
1110
1111 The Jordan multiplication is what we think it is::
1112
1113 sage: set_random_seed()
1114 sage: J = RealSymmetricEJA.random_instance()
1115 sage: x,y = J.random_elements(2)
1116 sage: actual = (x*y).natural_representation()
1117 sage: X = x.natural_representation()
1118 sage: Y = y.natural_representation()
1119 sage: expected = (X*Y + Y*X)/2
1120 sage: actual == expected
1121 True
1122 sage: J(expected) == x*y
1123 True
1124
1125 We can change the generator prefix::
1126
1127 sage: RealSymmetricEJA(3, prefix='q').gens()
1128 (q0, q1, q2, q3, q4, q5)
1129
1130 Our natural basis is normalized with respect to the natural inner
1131 product unless we specify otherwise::
1132
1133 sage: set_random_seed()
1134 sage: J = RealSymmetricEJA.random_instance()
1135 sage: all( b.norm() == 1 for b in J.gens() )
1136 True
1137
1138 Since our natural basis is normalized with respect to the natural
1139 inner product, and since we know that this algebra is an EJA, any
1140 left-multiplication operator's matrix will be symmetric because
1141 natural->EJA basis representation is an isometry and within the EJA
1142 the operator is self-adjoint by the Jordan axiom::
1143
1144 sage: set_random_seed()
1145 sage: x = RealSymmetricEJA.random_instance().random_element()
1146 sage: x.operator().matrix().is_symmetric()
1147 True
1148
1149 """
1150 @classmethod
1151 def _denormalized_basis(cls, n, field):
1152 """
1153 Return a basis for the space of real symmetric n-by-n matrices.
1154
1155 SETUP::
1156
1157 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1158
1159 TESTS::
1160
1161 sage: set_random_seed()
1162 sage: n = ZZ.random_element(1,5)
1163 sage: B = RealSymmetricEJA._denormalized_basis(n,QQ)
1164 sage: all( M.is_symmetric() for M in B)
1165 True
1166
1167 """
1168 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1169 # coordinates.
1170 S = []
1171 for i in xrange(n):
1172 for j in xrange(i+1):
1173 Eij = matrix(field, n, lambda k,l: k==i and l==j)
1174 if i == j:
1175 Sij = Eij
1176 else:
1177 Sij = Eij + Eij.transpose()
1178 S.append(Sij)
1179 return S
1180
1181
1182 @staticmethod
1183 def _max_test_case_size():
1184 return 4 # Dimension 10
1185
1186
1187 def __init__(self, n, field=QQ, **kwargs):
1188 basis = self._denormalized_basis(n, field)
1189 super(RealSymmetricEJA, self).__init__(field, basis, n, **kwargs)
1190
1191
1192 class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
1193 @staticmethod
1194 def real_embed(M):
1195 """
1196 Embed the n-by-n complex matrix ``M`` into the space of real
1197 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1198 bi` to the block matrix ``[[a,b],[-b,a]]``.
1199
1200 SETUP::
1201
1202 sage: from mjo.eja.eja_algebra import \
1203 ....: ComplexMatrixEuclideanJordanAlgebra
1204
1205 EXAMPLES::
1206
1207 sage: F = QuadraticField(-1, 'i')
1208 sage: x1 = F(4 - 2*i)
1209 sage: x2 = F(1 + 2*i)
1210 sage: x3 = F(-i)
1211 sage: x4 = F(6)
1212 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1213 sage: ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1214 [ 4 -2| 1 2]
1215 [ 2 4|-2 1]
1216 [-----+-----]
1217 [ 0 -1| 6 0]
1218 [ 1 0| 0 6]
1219
1220 TESTS:
1221
1222 Embedding is a homomorphism (isomorphism, in fact)::
1223
1224 sage: set_random_seed()
1225 sage: n_max = ComplexMatrixEuclideanJordanAlgebra._max_test_case_size()
1226 sage: n = ZZ.random_element(n_max)
1227 sage: F = QuadraticField(-1, 'i')
1228 sage: X = random_matrix(F, n)
1229 sage: Y = random_matrix(F, n)
1230 sage: Xe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X)
1231 sage: Ye = ComplexMatrixEuclideanJordanAlgebra.real_embed(Y)
1232 sage: XYe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1233 sage: Xe*Ye == XYe
1234 True
1235
1236 """
1237 n = M.nrows()
1238 if M.ncols() != n:
1239 raise ValueError("the matrix 'M' must be square")
1240
1241 # We don't need any adjoined elements...
1242 field = M.base_ring().base_ring()
1243
1244 blocks = []
1245 for z in M.list():
1246 a = z.list()[0] # real part, I guess
1247 b = z.list()[1] # imag part, I guess
1248 blocks.append(matrix(field, 2, [[a,b],[-b,a]]))
1249
1250 return matrix.block(field, n, blocks)
1251
1252
1253 @staticmethod
1254 def real_unembed(M):
1255 """
1256 The inverse of _embed_complex_matrix().
1257
1258 SETUP::
1259
1260 sage: from mjo.eja.eja_algebra import \
1261 ....: ComplexMatrixEuclideanJordanAlgebra
1262
1263 EXAMPLES::
1264
1265 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1266 ....: [-2, 1, -4, 3],
1267 ....: [ 9, 10, 11, 12],
1268 ....: [-10, 9, -12, 11] ])
1269 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(A)
1270 [ 2*i + 1 4*i + 3]
1271 [ 10*i + 9 12*i + 11]
1272
1273 TESTS:
1274
1275 Unembedding is the inverse of embedding::
1276
1277 sage: set_random_seed()
1278 sage: F = QuadraticField(-1, 'i')
1279 sage: M = random_matrix(F, 3)
1280 sage: Me = ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1281 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1282 True
1283
1284 """
1285 n = ZZ(M.nrows())
1286 if M.ncols() != n:
1287 raise ValueError("the matrix 'M' must be square")
1288 if not n.mod(2).is_zero():
1289 raise ValueError("the matrix 'M' must be a complex embedding")
1290
1291 # If "M" was normalized, its base ring might have roots
1292 # adjoined and they can stick around after unembedding.
1293 field = M.base_ring()
1294 R = PolynomialRing(field, 'z')
1295 z = R.gen()
1296 F = field.extension(z**2 + 1, 'i', embedding=CLF(-1).sqrt())
1297 i = F.gen()
1298
1299 # Go top-left to bottom-right (reading order), converting every
1300 # 2-by-2 block we see to a single complex element.
1301 elements = []
1302 for k in xrange(n/2):
1303 for j in xrange(n/2):
1304 submat = M[2*k:2*k+2,2*j:2*j+2]
1305 if submat[0,0] != submat[1,1]:
1306 raise ValueError('bad on-diagonal submatrix')
1307 if submat[0,1] != -submat[1,0]:
1308 raise ValueError('bad off-diagonal submatrix')
1309 z = submat[0,0] + submat[0,1]*i
1310 elements.append(z)
1311
1312 return matrix(F, n/2, elements)
1313
1314
1315 @classmethod
1316 def natural_inner_product(cls,X,Y):
1317 """
1318 Compute a natural inner product in this algebra directly from
1319 its real embedding.
1320
1321 SETUP::
1322
1323 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1324
1325 TESTS:
1326
1327 This gives the same answer as the slow, default method implemented
1328 in :class:`MatrixEuclideanJordanAlgebra`::
1329
1330 sage: set_random_seed()
1331 sage: J = ComplexHermitianEJA.random_instance()
1332 sage: x,y = J.random_elements(2)
1333 sage: Xe = x.natural_representation()
1334 sage: Ye = y.natural_representation()
1335 sage: X = ComplexHermitianEJA.real_unembed(Xe)
1336 sage: Y = ComplexHermitianEJA.real_unembed(Ye)
1337 sage: expected = (X*Y).trace().vector()[0]
1338 sage: actual = ComplexHermitianEJA.natural_inner_product(Xe,Ye)
1339 sage: actual == expected
1340 True
1341
1342 """
1343 return RealMatrixEuclideanJordanAlgebra.natural_inner_product(X,Y)/2
1344
1345
1346 class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra, KnownRankEJA):
1347 """
1348 The rank-n simple EJA consisting of complex Hermitian n-by-n
1349 matrices over the real numbers, the usual symmetric Jordan product,
1350 and the real-part-of-trace inner product. It has dimension `n^2` over
1351 the reals.
1352
1353 SETUP::
1354
1355 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1356
1357 EXAMPLES:
1358
1359 In theory, our "field" can be any subfield of the reals::
1360
1361 sage: ComplexHermitianEJA(2, AA)
1362 Euclidean Jordan algebra of dimension 4 over Algebraic Real Field
1363 sage: ComplexHermitianEJA(2, RR)
1364 Euclidean Jordan algebra of dimension 4 over Real Field with
1365 53 bits of precision
1366
1367 TESTS:
1368
1369 The dimension of this algebra is `n^2`::
1370
1371 sage: set_random_seed()
1372 sage: n_max = ComplexHermitianEJA._max_test_case_size()
1373 sage: n = ZZ.random_element(1, n_max)
1374 sage: J = ComplexHermitianEJA(n)
1375 sage: J.dimension() == n^2
1376 True
1377
1378 The Jordan multiplication is what we think it is::
1379
1380 sage: set_random_seed()
1381 sage: J = ComplexHermitianEJA.random_instance()
1382 sage: x,y = J.random_elements(2)
1383 sage: actual = (x*y).natural_representation()
1384 sage: X = x.natural_representation()
1385 sage: Y = y.natural_representation()
1386 sage: expected = (X*Y + Y*X)/2
1387 sage: actual == expected
1388 True
1389 sage: J(expected) == x*y
1390 True
1391
1392 We can change the generator prefix::
1393
1394 sage: ComplexHermitianEJA(2, prefix='z').gens()
1395 (z0, z1, z2, z3)
1396
1397 Our natural basis is normalized with respect to the natural inner
1398 product unless we specify otherwise::
1399
1400 sage: set_random_seed()
1401 sage: J = ComplexHermitianEJA.random_instance()
1402 sage: all( b.norm() == 1 for b in J.gens() )
1403 True
1404
1405 Since our natural basis is normalized with respect to the natural
1406 inner product, and since we know that this algebra is an EJA, any
1407 left-multiplication operator's matrix will be symmetric because
1408 natural->EJA basis representation is an isometry and within the EJA
1409 the operator is self-adjoint by the Jordan axiom::
1410
1411 sage: set_random_seed()
1412 sage: x = ComplexHermitianEJA.random_instance().random_element()
1413 sage: x.operator().matrix().is_symmetric()
1414 True
1415
1416 """
1417
1418 @classmethod
1419 def _denormalized_basis(cls, n, field):
1420 """
1421 Returns a basis for the space of complex Hermitian n-by-n matrices.
1422
1423 Why do we embed these? Basically, because all of numerical linear
1424 algebra assumes that you're working with vectors consisting of `n`
1425 entries from a field and scalars from the same field. There's no way
1426 to tell SageMath that (for example) the vectors contain complex
1427 numbers, while the scalar field is real.
1428
1429 SETUP::
1430
1431 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1432
1433 TESTS::
1434
1435 sage: set_random_seed()
1436 sage: n = ZZ.random_element(1,5)
1437 sage: field = QuadraticField(2, 'sqrt2')
1438 sage: B = ComplexHermitianEJA._denormalized_basis(n, field)
1439 sage: all( M.is_symmetric() for M in B)
1440 True
1441
1442 """
1443 R = PolynomialRing(field, 'z')
1444 z = R.gen()
1445 F = field.extension(z**2 + 1, 'I')
1446 I = F.gen()
1447
1448 # This is like the symmetric case, but we need to be careful:
1449 #
1450 # * We want conjugate-symmetry, not just symmetry.
1451 # * The diagonal will (as a result) be real.
1452 #
1453 S = []
1454 for i in xrange(n):
1455 for j in xrange(i+1):
1456 Eij = matrix(F, n, lambda k,l: k==i and l==j)
1457 if i == j:
1458 Sij = cls.real_embed(Eij)
1459 S.append(Sij)
1460 else:
1461 # The second one has a minus because it's conjugated.
1462 Sij_real = cls.real_embed(Eij + Eij.transpose())
1463 S.append(Sij_real)
1464 Sij_imag = cls.real_embed(I*Eij - I*Eij.transpose())
1465 S.append(Sij_imag)
1466
1467 # Since we embedded these, we can drop back to the "field" that we
1468 # started with instead of the complex extension "F".
1469 return ( s.change_ring(field) for s in S )
1470
1471
1472 def __init__(self, n, field=QQ, **kwargs):
1473 basis = self._denormalized_basis(n,field)
1474 super(ComplexHermitianEJA,self).__init__(field, basis, n, **kwargs)
1475
1476
1477 class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
1478 @staticmethod
1479 def real_embed(M):
1480 """
1481 Embed the n-by-n quaternion matrix ``M`` into the space of real
1482 matrices of size 4n-by-4n by first sending each quaternion entry `z
1483 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
1484 c+di],[-c + di, a-bi]]`, and then embedding those into a real
1485 matrix.
1486
1487 SETUP::
1488
1489 sage: from mjo.eja.eja_algebra import \
1490 ....: QuaternionMatrixEuclideanJordanAlgebra
1491
1492 EXAMPLES::
1493
1494 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1495 sage: i,j,k = Q.gens()
1496 sage: x = 1 + 2*i + 3*j + 4*k
1497 sage: M = matrix(Q, 1, [[x]])
1498 sage: QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1499 [ 1 2 3 4]
1500 [-2 1 -4 3]
1501 [-3 4 1 -2]
1502 [-4 -3 2 1]
1503
1504 Embedding is a homomorphism (isomorphism, in fact)::
1505
1506 sage: set_random_seed()
1507 sage: n_max = QuaternionMatrixEuclideanJordanAlgebra._max_test_case_size()
1508 sage: n = ZZ.random_element(n_max)
1509 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1510 sage: X = random_matrix(Q, n)
1511 sage: Y = random_matrix(Q, n)
1512 sage: Xe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X)
1513 sage: Ye = QuaternionMatrixEuclideanJordanAlgebra.real_embed(Y)
1514 sage: XYe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1515 sage: Xe*Ye == XYe
1516 True
1517
1518 """
1519 quaternions = M.base_ring()
1520 n = M.nrows()
1521 if M.ncols() != n:
1522 raise ValueError("the matrix 'M' must be square")
1523
1524 F = QuadraticField(-1, 'i')
1525 i = F.gen()
1526
1527 blocks = []
1528 for z in M.list():
1529 t = z.coefficient_tuple()
1530 a = t[0]
1531 b = t[1]
1532 c = t[2]
1533 d = t[3]
1534 cplxM = matrix(F, 2, [[ a + b*i, c + d*i],
1535 [-c + d*i, a - b*i]])
1536 realM = ComplexMatrixEuclideanJordanAlgebra.real_embed(cplxM)
1537 blocks.append(realM)
1538
1539 # We should have real entries by now, so use the realest field
1540 # we've got for the return value.
1541 return matrix.block(quaternions.base_ring(), n, blocks)
1542
1543
1544
1545 @staticmethod
1546 def real_unembed(M):
1547 """
1548 The inverse of _embed_quaternion_matrix().
1549
1550 SETUP::
1551
1552 sage: from mjo.eja.eja_algebra import \
1553 ....: QuaternionMatrixEuclideanJordanAlgebra
1554
1555 EXAMPLES::
1556
1557 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
1558 ....: [-2, 1, -4, 3],
1559 ....: [-3, 4, 1, -2],
1560 ....: [-4, -3, 2, 1]])
1561 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(M)
1562 [1 + 2*i + 3*j + 4*k]
1563
1564 TESTS:
1565
1566 Unembedding is the inverse of embedding::
1567
1568 sage: set_random_seed()
1569 sage: Q = QuaternionAlgebra(QQ, -1, -1)
1570 sage: M = random_matrix(Q, 3)
1571 sage: Me = QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1572 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1573 True
1574
1575 """
1576 n = ZZ(M.nrows())
1577 if M.ncols() != n:
1578 raise ValueError("the matrix 'M' must be square")
1579 if not n.mod(4).is_zero():
1580 raise ValueError("the matrix 'M' must be a quaternion embedding")
1581
1582 # Use the base ring of the matrix to ensure that its entries can be
1583 # multiplied by elements of the quaternion algebra.
1584 field = M.base_ring()
1585 Q = QuaternionAlgebra(field,-1,-1)
1586 i,j,k = Q.gens()
1587
1588 # Go top-left to bottom-right (reading order), converting every
1589 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
1590 # quaternion block.
1591 elements = []
1592 for l in xrange(n/4):
1593 for m in xrange(n/4):
1594 submat = ComplexMatrixEuclideanJordanAlgebra.real_unembed(
1595 M[4*l:4*l+4,4*m:4*m+4] )
1596 if submat[0,0] != submat[1,1].conjugate():
1597 raise ValueError('bad on-diagonal submatrix')
1598 if submat[0,1] != -submat[1,0].conjugate():
1599 raise ValueError('bad off-diagonal submatrix')
1600 z = submat[0,0].vector()[0] # real part
1601 z += submat[0,0].vector()[1]*i # imag part
1602 z += submat[0,1].vector()[0]*j # real part
1603 z += submat[0,1].vector()[1]*k # imag part
1604 elements.append(z)
1605
1606 return matrix(Q, n/4, elements)
1607
1608
1609 @classmethod
1610 def natural_inner_product(cls,X,Y):
1611 """
1612 Compute a natural inner product in this algebra directly from
1613 its real embedding.
1614
1615 SETUP::
1616
1617 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1618
1619 TESTS:
1620
1621 This gives the same answer as the slow, default method implemented
1622 in :class:`MatrixEuclideanJordanAlgebra`::
1623
1624 sage: set_random_seed()
1625 sage: J = QuaternionHermitianEJA.random_instance()
1626 sage: x,y = J.random_elements(2)
1627 sage: Xe = x.natural_representation()
1628 sage: Ye = y.natural_representation()
1629 sage: X = QuaternionHermitianEJA.real_unembed(Xe)
1630 sage: Y = QuaternionHermitianEJA.real_unembed(Ye)
1631 sage: expected = (X*Y).trace().coefficient_tuple()[0]
1632 sage: actual = QuaternionHermitianEJA.natural_inner_product(Xe,Ye)
1633 sage: actual == expected
1634 True
1635
1636 """
1637 return RealMatrixEuclideanJordanAlgebra.natural_inner_product(X,Y)/4
1638
1639
1640 class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra,
1641 KnownRankEJA):
1642 """
1643 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
1644 matrices, the usual symmetric Jordan product, and the
1645 real-part-of-trace inner product. It has dimension `2n^2 - n` over
1646 the reals.
1647
1648 SETUP::
1649
1650 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1651
1652 EXAMPLES:
1653
1654 In theory, our "field" can be any subfield of the reals::
1655
1656 sage: QuaternionHermitianEJA(2, AA)
1657 Euclidean Jordan algebra of dimension 6 over Algebraic Real Field
1658 sage: QuaternionHermitianEJA(2, RR)
1659 Euclidean Jordan algebra of dimension 6 over Real Field with
1660 53 bits of precision
1661
1662 TESTS:
1663
1664 The dimension of this algebra is `2*n^2 - n`::
1665
1666 sage: set_random_seed()
1667 sage: n_max = QuaternionHermitianEJA._max_test_case_size()
1668 sage: n = ZZ.random_element(1, n_max)
1669 sage: J = QuaternionHermitianEJA(n)
1670 sage: J.dimension() == 2*(n^2) - n
1671 True
1672
1673 The Jordan multiplication is what we think it is::
1674
1675 sage: set_random_seed()
1676 sage: J = QuaternionHermitianEJA.random_instance()
1677 sage: x,y = J.random_elements(2)
1678 sage: actual = (x*y).natural_representation()
1679 sage: X = x.natural_representation()
1680 sage: Y = y.natural_representation()
1681 sage: expected = (X*Y + Y*X)/2
1682 sage: actual == expected
1683 True
1684 sage: J(expected) == x*y
1685 True
1686
1687 We can change the generator prefix::
1688
1689 sage: QuaternionHermitianEJA(2, prefix='a').gens()
1690 (a0, a1, a2, a3, a4, a5)
1691
1692 Our natural basis is normalized with respect to the natural inner
1693 product unless we specify otherwise::
1694
1695 sage: set_random_seed()
1696 sage: J = QuaternionHermitianEJA.random_instance()
1697 sage: all( b.norm() == 1 for b in J.gens() )
1698 True
1699
1700 Since our natural basis is normalized with respect to the natural
1701 inner product, and since we know that this algebra is an EJA, any
1702 left-multiplication operator's matrix will be symmetric because
1703 natural->EJA basis representation is an isometry and within the EJA
1704 the operator is self-adjoint by the Jordan axiom::
1705
1706 sage: set_random_seed()
1707 sage: x = QuaternionHermitianEJA.random_instance().random_element()
1708 sage: x.operator().matrix().is_symmetric()
1709 True
1710
1711 """
1712 @classmethod
1713 def _denormalized_basis(cls, n, field):
1714 """
1715 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
1716
1717 Why do we embed these? Basically, because all of numerical
1718 linear algebra assumes that you're working with vectors consisting
1719 of `n` entries from a field and scalars from the same field. There's
1720 no way to tell SageMath that (for example) the vectors contain
1721 complex numbers, while the scalar field is real.
1722
1723 SETUP::
1724
1725 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1726
1727 TESTS::
1728
1729 sage: set_random_seed()
1730 sage: n = ZZ.random_element(1,5)
1731 sage: B = QuaternionHermitianEJA._denormalized_basis(n,QQ)
1732 sage: all( M.is_symmetric() for M in B )
1733 True
1734
1735 """
1736 Q = QuaternionAlgebra(QQ,-1,-1)
1737 I,J,K = Q.gens()
1738
1739 # This is like the symmetric case, but we need to be careful:
1740 #
1741 # * We want conjugate-symmetry, not just symmetry.
1742 # * The diagonal will (as a result) be real.
1743 #
1744 S = []
1745 for i in xrange(n):
1746 for j in xrange(i+1):
1747 Eij = matrix(Q, n, lambda k,l: k==i and l==j)
1748 if i == j:
1749 Sij = cls.real_embed(Eij)
1750 S.append(Sij)
1751 else:
1752 # The second, third, and fourth ones have a minus
1753 # because they're conjugated.
1754 Sij_real = cls.real_embed(Eij + Eij.transpose())
1755 S.append(Sij_real)
1756 Sij_I = cls.real_embed(I*Eij - I*Eij.transpose())
1757 S.append(Sij_I)
1758 Sij_J = cls.real_embed(J*Eij - J*Eij.transpose())
1759 S.append(Sij_J)
1760 Sij_K = cls.real_embed(K*Eij - K*Eij.transpose())
1761 S.append(Sij_K)
1762
1763 # Since we embedded these, we can drop back to the "field" that we
1764 # started with instead of the quaternion algebra "Q".
1765 return ( s.change_ring(field) for s in S )
1766
1767
1768 def __init__(self, n, field=QQ, **kwargs):
1769 basis = self._denormalized_basis(n,field)
1770 super(QuaternionHermitianEJA,self).__init__(field, basis, n, **kwargs)
1771
1772
1773 class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra, KnownRankEJA):
1774 """
1775 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
1776 with the usual inner product and jordan product ``x*y =
1777 (<x_bar,y_bar>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
1778 the reals.
1779
1780 SETUP::
1781
1782 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1783
1784 EXAMPLES:
1785
1786 This multiplication table can be verified by hand::
1787
1788 sage: J = JordanSpinEJA(4)
1789 sage: e0,e1,e2,e3 = J.gens()
1790 sage: e0*e0
1791 e0
1792 sage: e0*e1
1793 e1
1794 sage: e0*e2
1795 e2
1796 sage: e0*e3
1797 e3
1798 sage: e1*e2
1799 0
1800 sage: e1*e3
1801 0
1802 sage: e2*e3
1803 0
1804
1805 We can change the generator prefix::
1806
1807 sage: JordanSpinEJA(2, prefix='B').gens()
1808 (B0, B1)
1809
1810 """
1811 def __init__(self, n, field=QQ, **kwargs):
1812 V = VectorSpace(field, n)
1813 mult_table = [[V.zero() for j in xrange(n)] for i in xrange(n)]
1814 for i in xrange(n):
1815 for j in xrange(n):
1816 x = V.gen(i)
1817 y = V.gen(j)
1818 x0 = x[0]
1819 xbar = x[1:]
1820 y0 = y[0]
1821 ybar = y[1:]
1822 # z = x*y
1823 z0 = x.inner_product(y)
1824 zbar = y0*xbar + x0*ybar
1825 z = V([z0] + zbar.list())
1826 mult_table[i][j] = z
1827
1828 # The rank of the spin algebra is two, unless we're in a
1829 # one-dimensional ambient space (because the rank is bounded by
1830 # the ambient dimension).
1831 fdeja = super(JordanSpinEJA, self)
1832 return fdeja.__init__(field, mult_table, rank=min(n,2), **kwargs)
1833
1834 def inner_product(self, x, y):
1835 """
1836 Faster to reimplement than to use natural representations.
1837
1838 SETUP::
1839
1840 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1841
1842 TESTS:
1843
1844 Ensure that this is the usual inner product for the algebras
1845 over `R^n`::
1846
1847 sage: set_random_seed()
1848 sage: J = JordanSpinEJA.random_instance()
1849 sage: x,y = J.random_elements(2)
1850 sage: X = x.natural_representation()
1851 sage: Y = y.natural_representation()
1852 sage: x.inner_product(y) == J.natural_inner_product(X,Y)
1853 True
1854
1855 """
1856 return x.to_vector().inner_product(y.to_vector())