]> gitweb.michael.orlitzky.com - sage.d.git/blob - mjo/eja/eja_algebra.py
eja: add multiplication_table() for EJAs.
[sage.d.git] / mjo / eja / 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 sage.algebras.quatalg.quaternion_algebra import QuaternionAlgebra
9 from sage.categories.finite_dimensional_algebras_with_basis import FiniteDimensionalAlgebrasWithBasis
10 from sage.combinat.free_module import CombinatorialFreeModule
11 from sage.matrix.constructor import matrix
12 from sage.misc.cachefunc import cached_method
13 from sage.misc.prandom import choice
14 from sage.modules.free_module import VectorSpace
15 from sage.rings.integer_ring import ZZ
16 from sage.rings.number_field.number_field import QuadraticField
17 from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
18 from sage.rings.rational_field import QQ
19 from sage.structure.element import is_Matrix
20
21 from mjo.eja.eja_element import FiniteDimensionalEuclideanJordanAlgebraElement
22 from mjo.eja.eja_utils import _mat2vec
23
24 class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
25 def __init__(self,
26 field,
27 mult_table,
28 rank,
29 prefix='e',
30 category=None,
31 natural_basis=None):
32 """
33 SETUP::
34
35 sage: from mjo.eja.eja_algebra import random_eja
36
37 EXAMPLES:
38
39 By definition, Jordan multiplication commutes::
40
41 sage: set_random_seed()
42 sage: J = random_eja()
43 sage: x = J.random_element()
44 sage: y = J.random_element()
45 sage: x*y == y*x
46 True
47
48 """
49 self._rank = rank
50 self._natural_basis = natural_basis
51
52 if category is None:
53 category = FiniteDimensionalAlgebrasWithBasis(field).Unital()
54 fda = super(FiniteDimensionalEuclideanJordanAlgebra, self)
55 fda.__init__(field,
56 range(len(mult_table)),
57 prefix=prefix,
58 category=category)
59 self.print_options(bracket='')
60
61 # The multiplication table we're given is necessarily in terms
62 # of vectors, because we don't have an algebra yet for
63 # anything to be an element of. However, it's faster in the
64 # long run to have the multiplication table be in terms of
65 # algebra elements. We do this after calling the superclass
66 # constructor so that from_vector() knows what to do.
67 self._multiplication_table = [ map(lambda x: self.from_vector(x), ls)
68 for ls in mult_table ]
69
70
71 def _element_constructor_(self, elt):
72 """
73 Construct an element of this algebra from its natural
74 representation.
75
76 This gets called only after the parent element _call_ method
77 fails to find a coercion for the argument.
78
79 SETUP::
80
81 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
82 ....: RealCartesianProductEJA,
83 ....: RealSymmetricEJA)
84
85 EXAMPLES:
86
87 The identity in `S^n` is converted to the identity in the EJA::
88
89 sage: J = RealSymmetricEJA(3)
90 sage: I = matrix.identity(QQ,3)
91 sage: J(I) == J.one()
92 True
93
94 This skew-symmetric matrix can't be represented in the EJA::
95
96 sage: J = RealSymmetricEJA(3)
97 sage: A = matrix(QQ,3, lambda i,j: i-j)
98 sage: J(A)
99 Traceback (most recent call last):
100 ...
101 ArithmeticError: vector is not in free module
102
103 TESTS:
104
105 Ensure that we can convert any element of the two non-matrix
106 simple algebras (whose natural representations are their usual
107 vector representations) back and forth faithfully::
108
109 sage: set_random_seed()
110 sage: J = RealCartesianProductEJA(5)
111 sage: x = J.random_element()
112 sage: J(x.to_vector().column()) == x
113 True
114 sage: J = JordanSpinEJA(5)
115 sage: x = J.random_element()
116 sage: J(x.to_vector().column()) == x
117 True
118
119 """
120 natural_basis = self.natural_basis()
121 if elt not in natural_basis[0].matrix_space():
122 raise ValueError("not a naturally-represented algebra element")
123
124 # Thanks for nothing! Matrix spaces aren't vector
125 # spaces in Sage, so we have to figure out its
126 # natural-basis coordinates ourselves.
127 V = VectorSpace(elt.base_ring(), elt.nrows()*elt.ncols())
128 W = V.span_of_basis( _mat2vec(s) for s in natural_basis )
129 coords = W.coordinate_vector(_mat2vec(elt))
130 return self.from_vector(coords)
131
132
133 def _repr_(self):
134 """
135 Return a string representation of ``self``.
136
137 SETUP::
138
139 sage: from mjo.eja.eja_algebra import JordanSpinEJA
140
141 TESTS:
142
143 Ensure that it says what we think it says::
144
145 sage: JordanSpinEJA(2, field=QQ)
146 Euclidean Jordan algebra of degree 2 over Rational Field
147 sage: JordanSpinEJA(3, field=RDF)
148 Euclidean Jordan algebra of degree 3 over Real Double Field
149
150 """
151 # TODO: change this to say "dimension" and fix all the tests.
152 fmt = "Euclidean Jordan algebra of degree {} over {}"
153 return fmt.format(self.dimension(), self.base_ring())
154
155 def product_on_basis(self, i, j):
156 return self._multiplication_table[i][j]
157
158 def _a_regular_element(self):
159 """
160 Guess a regular element. Needed to compute the basis for our
161 characteristic polynomial coefficients.
162
163 SETUP::
164
165 sage: from mjo.eja.eja_algebra import random_eja
166
167 TESTS:
168
169 Ensure that this hacky method succeeds for every algebra that we
170 know how to construct::
171
172 sage: set_random_seed()
173 sage: J = random_eja()
174 sage: J._a_regular_element().is_regular()
175 True
176
177 """
178 gs = self.gens()
179 z = self.sum( (i+1)*gs[i] for i in range(len(gs)) )
180 if not z.is_regular():
181 raise ValueError("don't know a regular element")
182 return z
183
184
185 @cached_method
186 def _charpoly_basis_space(self):
187 """
188 Return the vector space spanned by the basis used in our
189 characteristic polynomial coefficients. This is used not only to
190 compute those coefficients, but also any time we need to
191 evaluate the coefficients (like when we compute the trace or
192 determinant).
193 """
194 z = self._a_regular_element()
195 V = self.vector_space()
196 V1 = V.span_of_basis( (z**k).to_vector() for k in range(self.rank()) )
197 b = (V1.basis() + V1.complement().basis())
198 return V.span_of_basis(b)
199
200
201 @cached_method
202 def _charpoly_coeff(self, i):
203 """
204 Return the coefficient polynomial "a_{i}" of this algebra's
205 general characteristic polynomial.
206
207 Having this be a separate cached method lets us compute and
208 store the trace/determinant (a_{r-1} and a_{0} respectively)
209 separate from the entire characteristic polynomial.
210 """
211 (A_of_x, x, xr, detA) = self._charpoly_matrix_system()
212 R = A_of_x.base_ring()
213 if i >= self.rank():
214 # Guaranteed by theory
215 return R.zero()
216
217 # Danger: the in-place modification is done for performance
218 # reasons (reconstructing a matrix with huge polynomial
219 # entries is slow), but I don't know how cached_method works,
220 # so it's highly possible that we're modifying some global
221 # list variable by reference, here. In other words, you
222 # probably shouldn't call this method twice on the same
223 # algebra, at the same time, in two threads
224 Ai_orig = A_of_x.column(i)
225 A_of_x.set_column(i,xr)
226 numerator = A_of_x.det()
227 A_of_x.set_column(i,Ai_orig)
228
229 # We're relying on the theory here to ensure that each a_i is
230 # indeed back in R, and the added negative signs are to make
231 # the whole charpoly expression sum to zero.
232 return R(-numerator/detA)
233
234
235 @cached_method
236 def _charpoly_matrix_system(self):
237 """
238 Compute the matrix whose entries A_ij are polynomials in
239 X1,...,XN, the vector ``x`` of variables X1,...,XN, the vector
240 corresponding to `x^r` and the determinent of the matrix A =
241 [A_ij]. In other words, all of the fixed (cachable) data needed
242 to compute the coefficients of the characteristic polynomial.
243 """
244 r = self.rank()
245 n = self.dimension()
246
247 # Construct a new algebra over a multivariate polynomial ring...
248 names = tuple('X' + str(i) for i in range(1,n+1))
249 R = PolynomialRing(self.base_ring(), names)
250 # Hack around the fact that our multiplication table is in terms of
251 # algebra elements but the constructor wants it in terms of vectors.
252 vmt = [ tuple(map(lambda x: x.to_vector(), ls))
253 for ls in self._multiplication_table ]
254 J = FiniteDimensionalEuclideanJordanAlgebra(R, tuple(vmt), r)
255
256 idmat = matrix.identity(J.base_ring(), n)
257
258 W = self._charpoly_basis_space()
259 W = W.change_ring(R.fraction_field())
260
261 # Starting with the standard coordinates x = (X1,X2,...,Xn)
262 # and then converting the entries to W-coordinates allows us
263 # to pass in the standard coordinates to the charpoly and get
264 # back the right answer. Specifically, with x = (X1,X2,...,Xn),
265 # we have
266 #
267 # W.coordinates(x^2) eval'd at (standard z-coords)
268 # =
269 # W-coords of (z^2)
270 # =
271 # W-coords of (standard coords of x^2 eval'd at std-coords of z)
272 #
273 # We want the middle equivalent thing in our matrix, but use
274 # the first equivalent thing instead so that we can pass in
275 # standard coordinates.
276 x = J.from_vector(W(R.gens()))
277
278 # Handle the zeroth power separately, because computing
279 # the unit element in J is mathematically suspect.
280 x0 = W.coordinate_vector(self.one().to_vector())
281 l1 = [ x0.column() ]
282 l1 += [ W.coordinate_vector((x**k).to_vector()).column()
283 for k in range(1,r) ]
284 l2 = [idmat.column(k-1).column() for k in range(r+1, n+1)]
285 A_of_x = matrix.block(R, 1, n, (l1 + l2))
286 xr = W.coordinate_vector((x**r).to_vector())
287 return (A_of_x, x, xr, A_of_x.det())
288
289
290 @cached_method
291 def characteristic_polynomial(self):
292 """
293 Return a characteristic polynomial that works for all elements
294 of this algebra.
295
296 The resulting polynomial has `n+1` variables, where `n` is the
297 dimension of this algebra. The first `n` variables correspond to
298 the coordinates of an algebra element: when evaluated at the
299 coordinates of an algebra element with respect to a certain
300 basis, the result is a univariate polynomial (in the one
301 remaining variable ``t``), namely the characteristic polynomial
302 of that element.
303
304 SETUP::
305
306 sage: from mjo.eja.eja_algebra import JordanSpinEJA
307
308 EXAMPLES:
309
310 The characteristic polynomial in the spin algebra is given in
311 Alizadeh, Example 11.11::
312
313 sage: J = JordanSpinEJA(3)
314 sage: p = J.characteristic_polynomial(); p
315 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
316 sage: xvec = J.one().to_vector()
317 sage: p(*xvec)
318 t^2 - 2*t + 1
319
320 """
321 r = self.rank()
322 n = self.dimension()
323
324 # The list of coefficient polynomials a_1, a_2, ..., a_n.
325 a = [ self._charpoly_coeff(i) for i in range(n) ]
326
327 # We go to a bit of trouble here to reorder the
328 # indeterminates, so that it's easier to evaluate the
329 # characteristic polynomial at x's coordinates and get back
330 # something in terms of t, which is what we want.
331 R = a[0].parent()
332 S = PolynomialRing(self.base_ring(),'t')
333 t = S.gen(0)
334 S = PolynomialRing(S, R.variable_names())
335 t = S(t)
336
337 # Note: all entries past the rth should be zero. The
338 # coefficient of the highest power (x^r) is 1, but it doesn't
339 # appear in the solution vector which contains coefficients
340 # for the other powers (to make them sum to x^r).
341 if (r < n):
342 a[r] = 1 # corresponds to x^r
343 else:
344 # When the rank is equal to the dimension, trying to
345 # assign a[r] goes out-of-bounds.
346 a.append(1) # corresponds to x^r
347
348 return sum( a[k]*(t**k) for k in range(len(a)) )
349
350
351 def inner_product(self, x, y):
352 """
353 The inner product associated with this Euclidean Jordan algebra.
354
355 Defaults to the trace inner product, but can be overridden by
356 subclasses if they are sure that the necessary properties are
357 satisfied.
358
359 SETUP::
360
361 sage: from mjo.eja.eja_algebra import random_eja
362
363 EXAMPLES:
364
365 The inner product must satisfy its axiom for this algebra to truly
366 be a Euclidean Jordan Algebra::
367
368 sage: set_random_seed()
369 sage: J = random_eja()
370 sage: x = J.random_element()
371 sage: y = J.random_element()
372 sage: z = J.random_element()
373 sage: (x*y).inner_product(z) == y.inner_product(x*z)
374 True
375
376 """
377 if (not x in self) or (not y in self):
378 raise TypeError("arguments must live in this algebra")
379 return x.trace_inner_product(y)
380
381
382 def multiplication_table(self):
383 """
384 Return a readable matrix representation of this algebra's
385 multiplication table. The (i,j)th entry in the matrix contains
386 the product of the ith basis element with the jth.
387
388 This is not extraordinarily useful, but it overrides a superclass
389 method that would otherwise just crash and complain about the
390 algebra being infinite.
391
392 EXAMPLES::
393
394 sage: J = RealCartesianProductEJA(3)
395 sage: J.multiplication_table()
396 [e0 0 0]
397 [ 0 e1 0]
398 [ 0 0 e2]
399
400 ::
401
402 sage: J = JordanSpinEJA(3)
403 sage: J.multiplication_table()
404 [e0 e1 e2]
405 [e1 e0 0]
406 [e2 0 e0]
407
408 """
409 return matrix(self._multiplication_table)
410
411
412 def natural_basis(self):
413 """
414 Return a more-natural representation of this algebra's basis.
415
416 Every finite-dimensional Euclidean Jordan Algebra is a direct
417 sum of five simple algebras, four of which comprise Hermitian
418 matrices. This method returns the original "natural" basis
419 for our underlying vector space. (Typically, the natural basis
420 is used to construct the multiplication table in the first place.)
421
422 Note that this will always return a matrix. The standard basis
423 in `R^n` will be returned as `n`-by-`1` column matrices.
424
425 SETUP::
426
427 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
428 ....: RealSymmetricEJA)
429
430 EXAMPLES::
431
432 sage: J = RealSymmetricEJA(2)
433 sage: J.basis()
434 Finite family {0: e0, 1: e1, 2: e2}
435 sage: J.natural_basis()
436 (
437 [1 0] [0 1] [0 0]
438 [0 0], [1 0], [0 1]
439 )
440
441 ::
442
443 sage: J = JordanSpinEJA(2)
444 sage: J.basis()
445 Finite family {0: e0, 1: e1}
446 sage: J.natural_basis()
447 (
448 [1] [0]
449 [0], [1]
450 )
451
452 """
453 if self._natural_basis is None:
454 return tuple( b.to_vector().column() for b in self.basis() )
455 else:
456 return self._natural_basis
457
458
459 @cached_method
460 def one(self):
461 """
462 Return the unit element of this algebra.
463
464 SETUP::
465
466 sage: from mjo.eja.eja_algebra import (RealCartesianProductEJA,
467 ....: random_eja)
468
469 EXAMPLES::
470
471 sage: J = RealCartesianProductEJA(5)
472 sage: J.one()
473 e0 + e1 + e2 + e3 + e4
474
475 TESTS::
476
477 The identity element acts like the identity::
478
479 sage: set_random_seed()
480 sage: J = random_eja()
481 sage: x = J.random_element()
482 sage: J.one()*x == x and x*J.one() == x
483 True
484
485 The matrix of the unit element's operator is the identity::
486
487 sage: set_random_seed()
488 sage: J = random_eja()
489 sage: actual = J.one().operator().matrix()
490 sage: expected = matrix.identity(J.base_ring(), J.dimension())
491 sage: actual == expected
492 True
493
494 """
495 # We can brute-force compute the matrices of the operators
496 # that correspond to the basis elements of this algebra.
497 # If some linear combination of those basis elements is the
498 # algebra identity, then the same linear combination of
499 # their matrices has to be the identity matrix.
500 #
501 # Of course, matrices aren't vectors in sage, so we have to
502 # appeal to the "long vectors" isometry.
503 oper_vecs = [ _mat2vec(g.operator().matrix()) for g in self.gens() ]
504
505 # Now we use basis linear algebra to find the coefficients,
506 # of the matrices-as-vectors-linear-combination, which should
507 # work for the original algebra basis too.
508 A = matrix.column(self.base_ring(), oper_vecs)
509
510 # We used the isometry on the left-hand side already, but we
511 # still need to do it for the right-hand side. Recall that we
512 # wanted something that summed to the identity matrix.
513 b = _mat2vec( matrix.identity(self.base_ring(), self.dimension()) )
514
515 # Now if there's an identity element in the algebra, this should work.
516 coeffs = A.solve_right(b)
517 return self.linear_combination(zip(self.gens(), coeffs))
518
519
520 def rank(self):
521 """
522 Return the rank of this EJA.
523
524 ALGORITHM:
525
526 The author knows of no algorithm to compute the rank of an EJA
527 where only the multiplication table is known. In lieu of one, we
528 require the rank to be specified when the algebra is created,
529 and simply pass along that number here.
530
531 SETUP::
532
533 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
534 ....: RealSymmetricEJA,
535 ....: ComplexHermitianEJA,
536 ....: QuaternionHermitianEJA,
537 ....: random_eja)
538
539 EXAMPLES:
540
541 The rank of the Jordan spin algebra is always two::
542
543 sage: JordanSpinEJA(2).rank()
544 2
545 sage: JordanSpinEJA(3).rank()
546 2
547 sage: JordanSpinEJA(4).rank()
548 2
549
550 The rank of the `n`-by-`n` Hermitian real, complex, or
551 quaternion matrices is `n`::
552
553 sage: RealSymmetricEJA(2).rank()
554 2
555 sage: ComplexHermitianEJA(2).rank()
556 2
557 sage: QuaternionHermitianEJA(2).rank()
558 2
559 sage: RealSymmetricEJA(5).rank()
560 5
561 sage: ComplexHermitianEJA(5).rank()
562 5
563 sage: QuaternionHermitianEJA(5).rank()
564 5
565
566 TESTS:
567
568 Ensure that every EJA that we know how to construct has a
569 positive integer rank::
570
571 sage: set_random_seed()
572 sage: r = random_eja().rank()
573 sage: r in ZZ and r > 0
574 True
575
576 """
577 return self._rank
578
579
580 def vector_space(self):
581 """
582 Return the vector space that underlies this algebra.
583
584 SETUP::
585
586 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
587
588 EXAMPLES::
589
590 sage: J = RealSymmetricEJA(2)
591 sage: J.vector_space()
592 Vector space of dimension 3 over Rational Field
593
594 """
595 return self.zero().to_vector().parent().ambient_vector_space()
596
597
598 Element = FiniteDimensionalEuclideanJordanAlgebraElement
599
600
601 class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra):
602 """
603 Return the Euclidean Jordan Algebra corresponding to the set
604 `R^n` under the Hadamard product.
605
606 Note: this is nothing more than the Cartesian product of ``n``
607 copies of the spin algebra. Once Cartesian product algebras
608 are implemented, this can go.
609
610 SETUP::
611
612 sage: from mjo.eja.eja_algebra import RealCartesianProductEJA
613
614 EXAMPLES:
615
616 This multiplication table can be verified by hand::
617
618 sage: J = RealCartesianProductEJA(3)
619 sage: e0,e1,e2 = J.gens()
620 sage: e0*e0
621 e0
622 sage: e0*e1
623 0
624 sage: e0*e2
625 0
626 sage: e1*e1
627 e1
628 sage: e1*e2
629 0
630 sage: e2*e2
631 e2
632
633 """
634 def __init__(self, n, field=QQ):
635 V = VectorSpace(field, n)
636 mult_table = [ [ V.basis()[i]*(i == j) for i in range(n) ]
637 for j in range(n) ]
638
639 fdeja = super(RealCartesianProductEJA, self)
640 return fdeja.__init__(field, mult_table, rank=n)
641
642 def inner_product(self, x, y):
643 return _usual_ip(x,y)
644
645
646 def random_eja():
647 """
648 Return a "random" finite-dimensional Euclidean Jordan Algebra.
649
650 ALGORITHM:
651
652 For now, we choose a random natural number ``n`` (greater than zero)
653 and then give you back one of the following:
654
655 * The cartesian product of the rational numbers ``n`` times; this is
656 ``QQ^n`` with the Hadamard product.
657
658 * The Jordan spin algebra on ``QQ^n``.
659
660 * The ``n``-by-``n`` rational symmetric matrices with the symmetric
661 product.
662
663 * The ``n``-by-``n`` complex-rational Hermitian matrices embedded
664 in the space of ``2n``-by-``2n`` real symmetric matrices.
665
666 * The ``n``-by-``n`` quaternion-rational Hermitian matrices embedded
667 in the space of ``4n``-by-``4n`` real symmetric matrices.
668
669 Later this might be extended to return Cartesian products of the
670 EJAs above.
671
672 SETUP::
673
674 sage: from mjo.eja.eja_algebra import random_eja
675
676 TESTS::
677
678 sage: random_eja()
679 Euclidean Jordan algebra of degree...
680
681 """
682
683 # The max_n component lets us choose different upper bounds on the
684 # value "n" that gets passed to the constructor. This is needed
685 # because e.g. R^{10} is reasonable to test, while the Hermitian
686 # 10-by-10 quaternion matrices are not.
687 (constructor, max_n) = choice([(RealCartesianProductEJA, 6),
688 (JordanSpinEJA, 6),
689 (RealSymmetricEJA, 5),
690 (ComplexHermitianEJA, 4),
691 (QuaternionHermitianEJA, 3)])
692 n = ZZ.random_element(1, max_n)
693 return constructor(n, field=QQ)
694
695
696
697 def _real_symmetric_basis(n, field=QQ):
698 """
699 Return a basis for the space of real symmetric n-by-n matrices.
700 """
701 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
702 # coordinates.
703 S = []
704 for i in xrange(n):
705 for j in xrange(i+1):
706 Eij = matrix(field, n, lambda k,l: k==i and l==j)
707 if i == j:
708 Sij = Eij
709 else:
710 # Beware, orthogonal but not normalized!
711 Sij = Eij + Eij.transpose()
712 S.append(Sij)
713 return tuple(S)
714
715
716 def _complex_hermitian_basis(n, field=QQ):
717 """
718 Returns a basis for the space of complex Hermitian n-by-n matrices.
719
720 SETUP::
721
722 sage: from mjo.eja.eja_algebra import _complex_hermitian_basis
723
724 TESTS::
725
726 sage: set_random_seed()
727 sage: n = ZZ.random_element(1,5)
728 sage: all( M.is_symmetric() for M in _complex_hermitian_basis(n) )
729 True
730
731 """
732 F = QuadraticField(-1, 'I')
733 I = F.gen()
734
735 # This is like the symmetric case, but we need to be careful:
736 #
737 # * We want conjugate-symmetry, not just symmetry.
738 # * The diagonal will (as a result) be real.
739 #
740 S = []
741 for i in xrange(n):
742 for j in xrange(i+1):
743 Eij = matrix(field, n, lambda k,l: k==i and l==j)
744 if i == j:
745 Sij = _embed_complex_matrix(Eij)
746 S.append(Sij)
747 else:
748 # Beware, orthogonal but not normalized! The second one
749 # has a minus because it's conjugated.
750 Sij_real = _embed_complex_matrix(Eij + Eij.transpose())
751 S.append(Sij_real)
752 Sij_imag = _embed_complex_matrix(I*Eij - I*Eij.transpose())
753 S.append(Sij_imag)
754 return tuple(S)
755
756
757 def _quaternion_hermitian_basis(n, field=QQ):
758 """
759 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
760
761 SETUP::
762
763 sage: from mjo.eja.eja_algebra import _quaternion_hermitian_basis
764
765 TESTS::
766
767 sage: set_random_seed()
768 sage: n = ZZ.random_element(1,5)
769 sage: all( M.is_symmetric() for M in _quaternion_hermitian_basis(n) )
770 True
771
772 """
773 Q = QuaternionAlgebra(QQ,-1,-1)
774 I,J,K = Q.gens()
775
776 # This is like the symmetric case, but we need to be careful:
777 #
778 # * We want conjugate-symmetry, not just symmetry.
779 # * The diagonal will (as a result) be real.
780 #
781 S = []
782 for i in xrange(n):
783 for j in xrange(i+1):
784 Eij = matrix(Q, n, lambda k,l: k==i and l==j)
785 if i == j:
786 Sij = _embed_quaternion_matrix(Eij)
787 S.append(Sij)
788 else:
789 # Beware, orthogonal but not normalized! The second,
790 # third, and fourth ones have a minus because they're
791 # conjugated.
792 Sij_real = _embed_quaternion_matrix(Eij + Eij.transpose())
793 S.append(Sij_real)
794 Sij_I = _embed_quaternion_matrix(I*Eij - I*Eij.transpose())
795 S.append(Sij_I)
796 Sij_J = _embed_quaternion_matrix(J*Eij - J*Eij.transpose())
797 S.append(Sij_J)
798 Sij_K = _embed_quaternion_matrix(K*Eij - K*Eij.transpose())
799 S.append(Sij_K)
800 return tuple(S)
801
802
803
804 def _multiplication_table_from_matrix_basis(basis):
805 """
806 At least three of the five simple Euclidean Jordan algebras have the
807 symmetric multiplication (A,B) |-> (AB + BA)/2, where the
808 multiplication on the right is matrix multiplication. Given a basis
809 for the underlying matrix space, this function returns a
810 multiplication table (obtained by looping through the basis
811 elements) for an algebra of those matrices.
812 """
813 # In S^2, for example, we nominally have four coordinates even
814 # though the space is of dimension three only. The vector space V
815 # is supposed to hold the entire long vector, and the subspace W
816 # of V will be spanned by the vectors that arise from symmetric
817 # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
818 field = basis[0].base_ring()
819 dimension = basis[0].nrows()
820
821 V = VectorSpace(field, dimension**2)
822 W = V.span_of_basis( _mat2vec(s) for s in basis )
823 n = len(basis)
824 mult_table = [[W.zero() for i in range(n)] for j in range(n)]
825 for i in range(n):
826 for j in range(n):
827 mat_entry = (basis[i]*basis[j] + basis[j]*basis[i])/2
828 mult_table[i][j] = W.coordinate_vector(_mat2vec(mat_entry))
829
830 return mult_table
831
832
833 def _embed_complex_matrix(M):
834 """
835 Embed the n-by-n complex matrix ``M`` into the space of real
836 matrices of size 2n-by-2n via the map the sends each entry `z = a +
837 bi` to the block matrix ``[[a,b],[-b,a]]``.
838
839 SETUP::
840
841 sage: from mjo.eja.eja_algebra import _embed_complex_matrix
842
843 EXAMPLES::
844
845 sage: F = QuadraticField(-1,'i')
846 sage: x1 = F(4 - 2*i)
847 sage: x2 = F(1 + 2*i)
848 sage: x3 = F(-i)
849 sage: x4 = F(6)
850 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
851 sage: _embed_complex_matrix(M)
852 [ 4 -2| 1 2]
853 [ 2 4|-2 1]
854 [-----+-----]
855 [ 0 -1| 6 0]
856 [ 1 0| 0 6]
857
858 TESTS:
859
860 Embedding is a homomorphism (isomorphism, in fact)::
861
862 sage: set_random_seed()
863 sage: n = ZZ.random_element(5)
864 sage: F = QuadraticField(-1, 'i')
865 sage: X = random_matrix(F, n)
866 sage: Y = random_matrix(F, n)
867 sage: actual = _embed_complex_matrix(X) * _embed_complex_matrix(Y)
868 sage: expected = _embed_complex_matrix(X*Y)
869 sage: actual == expected
870 True
871
872 """
873 n = M.nrows()
874 if M.ncols() != n:
875 raise ValueError("the matrix 'M' must be square")
876 field = M.base_ring()
877 blocks = []
878 for z in M.list():
879 a = z.real()
880 b = z.imag()
881 blocks.append(matrix(field, 2, [[a,b],[-b,a]]))
882
883 # We can drop the imaginaries here.
884 return matrix.block(field.base_ring(), n, blocks)
885
886
887 def _unembed_complex_matrix(M):
888 """
889 The inverse of _embed_complex_matrix().
890
891 SETUP::
892
893 sage: from mjo.eja.eja_algebra import (_embed_complex_matrix,
894 ....: _unembed_complex_matrix)
895
896 EXAMPLES::
897
898 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
899 ....: [-2, 1, -4, 3],
900 ....: [ 9, 10, 11, 12],
901 ....: [-10, 9, -12, 11] ])
902 sage: _unembed_complex_matrix(A)
903 [ 2*i + 1 4*i + 3]
904 [ 10*i + 9 12*i + 11]
905
906 TESTS:
907
908 Unembedding is the inverse of embedding::
909
910 sage: set_random_seed()
911 sage: F = QuadraticField(-1, 'i')
912 sage: M = random_matrix(F, 3)
913 sage: _unembed_complex_matrix(_embed_complex_matrix(M)) == M
914 True
915
916 """
917 n = ZZ(M.nrows())
918 if M.ncols() != n:
919 raise ValueError("the matrix 'M' must be square")
920 if not n.mod(2).is_zero():
921 raise ValueError("the matrix 'M' must be a complex embedding")
922
923 F = QuadraticField(-1, 'i')
924 i = F.gen()
925
926 # Go top-left to bottom-right (reading order), converting every
927 # 2-by-2 block we see to a single complex element.
928 elements = []
929 for k in xrange(n/2):
930 for j in xrange(n/2):
931 submat = M[2*k:2*k+2,2*j:2*j+2]
932 if submat[0,0] != submat[1,1]:
933 raise ValueError('bad on-diagonal submatrix')
934 if submat[0,1] != -submat[1,0]:
935 raise ValueError('bad off-diagonal submatrix')
936 z = submat[0,0] + submat[0,1]*i
937 elements.append(z)
938
939 return matrix(F, n/2, elements)
940
941
942 def _embed_quaternion_matrix(M):
943 """
944 Embed the n-by-n quaternion matrix ``M`` into the space of real
945 matrices of size 4n-by-4n by first sending each quaternion entry
946 `z = a + bi + cj + dk` to the block-complex matrix
947 ``[[a + bi, c+di],[-c + di, a-bi]]`, and then embedding those into
948 a real matrix.
949
950 SETUP::
951
952 sage: from mjo.eja.eja_algebra import _embed_quaternion_matrix
953
954 EXAMPLES::
955
956 sage: Q = QuaternionAlgebra(QQ,-1,-1)
957 sage: i,j,k = Q.gens()
958 sage: x = 1 + 2*i + 3*j + 4*k
959 sage: M = matrix(Q, 1, [[x]])
960 sage: _embed_quaternion_matrix(M)
961 [ 1 2 3 4]
962 [-2 1 -4 3]
963 [-3 4 1 -2]
964 [-4 -3 2 1]
965
966 Embedding is a homomorphism (isomorphism, in fact)::
967
968 sage: set_random_seed()
969 sage: n = ZZ.random_element(5)
970 sage: Q = QuaternionAlgebra(QQ,-1,-1)
971 sage: X = random_matrix(Q, n)
972 sage: Y = random_matrix(Q, n)
973 sage: actual = _embed_quaternion_matrix(X)*_embed_quaternion_matrix(Y)
974 sage: expected = _embed_quaternion_matrix(X*Y)
975 sage: actual == expected
976 True
977
978 """
979 quaternions = M.base_ring()
980 n = M.nrows()
981 if M.ncols() != n:
982 raise ValueError("the matrix 'M' must be square")
983
984 F = QuadraticField(-1, 'i')
985 i = F.gen()
986
987 blocks = []
988 for z in M.list():
989 t = z.coefficient_tuple()
990 a = t[0]
991 b = t[1]
992 c = t[2]
993 d = t[3]
994 cplx_matrix = matrix(F, 2, [[ a + b*i, c + d*i],
995 [-c + d*i, a - b*i]])
996 blocks.append(_embed_complex_matrix(cplx_matrix))
997
998 # We should have real entries by now, so use the realest field
999 # we've got for the return value.
1000 return matrix.block(quaternions.base_ring(), n, blocks)
1001
1002
1003 def _unembed_quaternion_matrix(M):
1004 """
1005 The inverse of _embed_quaternion_matrix().
1006
1007 SETUP::
1008
1009 sage: from mjo.eja.eja_algebra import (_embed_quaternion_matrix,
1010 ....: _unembed_quaternion_matrix)
1011
1012 EXAMPLES::
1013
1014 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
1015 ....: [-2, 1, -4, 3],
1016 ....: [-3, 4, 1, -2],
1017 ....: [-4, -3, 2, 1]])
1018 sage: _unembed_quaternion_matrix(M)
1019 [1 + 2*i + 3*j + 4*k]
1020
1021 TESTS:
1022
1023 Unembedding is the inverse of embedding::
1024
1025 sage: set_random_seed()
1026 sage: Q = QuaternionAlgebra(QQ, -1, -1)
1027 sage: M = random_matrix(Q, 3)
1028 sage: _unembed_quaternion_matrix(_embed_quaternion_matrix(M)) == M
1029 True
1030
1031 """
1032 n = ZZ(M.nrows())
1033 if M.ncols() != n:
1034 raise ValueError("the matrix 'M' must be square")
1035 if not n.mod(4).is_zero():
1036 raise ValueError("the matrix 'M' must be a complex embedding")
1037
1038 Q = QuaternionAlgebra(QQ,-1,-1)
1039 i,j,k = Q.gens()
1040
1041 # Go top-left to bottom-right (reading order), converting every
1042 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
1043 # quaternion block.
1044 elements = []
1045 for l in xrange(n/4):
1046 for m in xrange(n/4):
1047 submat = _unembed_complex_matrix(M[4*l:4*l+4,4*m:4*m+4])
1048 if submat[0,0] != submat[1,1].conjugate():
1049 raise ValueError('bad on-diagonal submatrix')
1050 if submat[0,1] != -submat[1,0].conjugate():
1051 raise ValueError('bad off-diagonal submatrix')
1052 z = submat[0,0].real() + submat[0,0].imag()*i
1053 z += submat[0,1].real()*j + submat[0,1].imag()*k
1054 elements.append(z)
1055
1056 return matrix(Q, n/4, elements)
1057
1058
1059 # The usual inner product on R^n.
1060 def _usual_ip(x,y):
1061 return x.to_vector().inner_product(y.to_vector())
1062
1063 # The inner product used for the real symmetric simple EJA.
1064 # We keep it as a separate function because e.g. the complex
1065 # algebra uses the same inner product, except divided by 2.
1066 def _matrix_ip(X,Y):
1067 X_mat = X.natural_representation()
1068 Y_mat = Y.natural_representation()
1069 return (X_mat*Y_mat).trace()
1070
1071
1072 class RealSymmetricEJA(FiniteDimensionalEuclideanJordanAlgebra):
1073 """
1074 The rank-n simple EJA consisting of real symmetric n-by-n
1075 matrices, the usual symmetric Jordan product, and the trace inner
1076 product. It has dimension `(n^2 + n)/2` over the reals.
1077
1078 SETUP::
1079
1080 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1081
1082 EXAMPLES::
1083
1084 sage: J = RealSymmetricEJA(2)
1085 sage: e0, e1, e2 = J.gens()
1086 sage: e0*e0
1087 e0
1088 sage: e1*e1
1089 e0 + e2
1090 sage: e2*e2
1091 e2
1092
1093 TESTS:
1094
1095 The dimension of this algebra is `(n^2 + n) / 2`::
1096
1097 sage: set_random_seed()
1098 sage: n = ZZ.random_element(1,5)
1099 sage: J = RealSymmetricEJA(n)
1100 sage: J.dimension() == (n^2 + n)/2
1101 True
1102
1103 The Jordan multiplication is what we think it is::
1104
1105 sage: set_random_seed()
1106 sage: n = ZZ.random_element(1,5)
1107 sage: J = RealSymmetricEJA(n)
1108 sage: x = J.random_element()
1109 sage: y = J.random_element()
1110 sage: actual = (x*y).natural_representation()
1111 sage: X = x.natural_representation()
1112 sage: Y = y.natural_representation()
1113 sage: expected = (X*Y + Y*X)/2
1114 sage: actual == expected
1115 True
1116 sage: J(expected) == x*y
1117 True
1118
1119 """
1120 def __init__(self, n, field=QQ):
1121 S = _real_symmetric_basis(n, field=field)
1122 Qs = _multiplication_table_from_matrix_basis(S)
1123
1124 fdeja = super(RealSymmetricEJA, self)
1125 return fdeja.__init__(field,
1126 Qs,
1127 rank=n,
1128 natural_basis=S)
1129
1130 def inner_product(self, x, y):
1131 return _matrix_ip(x,y)
1132
1133
1134 class ComplexHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra):
1135 """
1136 The rank-n simple EJA consisting of complex Hermitian n-by-n
1137 matrices over the real numbers, the usual symmetric Jordan product,
1138 and the real-part-of-trace inner product. It has dimension `n^2` over
1139 the reals.
1140
1141 SETUP::
1142
1143 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1144
1145 TESTS:
1146
1147 The dimension of this algebra is `n^2`::
1148
1149 sage: set_random_seed()
1150 sage: n = ZZ.random_element(1,5)
1151 sage: J = ComplexHermitianEJA(n)
1152 sage: J.dimension() == n^2
1153 True
1154
1155 The Jordan multiplication is what we think it is::
1156
1157 sage: set_random_seed()
1158 sage: n = ZZ.random_element(1,5)
1159 sage: J = ComplexHermitianEJA(n)
1160 sage: x = J.random_element()
1161 sage: y = J.random_element()
1162 sage: actual = (x*y).natural_representation()
1163 sage: X = x.natural_representation()
1164 sage: Y = y.natural_representation()
1165 sage: expected = (X*Y + Y*X)/2
1166 sage: actual == expected
1167 True
1168 sage: J(expected) == x*y
1169 True
1170
1171 """
1172 def __init__(self, n, field=QQ):
1173 S = _complex_hermitian_basis(n)
1174 Qs = _multiplication_table_from_matrix_basis(S)
1175
1176 fdeja = super(ComplexHermitianEJA, self)
1177 return fdeja.__init__(field,
1178 Qs,
1179 rank=n,
1180 natural_basis=S)
1181
1182
1183 def inner_product(self, x, y):
1184 # Since a+bi on the diagonal is represented as
1185 #
1186 # a + bi = [ a b ]
1187 # [ -b a ],
1188 #
1189 # we'll double-count the "a" entries if we take the trace of
1190 # the embedding.
1191 return _matrix_ip(x,y)/2
1192
1193
1194 class QuaternionHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra):
1195 """
1196 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
1197 matrices, the usual symmetric Jordan product, and the
1198 real-part-of-trace inner product. It has dimension `2n^2 - n` over
1199 the reals.
1200
1201 SETUP::
1202
1203 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1204
1205 TESTS:
1206
1207 The dimension of this algebra is `n^2`::
1208
1209 sage: set_random_seed()
1210 sage: n = ZZ.random_element(1,5)
1211 sage: J = QuaternionHermitianEJA(n)
1212 sage: J.dimension() == 2*(n^2) - n
1213 True
1214
1215 The Jordan multiplication is what we think it is::
1216
1217 sage: set_random_seed()
1218 sage: n = ZZ.random_element(1,5)
1219 sage: J = QuaternionHermitianEJA(n)
1220 sage: x = J.random_element()
1221 sage: y = J.random_element()
1222 sage: actual = (x*y).natural_representation()
1223 sage: X = x.natural_representation()
1224 sage: Y = y.natural_representation()
1225 sage: expected = (X*Y + Y*X)/2
1226 sage: actual == expected
1227 True
1228 sage: J(expected) == x*y
1229 True
1230
1231 """
1232 def __init__(self, n, field=QQ):
1233 S = _quaternion_hermitian_basis(n)
1234 Qs = _multiplication_table_from_matrix_basis(S)
1235
1236 fdeja = super(QuaternionHermitianEJA, self)
1237 return fdeja.__init__(field,
1238 Qs,
1239 rank=n,
1240 natural_basis=S)
1241
1242 def inner_product(self, x, y):
1243 # Since a+bi+cj+dk on the diagonal is represented as
1244 #
1245 # a + bi +cj + dk = [ a b c d]
1246 # [ -b a -d c]
1247 # [ -c d a -b]
1248 # [ -d -c b a],
1249 #
1250 # we'll quadruple-count the "a" entries if we take the trace of
1251 # the embedding.
1252 return _matrix_ip(x,y)/4
1253
1254
1255 class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra):
1256 """
1257 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
1258 with the usual inner product and jordan product ``x*y =
1259 (<x_bar,y_bar>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
1260 the reals.
1261
1262 SETUP::
1263
1264 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1265
1266 EXAMPLES:
1267
1268 This multiplication table can be verified by hand::
1269
1270 sage: J = JordanSpinEJA(4)
1271 sage: e0,e1,e2,e3 = J.gens()
1272 sage: e0*e0
1273 e0
1274 sage: e0*e1
1275 e1
1276 sage: e0*e2
1277 e2
1278 sage: e0*e3
1279 e3
1280 sage: e1*e2
1281 0
1282 sage: e1*e3
1283 0
1284 sage: e2*e3
1285 0
1286
1287 """
1288 def __init__(self, n, field=QQ):
1289 V = VectorSpace(field, n)
1290 mult_table = [[V.zero() for i in range(n)] for j in range(n)]
1291 for i in range(n):
1292 for j in range(n):
1293 x = V.basis()[i]
1294 y = V.basis()[j]
1295 x0 = x[0]
1296 xbar = x[1:]
1297 y0 = y[0]
1298 ybar = y[1:]
1299 # z = x*y
1300 z0 = x.inner_product(y)
1301 zbar = y0*xbar + x0*ybar
1302 z = V([z0] + zbar.list())
1303 mult_table[i][j] = z
1304
1305 # The rank of the spin algebra is two, unless we're in a
1306 # one-dimensional ambient space (because the rank is bounded by
1307 # the ambient dimension).
1308 fdeja = super(JordanSpinEJA, self)
1309 return fdeja.__init__(field, mult_table, rank=min(n,2))
1310
1311 def inner_product(self, x, y):
1312 return _usual_ip(x,y)