2 Euclidean Jordan Algebras. These are formally-real Jordan Algebras;
3 specifically those where u^2 + v^2 = 0 implies that u = v = 0. They
4 are used in optimization, and have some additional nice methods beyond
5 what can be supported in a general Jordan Algebra.
8 from itertools
import repeat
10 from sage
.algebras
.quatalg
.quaternion_algebra
import QuaternionAlgebra
11 from sage
.categories
.magmatic_algebras
import MagmaticAlgebras
12 from sage
.combinat
.free_module
import CombinatorialFreeModule
13 from sage
.matrix
.constructor
import matrix
14 from sage
.matrix
.matrix_space
import MatrixSpace
15 from sage
.misc
.cachefunc
import cached_method
16 from sage
.misc
.lazy_import
import lazy_import
17 from sage
.misc
.prandom
import choice
18 from sage
.misc
.table
import table
19 from sage
.modules
.free_module
import FreeModule
, VectorSpace
20 from sage
.rings
.all
import (ZZ
, QQ
, AA
, QQbar
, RR
, RLF
, CLF
,
23 from mjo
.eja
.eja_element
import FiniteDimensionalEuclideanJordanAlgebraElement
24 lazy_import('mjo.eja.eja_subalgebra',
25 'FiniteDimensionalEuclideanJordanSubalgebra')
26 from mjo
.eja
.eja_utils
import _mat2vec
28 class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule
):
30 def _coerce_map_from_base_ring(self
):
32 Disable the map from the base ring into the algebra.
34 Performing a nonsense conversion like this automatically
35 is counterpedagogical. The fallback is to try the usual
36 element constructor, which should also fail.
40 sage: from mjo.eja.eja_algebra import random_eja
44 sage: set_random_seed()
45 sage: J = random_eja()
47 Traceback (most recent call last):
49 ValueError: not a naturally-represented algebra element
65 sage: from mjo.eja.eja_algebra import (
66 ....: FiniteDimensionalEuclideanJordanAlgebra,
72 By definition, Jordan multiplication commutes::
74 sage: set_random_seed()
75 sage: J = random_eja()
76 sage: x,y = J.random_elements(2)
82 The ``field`` we're given must be real with ``check_field=True``::
84 sage: JordanSpinEJA(2,QQbar)
85 Traceback (most recent call last):
87 ValueError: scalar field is not real
89 The multiplication table must be square with ``check_axioms=True``::
91 sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,((),()))
92 Traceback (most recent call last):
94 ValueError: multiplication table is not square
98 if not field
.is_subring(RR
):
99 # Note: this does return true for the real algebraic
100 # field, the rationals, and any quadratic field where
101 # we've specified a real embedding.
102 raise ValueError("scalar field is not real")
104 # The multiplication table had better be square
107 if not all( len(l
) == n
for l
in mult_table
):
108 raise ValueError("multiplication table is not square")
110 self
._natural
_basis
= natural_basis
113 category
= MagmaticAlgebras(field
).FiniteDimensional()
114 category
= category
.WithBasis().Unital()
116 fda
= super(FiniteDimensionalEuclideanJordanAlgebra
, self
)
121 self
.print_options(bracket
='')
123 # The multiplication table we're given is necessarily in terms
124 # of vectors, because we don't have an algebra yet for
125 # anything to be an element of. However, it's faster in the
126 # long run to have the multiplication table be in terms of
127 # algebra elements. We do this after calling the superclass
128 # constructor so that from_vector() knows what to do.
129 self
._multiplication
_table
= [
130 list(map(lambda x
: self
.from_vector(x
), ls
))
135 if not self
._is
_commutative
():
136 raise ValueError("algebra is not commutative")
137 if not self
._is
_jordanian
():
138 raise ValueError("Jordan identity does not hold")
139 if not self
._inner
_product
_is
_associative
():
140 raise ValueError("inner product is not associative")
142 def _element_constructor_(self
, elt
):
144 Construct an element of this algebra from its natural
147 This gets called only after the parent element _call_ method
148 fails to find a coercion for the argument.
152 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
154 ....: RealSymmetricEJA)
158 The identity in `S^n` is converted to the identity in the EJA::
160 sage: J = RealSymmetricEJA(3)
161 sage: I = matrix.identity(QQ,3)
162 sage: J(I) == J.one()
165 This skew-symmetric matrix can't be represented in the EJA::
167 sage: J = RealSymmetricEJA(3)
168 sage: A = matrix(QQ,3, lambda i,j: i-j)
170 Traceback (most recent call last):
172 ArithmeticError: vector is not in free module
176 Ensure that we can convert any element of the two non-matrix
177 simple algebras (whose natural representations are their usual
178 vector representations) back and forth faithfully::
180 sage: set_random_seed()
181 sage: J = HadamardEJA.random_instance()
182 sage: x = J.random_element()
183 sage: J(x.to_vector().column()) == x
185 sage: J = JordanSpinEJA.random_instance()
186 sage: x = J.random_element()
187 sage: J(x.to_vector().column()) == x
191 msg
= "not a naturally-represented algebra element"
193 # The superclass implementation of random_element()
194 # needs to be able to coerce "0" into the algebra.
196 elif elt
in self
.base_ring():
197 # Ensure that no base ring -> algebra coercion is performed
198 # by this method. There's some stupidity in sage that would
199 # otherwise propagate to this method; for example, sage thinks
200 # that the integer 3 belongs to the space of 2-by-2 matrices.
201 raise ValueError(msg
)
203 natural_basis
= self
.natural_basis()
204 basis_space
= natural_basis
[0].matrix_space()
205 if elt
not in basis_space
:
206 raise ValueError(msg
)
208 # Thanks for nothing! Matrix spaces aren't vector spaces in
209 # Sage, so we have to figure out its natural-basis coordinates
210 # ourselves. We use the basis space's ring instead of the
211 # element's ring because the basis space might be an algebraic
212 # closure whereas the base ring of the 3-by-3 identity matrix
213 # could be QQ instead of QQbar.
214 V
= VectorSpace(basis_space
.base_ring(), elt
.nrows()*elt
.ncols())
215 W
= V
.span_of_basis( _mat2vec(s
) for s
in natural_basis
)
216 coords
= W
.coordinate_vector(_mat2vec(elt
))
217 return self
.from_vector(coords
)
220 def _max_random_instance_size():
222 Return an integer "size" that is an upper bound on the size of
223 this algebra when it is used in a random test
224 case. Unfortunately, the term "size" is quite vague -- when
225 dealing with `R^n` under either the Hadamard or Jordan spin
226 product, the "size" refers to the dimension `n`. When dealing
227 with a matrix algebra (real symmetric or complex/quaternion
228 Hermitian), it refers to the size of the matrix, which is
229 far less than the dimension of the underlying vector space.
231 We default to five in this class, which is safe in `R^n`. The
232 matrix algebra subclasses (or any class where the "size" is
233 interpreted to be far less than the dimension) should override
234 with a smaller number.
236 raise NotImplementedError
240 Return a string representation of ``self``.
244 sage: from mjo.eja.eja_algebra import JordanSpinEJA
248 Ensure that it says what we think it says::
250 sage: JordanSpinEJA(2, field=AA)
251 Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
252 sage: JordanSpinEJA(3, field=RDF)
253 Euclidean Jordan algebra of dimension 3 over Real Double Field
256 fmt
= "Euclidean Jordan algebra of dimension {} over {}"
257 return fmt
.format(self
.dimension(), self
.base_ring())
259 def product_on_basis(self
, i
, j
):
260 return self
._multiplication
_table
[i
][j
]
262 def _is_commutative(self
):
264 Whether or not this algebra's multiplication table is commutative.
266 This method should of course always return ``True``, unless
267 this algebra was constructed with ``check_axioms=False`` and
268 passed an invalid multiplication table.
270 return all( self
.product_on_basis(i
,j
) == self
.product_on_basis(i
,j
)
271 for i
in range(self
.dimension())
272 for j
in range(self
.dimension()) )
274 def _is_jordanian(self
):
276 Whether or not this algebra's multiplication table respects the
277 Jordan identity `(x^{2})(xy) = x(x^{2}y)`.
279 We only check one arrangement of `x` and `y`, so for a
280 ``True`` result to be truly true, you should also check
281 :meth:`_is_commutative`. This method should of course always
282 return ``True``, unless this algebra was constructed with
283 ``check_axioms=False`` and passed an invalid multiplication table.
285 return all( (self
.monomial(i
)**2)*(self
.monomial(i
)*self
.monomial(j
))
287 (self
.monomial(i
))*((self
.monomial(i
)**2)*self
.monomial(j
))
288 for i
in range(self
.dimension())
289 for j
in range(self
.dimension()) )
291 def _inner_product_is_associative(self
):
293 Return whether or not this algebra's inner product `B` is
294 associative; that is, whether or not `B(xy,z) = B(x,yz)`.
296 This method should of course always return ``True``, unless
297 this algebra was constructed with ``check_axioms=False`` and
298 passed an invalid multiplication table.
301 # Used to check whether or not something is zero in an inexact
302 # ring. This number is sufficient to allow the construction of
303 # QuaternionHermitianEJA(2, RDF) with check_axioms=True.
306 for i
in range(self
.dimension()):
307 for j
in range(self
.dimension()):
308 for k
in range(self
.dimension()):
312 diff
= (x
*y
).inner_product(z
) - x
.inner_product(y
*z
)
314 if self
.base_ring().is_exact():
318 if diff
.abs() > epsilon
:
324 def characteristic_polynomial_of(self
):
326 Return the algebra's "characteristic polynomial of" function,
327 which is itself a multivariate polynomial that, when evaluated
328 at the coordinates of some algebra element, returns that
329 element's characteristic polynomial.
331 The resulting polynomial has `n+1` variables, where `n` is the
332 dimension of this algebra. The first `n` variables correspond to
333 the coordinates of an algebra element: when evaluated at the
334 coordinates of an algebra element with respect to a certain
335 basis, the result is a univariate polynomial (in the one
336 remaining variable ``t``), namely the characteristic polynomial
341 sage: from mjo.eja.eja_algebra import JordanSpinEJA, TrivialEJA
345 The characteristic polynomial in the spin algebra is given in
346 Alizadeh, Example 11.11::
348 sage: J = JordanSpinEJA(3)
349 sage: p = J.characteristic_polynomial_of(); p
350 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
351 sage: xvec = J.one().to_vector()
355 By definition, the characteristic polynomial is a monic
356 degree-zero polynomial in a rank-zero algebra. Note that
357 Cayley-Hamilton is indeed satisfied since the polynomial
358 ``1`` evaluates to the identity element of the algebra on
361 sage: J = TrivialEJA()
362 sage: J.characteristic_polynomial_of()
369 # The list of coefficient polynomials a_0, a_1, a_2, ..., a_(r-1).
370 a
= self
._charpoly
_coefficients
()
372 # We go to a bit of trouble here to reorder the
373 # indeterminates, so that it's easier to evaluate the
374 # characteristic polynomial at x's coordinates and get back
375 # something in terms of t, which is what we want.
376 S
= PolynomialRing(self
.base_ring(),'t')
380 S
= PolynomialRing(S
, R
.variable_names())
383 return (t
**r
+ sum( a
[k
]*(t
**k
) for k
in range(r
) ))
386 def inner_product(self
, x
, y
):
388 The inner product associated with this Euclidean Jordan algebra.
390 Defaults to the trace inner product, but can be overridden by
391 subclasses if they are sure that the necessary properties are
396 sage: from mjo.eja.eja_algebra import random_eja
400 Our inner product is "associative," which means the following for
401 a symmetric bilinear form::
403 sage: set_random_seed()
404 sage: J = random_eja()
405 sage: x,y,z = J.random_elements(3)
406 sage: (x*y).inner_product(z) == y.inner_product(x*z)
410 X
= x
.natural_representation()
411 Y
= y
.natural_representation()
412 return self
.natural_inner_product(X
,Y
)
415 def is_trivial(self
):
417 Return whether or not this algebra is trivial.
419 A trivial algebra contains only the zero element.
423 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
428 sage: J = ComplexHermitianEJA(3)
434 sage: J = TrivialEJA()
439 return self
.dimension() == 0
442 def multiplication_table(self
):
444 Return a visual representation of this algebra's multiplication
445 table (on basis elements).
449 sage: from mjo.eja.eja_algebra import JordanSpinEJA
453 sage: J = JordanSpinEJA(4)
454 sage: J.multiplication_table()
455 +----++----+----+----+----+
456 | * || e0 | e1 | e2 | e3 |
457 +====++====+====+====+====+
458 | e0 || e0 | e1 | e2 | e3 |
459 +----++----+----+----+----+
460 | e1 || e1 | e0 | 0 | 0 |
461 +----++----+----+----+----+
462 | e2 || e2 | 0 | e0 | 0 |
463 +----++----+----+----+----+
464 | e3 || e3 | 0 | 0 | e0 |
465 +----++----+----+----+----+
468 M
= list(self
._multiplication
_table
) # copy
469 for i
in range(len(M
)):
470 # M had better be "square"
471 M
[i
] = [self
.monomial(i
)] + M
[i
]
472 M
= [["*"] + list(self
.gens())] + M
473 return table(M
, header_row
=True, header_column
=True, frame
=True)
476 def natural_basis(self
):
478 Return a more-natural representation of this algebra's basis.
480 Every finite-dimensional Euclidean Jordan Algebra is a direct
481 sum of five simple algebras, four of which comprise Hermitian
482 matrices. This method returns the original "natural" basis
483 for our underlying vector space. (Typically, the natural basis
484 is used to construct the multiplication table in the first place.)
486 Note that this will always return a matrix. The standard basis
487 in `R^n` will be returned as `n`-by-`1` column matrices.
491 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
492 ....: RealSymmetricEJA)
496 sage: J = RealSymmetricEJA(2)
498 Finite family {0: e0, 1: e1, 2: e2}
499 sage: J.natural_basis()
501 [1 0] [ 0 0.7071067811865475?] [0 0]
502 [0 0], [0.7071067811865475? 0], [0 1]
507 sage: J = JordanSpinEJA(2)
509 Finite family {0: e0, 1: e1}
510 sage: J.natural_basis()
517 if self
._natural
_basis
is None:
518 M
= self
.natural_basis_space()
519 return tuple( M(b
.to_vector()) for b
in self
.basis() )
521 return self
._natural
_basis
524 def natural_basis_space(self
):
526 Return the matrix space in which this algebra's natural basis
529 Generally this will be an `n`-by-`1` column-vector space,
530 except when the algebra is trivial. There it's `n`-by-`n`
531 (where `n` is zero), to ensure that two elements of the
532 natural basis space (empty matrices) can be multiplied.
534 if self
.is_trivial():
535 return MatrixSpace(self
.base_ring(), 0)
536 elif self
._natural
_basis
is None or len(self
._natural
_basis
) == 0:
537 return MatrixSpace(self
.base_ring(), self
.dimension(), 1)
539 return self
._natural
_basis
[0].matrix_space()
543 def natural_inner_product(X
,Y
):
545 Compute the inner product of two naturally-represented elements.
547 For example in the real symmetric matrix EJA, this will compute
548 the trace inner-product of two n-by-n symmetric matrices. The
549 default should work for the real cartesian product EJA, the
550 Jordan spin EJA, and the real symmetric matrices. The others
551 will have to be overridden.
553 return (X
.conjugate_transpose()*Y
).trace()
559 Return the unit element of this algebra.
563 sage: from mjo.eja.eja_algebra import (HadamardEJA,
568 sage: J = HadamardEJA(5)
570 e0 + e1 + e2 + e3 + e4
574 The identity element acts like the identity::
576 sage: set_random_seed()
577 sage: J = random_eja()
578 sage: x = J.random_element()
579 sage: J.one()*x == x and x*J.one() == x
582 The matrix of the unit element's operator is the identity::
584 sage: set_random_seed()
585 sage: J = random_eja()
586 sage: actual = J.one().operator().matrix()
587 sage: expected = matrix.identity(J.base_ring(), J.dimension())
588 sage: actual == expected
591 Ensure that the cached unit element (often precomputed by
592 hand) agrees with the computed one::
594 sage: set_random_seed()
595 sage: J = random_eja()
596 sage: cached = J.one()
597 sage: J.one.clear_cache()
598 sage: J.one() == cached
602 # We can brute-force compute the matrices of the operators
603 # that correspond to the basis elements of this algebra.
604 # If some linear combination of those basis elements is the
605 # algebra identity, then the same linear combination of
606 # their matrices has to be the identity matrix.
608 # Of course, matrices aren't vectors in sage, so we have to
609 # appeal to the "long vectors" isometry.
610 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
612 # Now we use basic linear algebra to find the coefficients,
613 # of the matrices-as-vectors-linear-combination, which should
614 # work for the original algebra basis too.
615 A
= matrix(self
.base_ring(), oper_vecs
)
617 # We used the isometry on the left-hand side already, but we
618 # still need to do it for the right-hand side. Recall that we
619 # wanted something that summed to the identity matrix.
620 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
622 # Now if there's an identity element in the algebra, this
623 # should work. We solve on the left to avoid having to
624 # transpose the matrix "A".
625 return self
.from_vector(A
.solve_left(b
))
628 def peirce_decomposition(self
, c
):
630 The Peirce decomposition of this algebra relative to the
633 In the future, this can be extended to a complete system of
634 orthogonal idempotents.
638 - ``c`` -- an idempotent of this algebra.
642 A triple (J0, J5, J1) containing two subalgebras and one subspace
645 - ``J0`` -- the algebra on the eigenspace of ``c.operator()``
646 corresponding to the eigenvalue zero.
648 - ``J5`` -- the eigenspace (NOT a subalgebra) of ``c.operator()``
649 corresponding to the eigenvalue one-half.
651 - ``J1`` -- the algebra on the eigenspace of ``c.operator()``
652 corresponding to the eigenvalue one.
654 These are the only possible eigenspaces for that operator, and this
655 algebra is a direct sum of them. The spaces ``J0`` and ``J1`` are
656 orthogonal, and are subalgebras of this algebra with the appropriate
661 sage: from mjo.eja.eja_algebra import random_eja, RealSymmetricEJA
665 The canonical example comes from the symmetric matrices, which
666 decompose into diagonal and off-diagonal parts::
668 sage: J = RealSymmetricEJA(3)
669 sage: C = matrix(QQ, [ [1,0,0],
673 sage: J0,J5,J1 = J.peirce_decomposition(c)
675 Euclidean Jordan algebra of dimension 1...
677 Vector space of degree 6 and dimension 2...
679 Euclidean Jordan algebra of dimension 3...
680 sage: J0.one().natural_representation()
684 sage: orig_df = AA.options.display_format
685 sage: AA.options.display_format = 'radical'
686 sage: J.from_vector(J5.basis()[0]).natural_representation()
690 sage: J.from_vector(J5.basis()[1]).natural_representation()
694 sage: AA.options.display_format = orig_df
695 sage: J1.one().natural_representation()
702 Every algebra decomposes trivially with respect to its identity
705 sage: set_random_seed()
706 sage: J = random_eja()
707 sage: J0,J5,J1 = J.peirce_decomposition(J.one())
708 sage: J0.dimension() == 0 and J5.dimension() == 0
710 sage: J1.superalgebra() == J and J1.dimension() == J.dimension()
713 The decomposition is into eigenspaces, and its components are
714 therefore necessarily orthogonal. Moreover, the identity
715 elements in the two subalgebras are the projections onto their
716 respective subspaces of the superalgebra's identity element::
718 sage: set_random_seed()
719 sage: J = random_eja()
720 sage: x = J.random_element()
721 sage: if not J.is_trivial():
722 ....: while x.is_nilpotent():
723 ....: x = J.random_element()
724 sage: c = x.subalgebra_idempotent()
725 sage: J0,J5,J1 = J.peirce_decomposition(c)
727 sage: for (w,y,z) in zip(J0.basis(), J5.basis(), J1.basis()):
728 ....: w = w.superalgebra_element()
729 ....: y = J.from_vector(y)
730 ....: z = z.superalgebra_element()
731 ....: ipsum += w.inner_product(y).abs()
732 ....: ipsum += w.inner_product(z).abs()
733 ....: ipsum += y.inner_product(z).abs()
736 sage: J1(c) == J1.one()
738 sage: J0(J.one() - c) == J0.one()
742 if not c
.is_idempotent():
743 raise ValueError("element is not idempotent: %s" % c
)
745 # Default these to what they should be if they turn out to be
746 # trivial, because eigenspaces_left() won't return eigenvalues
747 # corresponding to trivial spaces (e.g. it returns only the
748 # eigenspace corresponding to lambda=1 if you take the
749 # decomposition relative to the identity element).
750 trivial
= FiniteDimensionalEuclideanJordanSubalgebra(self
, ())
751 J0
= trivial
# eigenvalue zero
752 J5
= VectorSpace(self
.base_ring(), 0) # eigenvalue one-half
753 J1
= trivial
# eigenvalue one
755 for (eigval
, eigspace
) in c
.operator().matrix().right_eigenspaces():
756 if eigval
== ~
(self
.base_ring()(2)):
759 gens
= tuple( self
.from_vector(b
) for b
in eigspace
.basis() )
760 subalg
= FiniteDimensionalEuclideanJordanSubalgebra(self
,
768 raise ValueError("unexpected eigenvalue: %s" % eigval
)
773 def random_element(self
, thorough
=False):
775 Return a random element of this algebra.
777 Our algebra superclass method only returns a linear
778 combination of at most two basis elements. We instead
779 want the vector space "random element" method that
780 returns a more diverse selection.
784 - ``thorough`` -- (boolean; default False) whether or not we
785 should generate irrational coefficients for the random
786 element when our base ring is irrational; this slows the
787 algebra operations to a crawl, but any truly random method
791 # For a general base ring... maybe we can trust this to do the
792 # right thing? Unlikely, but.
793 V
= self
.vector_space()
794 v
= V
.random_element()
796 if self
.base_ring() is AA
:
797 # The "random element" method of the algebraic reals is
798 # stupid at the moment, and only returns integers between
799 # -2 and 2, inclusive:
801 # https://trac.sagemath.org/ticket/30875
803 # Instead, we implement our own "random vector" method,
804 # and then coerce that into the algebra. We use the vector
805 # space degree here instead of the dimension because a
806 # subalgebra could (for example) be spanned by only two
807 # vectors, each with five coordinates. We need to
808 # generate all five coordinates.
810 v
*= QQbar
.random_element().real()
812 v
*= QQ
.random_element()
814 return self
.from_vector(V
.coordinate_vector(v
))
816 def random_elements(self
, count
, thorough
=False):
818 Return ``count`` random elements as a tuple.
822 - ``thorough`` -- (boolean; default False) whether or not we
823 should generate irrational coefficients for the random
824 elements when our base ring is irrational; this slows the
825 algebra operations to a crawl, but any truly random method
830 sage: from mjo.eja.eja_algebra import JordanSpinEJA
834 sage: J = JordanSpinEJA(3)
835 sage: x,y,z = J.random_elements(3)
836 sage: all( [ x in J, y in J, z in J ])
838 sage: len( J.random_elements(10) ) == 10
842 return tuple( self
.random_element(thorough
)
843 for idx
in range(count
) )
846 def random_instance(cls
, field
=AA
, **kwargs
):
848 Return a random instance of this type of algebra.
850 Beware, this will crash for "most instances" because the
851 constructor below looks wrong.
853 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
854 return cls(n
, field
, **kwargs
)
857 def _charpoly_coefficients(self
):
859 The `r` polynomial coefficients of the "characteristic polynomial
863 var_names
= [ "X" + str(z
) for z
in range(1,n
+1) ]
864 R
= PolynomialRing(self
.base_ring(), var_names
)
866 F
= R
.fraction_field()
869 # From a result in my book, these are the entries of the
870 # basis representation of L_x.
871 return sum( vars[k
]*self
.monomial(k
).operator().matrix()[i
,j
]
874 L_x
= matrix(F
, n
, n
, L_x_i_j
)
877 if self
.rank
.is_in_cache():
879 # There's no need to pad the system with redundant
880 # columns if we *know* they'll be redundant.
883 # Compute an extra power in case the rank is equal to
884 # the dimension (otherwise, we would stop at x^(r-1)).
885 x_powers
= [ (L_x
**k
)*self
.one().to_vector()
886 for k
in range(n
+1) ]
887 A
= matrix
.column(F
, x_powers
[:n
])
888 AE
= A
.extended_echelon_form()
895 # The theory says that only the first "r" coefficients are
896 # nonzero, and they actually live in the original polynomial
897 # ring and not the fraction field. We negate them because
898 # in the actual characteristic polynomial, they get moved
899 # to the other side where x^r lives.
900 return -A_rref
.solve_right(E
*b
).change_ring(R
)[:r
]
905 Return the rank of this EJA.
907 This is a cached method because we know the rank a priori for
908 all of the algebras we can construct. Thus we can avoid the
909 expensive ``_charpoly_coefficients()`` call unless we truly
910 need to compute the whole characteristic polynomial.
914 sage: from mjo.eja.eja_algebra import (HadamardEJA,
916 ....: RealSymmetricEJA,
917 ....: ComplexHermitianEJA,
918 ....: QuaternionHermitianEJA,
923 The rank of the Jordan spin algebra is always two::
925 sage: JordanSpinEJA(2).rank()
927 sage: JordanSpinEJA(3).rank()
929 sage: JordanSpinEJA(4).rank()
932 The rank of the `n`-by-`n` Hermitian real, complex, or
933 quaternion matrices is `n`::
935 sage: RealSymmetricEJA(4).rank()
937 sage: ComplexHermitianEJA(3).rank()
939 sage: QuaternionHermitianEJA(2).rank()
944 Ensure that every EJA that we know how to construct has a
945 positive integer rank, unless the algebra is trivial in
946 which case its rank will be zero::
948 sage: set_random_seed()
949 sage: J = random_eja()
953 sage: r > 0 or (r == 0 and J.is_trivial())
956 Ensure that computing the rank actually works, since the ranks
957 of all simple algebras are known and will be cached by default::
959 sage: J = HadamardEJA(4)
960 sage: J.rank.clear_cache()
966 sage: J = JordanSpinEJA(4)
967 sage: J.rank.clear_cache()
973 sage: J = RealSymmetricEJA(3)
974 sage: J.rank.clear_cache()
980 sage: J = ComplexHermitianEJA(2)
981 sage: J.rank.clear_cache()
987 sage: J = QuaternionHermitianEJA(2)
988 sage: J.rank.clear_cache()
992 return len(self
._charpoly
_coefficients
())
995 def vector_space(self
):
997 Return the vector space that underlies this algebra.
1001 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1005 sage: J = RealSymmetricEJA(2)
1006 sage: J.vector_space()
1007 Vector space of dimension 3 over...
1010 return self
.zero().to_vector().parent().ambient_vector_space()
1013 Element
= FiniteDimensionalEuclideanJordanAlgebraElement
1017 def random_eja(field
=AA
):
1019 Return a "random" finite-dimensional Euclidean Jordan Algebra.
1023 sage: from mjo.eja.eja_algebra import random_eja
1028 Euclidean Jordan algebra of dimension...
1031 classname
= choice([TrivialEJA
,
1035 ComplexHermitianEJA
,
1036 QuaternionHermitianEJA
])
1037 return classname
.random_instance(field
=field
)
1042 class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra
):
1044 Algebras whose basis consists of vectors with rational
1045 entries. Equivalently, algebras whose multiplication tables
1046 contain only rational coefficients.
1048 When an EJA has a basis that can be made rational, we can speed up
1049 the computation of its characteristic polynomial by doing it over
1050 ``QQ``. All of the named EJA constructors that we provide fall
1054 def _charpoly_coefficients(self
):
1056 Override the parent method with something that tries to compute
1057 over a faster (non-extension) field.
1061 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1065 The base ring of the resulting polynomial coefficients is what
1066 it should be, and not the rationals (unless the algebra was
1067 already over the rationals)::
1069 sage: J = JordanSpinEJA(3)
1070 sage: J._charpoly_coefficients()
1071 (X1^2 - X2^2 - X3^2, -2*X1)
1072 sage: a0 = J._charpoly_coefficients()[0]
1074 Algebraic Real Field
1075 sage: a0.base_ring()
1076 Algebraic Real Field
1079 if self
.base_ring() is QQ
:
1080 # There's no need to construct *another* algebra over the
1081 # rationals if this one is already over the rationals.
1082 superclass
= super(RationalBasisEuclideanJordanAlgebra
, self
)
1083 return superclass
._charpoly
_coefficients
()
1086 map(lambda x
: x
.to_vector(), ls
)
1087 for ls
in self
._multiplication
_table
1090 # Do the computation over the rationals. The answer will be
1091 # the same, because our basis coordinates are (essentially)
1093 J
= FiniteDimensionalEuclideanJordanAlgebra(QQ
,
1097 a
= J
._charpoly
_coefficients
()
1098 return tuple(map(lambda x
: x
.change_ring(self
.base_ring()), a
))
1101 class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra
):
1103 def _max_random_instance_size():
1104 # Play it safe, since this will be squared and the underlying
1105 # field can have dimension 4 (quaternions) too.
1108 def __init__(self
, field
, basis
, normalize_basis
=True, **kwargs
):
1110 Compared to the superclass constructor, we take a basis instead of
1111 a multiplication table because the latter can be computed in terms
1112 of the former when the product is known (like it is here).
1114 # Used in this class's fast _charpoly_coefficients() override.
1115 self
._basis
_normalizers
= None
1117 # We're going to loop through this a few times, so now's a good
1118 # time to ensure that it isn't a generator expression.
1119 basis
= tuple(basis
)
1121 algebra_dim
= len(basis
)
1122 if algebra_dim
> 1 and normalize_basis
:
1123 # We'll need sqrt(2) to normalize the basis, and this
1124 # winds up in the multiplication table, so the whole
1125 # algebra needs to be over the field extension.
1126 R
= PolynomialRing(field
, 'z')
1129 if p
.is_irreducible():
1130 field
= field
.extension(p
, 'sqrt2', embedding
=RLF(2).sqrt())
1131 basis
= tuple( s
.change_ring(field
) for s
in basis
)
1132 self
._basis
_normalizers
= tuple(
1133 ~
(self
.natural_inner_product(s
,s
).sqrt()) for s
in basis
)
1134 basis
= tuple(s
*c
for (s
,c
) in zip(basis
,self
._basis
_normalizers
))
1136 Qs
= self
.multiplication_table_from_matrix_basis(basis
)
1138 super(MatrixEuclideanJordanAlgebra
, self
).__init
__(field
,
1140 natural_basis
=basis
,
1143 if algebra_dim
== 0:
1144 self
.one
.set_cache(self
.zero())
1146 n
= basis
[0].nrows()
1147 # The identity wrt (A,B) -> (AB + BA)/2 is independent of the
1148 # details of this algebra.
1149 self
.one
.set_cache(self(matrix
.identity(field
,n
)))
1153 def _charpoly_coefficients(self
):
1155 Override the parent method with something that tries to compute
1156 over a faster (non-extension) field.
1158 if self
._basis
_normalizers
is None or self
.base_ring() is QQ
:
1159 # We didn't normalize, or the basis we started with had
1160 # entries in a nice field already. Just compute the thing.
1161 return super(MatrixEuclideanJordanAlgebra
, self
)._charpoly
_coefficients
()
1163 basis
= ( (b
/n
) for (b
,n
) in zip(self
.natural_basis(),
1164 self
._basis
_normalizers
) )
1166 # Do this over the rationals and convert back at the end.
1167 # Only works because we know the entries of the basis are
1168 # integers. The argument ``check_axioms=False`` is required
1169 # because the trace inner-product method for this
1170 # class is a stub and can't actually be checked.
1171 J
= MatrixEuclideanJordanAlgebra(QQ
,
1173 normalize_basis
=False,
1176 a
= J
._charpoly
_coefficients
()
1178 # Unfortunately, changing the basis does change the
1179 # coefficients of the characteristic polynomial, but since
1180 # these are really the coefficients of the "characteristic
1181 # polynomial of" function, everything is still nice and
1182 # unevaluated. It's therefore "obvious" how scaling the
1183 # basis affects the coordinate variables X1, X2, et
1184 # cetera. Scaling the first basis vector up by "n" adds a
1185 # factor of 1/n into every "X1" term, for example. So here
1186 # we simply undo the basis_normalizer scaling that we
1187 # performed earlier.
1189 # The a[0] access here is safe because trivial algebras
1190 # won't have any basis normalizers and therefore won't
1191 # make it to this "else" branch.
1192 XS
= a
[0].parent().gens()
1193 subs_dict
= { XS
[i
]: self
._basis
_normalizers
[i
]*XS
[i
]
1194 for i
in range(len(XS
)) }
1195 return tuple( a_i
.subs(subs_dict
) for a_i
in a
)
1199 def multiplication_table_from_matrix_basis(basis
):
1201 At least three of the five simple Euclidean Jordan algebras have the
1202 symmetric multiplication (A,B) |-> (AB + BA)/2, where the
1203 multiplication on the right is matrix multiplication. Given a basis
1204 for the underlying matrix space, this function returns a
1205 multiplication table (obtained by looping through the basis
1206 elements) for an algebra of those matrices.
1208 # In S^2, for example, we nominally have four coordinates even
1209 # though the space is of dimension three only. The vector space V
1210 # is supposed to hold the entire long vector, and the subspace W
1211 # of V will be spanned by the vectors that arise from symmetric
1212 # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
1216 field
= basis
[0].base_ring()
1217 dimension
= basis
[0].nrows()
1219 V
= VectorSpace(field
, dimension
**2)
1220 W
= V
.span_of_basis( _mat2vec(s
) for s
in basis
)
1222 mult_table
= [[W
.zero() for j
in range(n
)] for i
in range(n
)]
1225 mat_entry
= (basis
[i
]*basis
[j
] + basis
[j
]*basis
[i
])/2
1226 mult_table
[i
][j
] = W
.coordinate_vector(_mat2vec(mat_entry
))
1234 Embed the matrix ``M`` into a space of real matrices.
1236 The matrix ``M`` can have entries in any field at the moment:
1237 the real numbers, complex numbers, or quaternions. And although
1238 they are not a field, we can probably support octonions at some
1239 point, too. This function returns a real matrix that "acts like"
1240 the original with respect to matrix multiplication; i.e.
1242 real_embed(M*N) = real_embed(M)*real_embed(N)
1245 raise NotImplementedError
1249 def real_unembed(M
):
1251 The inverse of :meth:`real_embed`.
1253 raise NotImplementedError
1257 def natural_inner_product(cls
,X
,Y
):
1258 Xu
= cls
.real_unembed(X
)
1259 Yu
= cls
.real_unembed(Y
)
1260 tr
= (Xu
*Yu
).trace()
1263 # Works in QQ, AA, RDF, et cetera.
1265 except AttributeError:
1266 # A quaternion doesn't have a real() method, but does
1267 # have coefficient_tuple() method that returns the
1268 # coefficients of 1, i, j, and k -- in that order.
1269 return tr
.coefficient_tuple()[0]
1272 class RealMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1276 The identity function, for embedding real matrices into real
1282 def real_unembed(M
):
1284 The identity function, for unembedding real matrices from real
1290 class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra
):
1292 The rank-n simple EJA consisting of real symmetric n-by-n
1293 matrices, the usual symmetric Jordan product, and the trace inner
1294 product. It has dimension `(n^2 + n)/2` over the reals.
1298 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1302 sage: J = RealSymmetricEJA(2)
1303 sage: e0, e1, e2 = J.gens()
1311 In theory, our "field" can be any subfield of the reals::
1313 sage: RealSymmetricEJA(2, RDF)
1314 Euclidean Jordan algebra of dimension 3 over Real Double Field
1315 sage: RealSymmetricEJA(2, RR)
1316 Euclidean Jordan algebra of dimension 3 over Real Field with
1317 53 bits of precision
1321 The dimension of this algebra is `(n^2 + n) / 2`::
1323 sage: set_random_seed()
1324 sage: n_max = RealSymmetricEJA._max_random_instance_size()
1325 sage: n = ZZ.random_element(1, n_max)
1326 sage: J = RealSymmetricEJA(n)
1327 sage: J.dimension() == (n^2 + n)/2
1330 The Jordan multiplication is what we think it is::
1332 sage: set_random_seed()
1333 sage: J = RealSymmetricEJA.random_instance()
1334 sage: x,y = J.random_elements(2)
1335 sage: actual = (x*y).natural_representation()
1336 sage: X = x.natural_representation()
1337 sage: Y = y.natural_representation()
1338 sage: expected = (X*Y + Y*X)/2
1339 sage: actual == expected
1341 sage: J(expected) == x*y
1344 We can change the generator prefix::
1346 sage: RealSymmetricEJA(3, prefix='q').gens()
1347 (q0, q1, q2, q3, q4, q5)
1349 Our natural basis is normalized with respect to the natural inner
1350 product unless we specify otherwise::
1352 sage: set_random_seed()
1353 sage: J = RealSymmetricEJA.random_instance()
1354 sage: all( b.norm() == 1 for b in J.gens() )
1357 Since our natural basis is normalized with respect to the natural
1358 inner product, and since we know that this algebra is an EJA, any
1359 left-multiplication operator's matrix will be symmetric because
1360 natural->EJA basis representation is an isometry and within the EJA
1361 the operator is self-adjoint by the Jordan axiom::
1363 sage: set_random_seed()
1364 sage: x = RealSymmetricEJA.random_instance().random_element()
1365 sage: x.operator().matrix().is_symmetric()
1368 We can construct the (trivial) algebra of rank zero::
1370 sage: RealSymmetricEJA(0)
1371 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1375 def _denormalized_basis(cls
, n
, field
):
1377 Return a basis for the space of real symmetric n-by-n matrices.
1381 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1385 sage: set_random_seed()
1386 sage: n = ZZ.random_element(1,5)
1387 sage: B = RealSymmetricEJA._denormalized_basis(n,QQ)
1388 sage: all( M.is_symmetric() for M in B)
1392 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1396 for j
in range(i
+1):
1397 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
1401 Sij
= Eij
+ Eij
.transpose()
1407 def _max_random_instance_size():
1408 return 4 # Dimension 10
1411 def __init__(self
, n
, field
=AA
, **kwargs
):
1412 basis
= self
._denormalized
_basis
(n
, field
)
1413 super(RealSymmetricEJA
, self
).__init
__(field
,
1417 self
.rank
.set_cache(n
)
1420 class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1424 Embed the n-by-n complex matrix ``M`` into the space of real
1425 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1426 bi` to the block matrix ``[[a,b],[-b,a]]``.
1430 sage: from mjo.eja.eja_algebra import \
1431 ....: ComplexMatrixEuclideanJordanAlgebra
1435 sage: F = QuadraticField(-1, 'I')
1436 sage: x1 = F(4 - 2*i)
1437 sage: x2 = F(1 + 2*i)
1440 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1441 sage: ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1450 Embedding is a homomorphism (isomorphism, in fact)::
1452 sage: set_random_seed()
1453 sage: n_max = ComplexMatrixEuclideanJordanAlgebra._max_random_instance_size()
1454 sage: n = ZZ.random_element(n_max)
1455 sage: F = QuadraticField(-1, 'I')
1456 sage: X = random_matrix(F, n)
1457 sage: Y = random_matrix(F, n)
1458 sage: Xe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X)
1459 sage: Ye = ComplexMatrixEuclideanJordanAlgebra.real_embed(Y)
1460 sage: XYe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1467 raise ValueError("the matrix 'M' must be square")
1469 # We don't need any adjoined elements...
1470 field
= M
.base_ring().base_ring()
1474 a
= z
.list()[0] # real part, I guess
1475 b
= z
.list()[1] # imag part, I guess
1476 blocks
.append(matrix(field
, 2, [[a
,b
],[-b
,a
]]))
1478 return matrix
.block(field
, n
, blocks
)
1482 def real_unembed(M
):
1484 The inverse of _embed_complex_matrix().
1488 sage: from mjo.eja.eja_algebra import \
1489 ....: ComplexMatrixEuclideanJordanAlgebra
1493 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1494 ....: [-2, 1, -4, 3],
1495 ....: [ 9, 10, 11, 12],
1496 ....: [-10, 9, -12, 11] ])
1497 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(A)
1499 [ 10*I + 9 12*I + 11]
1503 Unembedding is the inverse of embedding::
1505 sage: set_random_seed()
1506 sage: F = QuadraticField(-1, 'I')
1507 sage: M = random_matrix(F, 3)
1508 sage: Me = ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1509 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1515 raise ValueError("the matrix 'M' must be square")
1516 if not n
.mod(2).is_zero():
1517 raise ValueError("the matrix 'M' must be a complex embedding")
1519 # If "M" was normalized, its base ring might have roots
1520 # adjoined and they can stick around after unembedding.
1521 field
= M
.base_ring()
1522 R
= PolynomialRing(field
, 'z')
1525 # Sage doesn't know how to embed AA into QQbar, i.e. how
1526 # to adjoin sqrt(-1) to AA.
1529 F
= field
.extension(z
**2 + 1, 'I', embedding
=CLF(-1).sqrt())
1532 # Go top-left to bottom-right (reading order), converting every
1533 # 2-by-2 block we see to a single complex element.
1535 for k
in range(n
/2):
1536 for j
in range(n
/2):
1537 submat
= M
[2*k
:2*k
+2,2*j
:2*j
+2]
1538 if submat
[0,0] != submat
[1,1]:
1539 raise ValueError('bad on-diagonal submatrix')
1540 if submat
[0,1] != -submat
[1,0]:
1541 raise ValueError('bad off-diagonal submatrix')
1542 z
= submat
[0,0] + submat
[0,1]*i
1545 return matrix(F
, n
/2, elements
)
1549 def natural_inner_product(cls
,X
,Y
):
1551 Compute a natural inner product in this algebra directly from
1556 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1560 This gives the same answer as the slow, default method implemented
1561 in :class:`MatrixEuclideanJordanAlgebra`::
1563 sage: set_random_seed()
1564 sage: J = ComplexHermitianEJA.random_instance()
1565 sage: x,y = J.random_elements(2)
1566 sage: Xe = x.natural_representation()
1567 sage: Ye = y.natural_representation()
1568 sage: X = ComplexHermitianEJA.real_unembed(Xe)
1569 sage: Y = ComplexHermitianEJA.real_unembed(Ye)
1570 sage: expected = (X*Y).trace().real()
1571 sage: actual = ComplexHermitianEJA.natural_inner_product(Xe,Ye)
1572 sage: actual == expected
1576 return RealMatrixEuclideanJordanAlgebra
.natural_inner_product(X
,Y
)/2
1579 class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra
):
1581 The rank-n simple EJA consisting of complex Hermitian n-by-n
1582 matrices over the real numbers, the usual symmetric Jordan product,
1583 and the real-part-of-trace inner product. It has dimension `n^2` over
1588 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1592 In theory, our "field" can be any subfield of the reals::
1594 sage: ComplexHermitianEJA(2, RDF)
1595 Euclidean Jordan algebra of dimension 4 over Real Double Field
1596 sage: ComplexHermitianEJA(2, RR)
1597 Euclidean Jordan algebra of dimension 4 over Real Field with
1598 53 bits of precision
1602 The dimension of this algebra is `n^2`::
1604 sage: set_random_seed()
1605 sage: n_max = ComplexHermitianEJA._max_random_instance_size()
1606 sage: n = ZZ.random_element(1, n_max)
1607 sage: J = ComplexHermitianEJA(n)
1608 sage: J.dimension() == n^2
1611 The Jordan multiplication is what we think it is::
1613 sage: set_random_seed()
1614 sage: J = ComplexHermitianEJA.random_instance()
1615 sage: x,y = J.random_elements(2)
1616 sage: actual = (x*y).natural_representation()
1617 sage: X = x.natural_representation()
1618 sage: Y = y.natural_representation()
1619 sage: expected = (X*Y + Y*X)/2
1620 sage: actual == expected
1622 sage: J(expected) == x*y
1625 We can change the generator prefix::
1627 sage: ComplexHermitianEJA(2, prefix='z').gens()
1630 Our natural basis is normalized with respect to the natural inner
1631 product unless we specify otherwise::
1633 sage: set_random_seed()
1634 sage: J = ComplexHermitianEJA.random_instance()
1635 sage: all( b.norm() == 1 for b in J.gens() )
1638 Since our natural basis is normalized with respect to the natural
1639 inner product, and since we know that this algebra is an EJA, any
1640 left-multiplication operator's matrix will be symmetric because
1641 natural->EJA basis representation is an isometry and within the EJA
1642 the operator is self-adjoint by the Jordan axiom::
1644 sage: set_random_seed()
1645 sage: x = ComplexHermitianEJA.random_instance().random_element()
1646 sage: x.operator().matrix().is_symmetric()
1649 We can construct the (trivial) algebra of rank zero::
1651 sage: ComplexHermitianEJA(0)
1652 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1657 def _denormalized_basis(cls
, n
, field
):
1659 Returns a basis for the space of complex Hermitian n-by-n matrices.
1661 Why do we embed these? Basically, because all of numerical linear
1662 algebra assumes that you're working with vectors consisting of `n`
1663 entries from a field and scalars from the same field. There's no way
1664 to tell SageMath that (for example) the vectors contain complex
1665 numbers, while the scalar field is real.
1669 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1673 sage: set_random_seed()
1674 sage: n = ZZ.random_element(1,5)
1675 sage: field = QuadraticField(2, 'sqrt2')
1676 sage: B = ComplexHermitianEJA._denormalized_basis(n, field)
1677 sage: all( M.is_symmetric() for M in B)
1681 R
= PolynomialRing(field
, 'z')
1683 F
= field
.extension(z
**2 + 1, 'I')
1686 # This is like the symmetric case, but we need to be careful:
1688 # * We want conjugate-symmetry, not just symmetry.
1689 # * The diagonal will (as a result) be real.
1693 for j
in range(i
+1):
1694 Eij
= matrix(F
, n
, lambda k
,l
: k
==i
and l
==j
)
1696 Sij
= cls
.real_embed(Eij
)
1699 # The second one has a minus because it's conjugated.
1700 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
1702 Sij_imag
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
1705 # Since we embedded these, we can drop back to the "field" that we
1706 # started with instead of the complex extension "F".
1707 return ( s
.change_ring(field
) for s
in S
)
1710 def __init__(self
, n
, field
=AA
, **kwargs
):
1711 basis
= self
._denormalized
_basis
(n
,field
)
1712 super(ComplexHermitianEJA
,self
).__init
__(field
,
1716 self
.rank
.set_cache(n
)
1719 class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1723 Embed the n-by-n quaternion matrix ``M`` into the space of real
1724 matrices of size 4n-by-4n by first sending each quaternion entry `z
1725 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
1726 c+di],[-c + di, a-bi]]`, and then embedding those into a real
1731 sage: from mjo.eja.eja_algebra import \
1732 ....: QuaternionMatrixEuclideanJordanAlgebra
1736 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1737 sage: i,j,k = Q.gens()
1738 sage: x = 1 + 2*i + 3*j + 4*k
1739 sage: M = matrix(Q, 1, [[x]])
1740 sage: QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1746 Embedding is a homomorphism (isomorphism, in fact)::
1748 sage: set_random_seed()
1749 sage: n_max = QuaternionMatrixEuclideanJordanAlgebra._max_random_instance_size()
1750 sage: n = ZZ.random_element(n_max)
1751 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1752 sage: X = random_matrix(Q, n)
1753 sage: Y = random_matrix(Q, n)
1754 sage: Xe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X)
1755 sage: Ye = QuaternionMatrixEuclideanJordanAlgebra.real_embed(Y)
1756 sage: XYe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1761 quaternions
= M
.base_ring()
1764 raise ValueError("the matrix 'M' must be square")
1766 F
= QuadraticField(-1, 'I')
1771 t
= z
.coefficient_tuple()
1776 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
1777 [-c
+ d
*i
, a
- b
*i
]])
1778 realM
= ComplexMatrixEuclideanJordanAlgebra
.real_embed(cplxM
)
1779 blocks
.append(realM
)
1781 # We should have real entries by now, so use the realest field
1782 # we've got for the return value.
1783 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
1788 def real_unembed(M
):
1790 The inverse of _embed_quaternion_matrix().
1794 sage: from mjo.eja.eja_algebra import \
1795 ....: QuaternionMatrixEuclideanJordanAlgebra
1799 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
1800 ....: [-2, 1, -4, 3],
1801 ....: [-3, 4, 1, -2],
1802 ....: [-4, -3, 2, 1]])
1803 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(M)
1804 [1 + 2*i + 3*j + 4*k]
1808 Unembedding is the inverse of embedding::
1810 sage: set_random_seed()
1811 sage: Q = QuaternionAlgebra(QQ, -1, -1)
1812 sage: M = random_matrix(Q, 3)
1813 sage: Me = QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1814 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1820 raise ValueError("the matrix 'M' must be square")
1821 if not n
.mod(4).is_zero():
1822 raise ValueError("the matrix 'M' must be a quaternion embedding")
1824 # Use the base ring of the matrix to ensure that its entries can be
1825 # multiplied by elements of the quaternion algebra.
1826 field
= M
.base_ring()
1827 Q
= QuaternionAlgebra(field
,-1,-1)
1830 # Go top-left to bottom-right (reading order), converting every
1831 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
1834 for l
in range(n
/4):
1835 for m
in range(n
/4):
1836 submat
= ComplexMatrixEuclideanJordanAlgebra
.real_unembed(
1837 M
[4*l
:4*l
+4,4*m
:4*m
+4] )
1838 if submat
[0,0] != submat
[1,1].conjugate():
1839 raise ValueError('bad on-diagonal submatrix')
1840 if submat
[0,1] != -submat
[1,0].conjugate():
1841 raise ValueError('bad off-diagonal submatrix')
1842 z
= submat
[0,0].real()
1843 z
+= submat
[0,0].imag()*i
1844 z
+= submat
[0,1].real()*j
1845 z
+= submat
[0,1].imag()*k
1848 return matrix(Q
, n
/4, elements
)
1852 def natural_inner_product(cls
,X
,Y
):
1854 Compute a natural inner product in this algebra directly from
1859 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1863 This gives the same answer as the slow, default method implemented
1864 in :class:`MatrixEuclideanJordanAlgebra`::
1866 sage: set_random_seed()
1867 sage: J = QuaternionHermitianEJA.random_instance()
1868 sage: x,y = J.random_elements(2)
1869 sage: Xe = x.natural_representation()
1870 sage: Ye = y.natural_representation()
1871 sage: X = QuaternionHermitianEJA.real_unembed(Xe)
1872 sage: Y = QuaternionHermitianEJA.real_unembed(Ye)
1873 sage: expected = (X*Y).trace().coefficient_tuple()[0]
1874 sage: actual = QuaternionHermitianEJA.natural_inner_product(Xe,Ye)
1875 sage: actual == expected
1879 return RealMatrixEuclideanJordanAlgebra
.natural_inner_product(X
,Y
)/4
1882 class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra
):
1884 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
1885 matrices, the usual symmetric Jordan product, and the
1886 real-part-of-trace inner product. It has dimension `2n^2 - n` over
1891 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1895 In theory, our "field" can be any subfield of the reals::
1897 sage: QuaternionHermitianEJA(2, RDF)
1898 Euclidean Jordan algebra of dimension 6 over Real Double Field
1899 sage: QuaternionHermitianEJA(2, RR)
1900 Euclidean Jordan algebra of dimension 6 over Real Field with
1901 53 bits of precision
1905 The dimension of this algebra is `2*n^2 - n`::
1907 sage: set_random_seed()
1908 sage: n_max = QuaternionHermitianEJA._max_random_instance_size()
1909 sage: n = ZZ.random_element(1, n_max)
1910 sage: J = QuaternionHermitianEJA(n)
1911 sage: J.dimension() == 2*(n^2) - n
1914 The Jordan multiplication is what we think it is::
1916 sage: set_random_seed()
1917 sage: J = QuaternionHermitianEJA.random_instance()
1918 sage: x,y = J.random_elements(2)
1919 sage: actual = (x*y).natural_representation()
1920 sage: X = x.natural_representation()
1921 sage: Y = y.natural_representation()
1922 sage: expected = (X*Y + Y*X)/2
1923 sage: actual == expected
1925 sage: J(expected) == x*y
1928 We can change the generator prefix::
1930 sage: QuaternionHermitianEJA(2, prefix='a').gens()
1931 (a0, a1, a2, a3, a4, a5)
1933 Our natural basis is normalized with respect to the natural inner
1934 product unless we specify otherwise::
1936 sage: set_random_seed()
1937 sage: J = QuaternionHermitianEJA.random_instance()
1938 sage: all( b.norm() == 1 for b in J.gens() )
1941 Since our natural basis is normalized with respect to the natural
1942 inner product, and since we know that this algebra is an EJA, any
1943 left-multiplication operator's matrix will be symmetric because
1944 natural->EJA basis representation is an isometry and within the EJA
1945 the operator is self-adjoint by the Jordan axiom::
1947 sage: set_random_seed()
1948 sage: x = QuaternionHermitianEJA.random_instance().random_element()
1949 sage: x.operator().matrix().is_symmetric()
1952 We can construct the (trivial) algebra of rank zero::
1954 sage: QuaternionHermitianEJA(0)
1955 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1959 def _denormalized_basis(cls
, n
, field
):
1961 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
1963 Why do we embed these? Basically, because all of numerical
1964 linear algebra assumes that you're working with vectors consisting
1965 of `n` entries from a field and scalars from the same field. There's
1966 no way to tell SageMath that (for example) the vectors contain
1967 complex numbers, while the scalar field is real.
1971 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1975 sage: set_random_seed()
1976 sage: n = ZZ.random_element(1,5)
1977 sage: B = QuaternionHermitianEJA._denormalized_basis(n,QQ)
1978 sage: all( M.is_symmetric() for M in B )
1982 Q
= QuaternionAlgebra(QQ
,-1,-1)
1985 # This is like the symmetric case, but we need to be careful:
1987 # * We want conjugate-symmetry, not just symmetry.
1988 # * The diagonal will (as a result) be real.
1992 for j
in range(i
+1):
1993 Eij
= matrix(Q
, n
, lambda k
,l
: k
==i
and l
==j
)
1995 Sij
= cls
.real_embed(Eij
)
1998 # The second, third, and fourth ones have a minus
1999 # because they're conjugated.
2000 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
2002 Sij_I
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
2004 Sij_J
= cls
.real_embed(J
*Eij
- J
*Eij
.transpose())
2006 Sij_K
= cls
.real_embed(K
*Eij
- K
*Eij
.transpose())
2009 # Since we embedded these, we can drop back to the "field" that we
2010 # started with instead of the quaternion algebra "Q".
2011 return ( s
.change_ring(field
) for s
in S
)
2014 def __init__(self
, n
, field
=AA
, **kwargs
):
2015 basis
= self
._denormalized
_basis
(n
,field
)
2016 super(QuaternionHermitianEJA
,self
).__init
__(field
,
2020 self
.rank
.set_cache(n
)
2023 class HadamardEJA(RationalBasisEuclideanJordanAlgebra
):
2025 Return the Euclidean Jordan Algebra corresponding to the set
2026 `R^n` under the Hadamard product.
2028 Note: this is nothing more than the Cartesian product of ``n``
2029 copies of the spin algebra. Once Cartesian product algebras
2030 are implemented, this can go.
2034 sage: from mjo.eja.eja_algebra import HadamardEJA
2038 This multiplication table can be verified by hand::
2040 sage: J = HadamardEJA(3)
2041 sage: e0,e1,e2 = J.gens()
2057 We can change the generator prefix::
2059 sage: HadamardEJA(3, prefix='r').gens()
2063 def __init__(self
, n
, field
=AA
, **kwargs
):
2064 V
= VectorSpace(field
, n
)
2065 mult_table
= [ [ V
.gen(i
)*(i
== j
) for j
in range(n
) ]
2068 super(HadamardEJA
, self
).__init
__(field
,
2072 self
.rank
.set_cache(n
)
2073 self
.one
.set_cache( sum(self
.gens()) )
2076 def _max_random_instance_size():
2079 def inner_product(self
, x
, y
):
2081 Faster to reimplement than to use natural representations.
2085 sage: from mjo.eja.eja_algebra import HadamardEJA
2089 Ensure that this is the usual inner product for the algebras
2092 sage: set_random_seed()
2093 sage: J = HadamardEJA.random_instance()
2094 sage: x,y = J.random_elements(2)
2095 sage: X = x.natural_representation()
2096 sage: Y = y.natural_representation()
2097 sage: x.inner_product(y) == J.natural_inner_product(X,Y)
2101 return x
.to_vector().inner_product(y
.to_vector())
2104 class BilinearFormEJA(RationalBasisEuclideanJordanAlgebra
):
2106 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2107 with the half-trace inner product and jordan product ``x*y =
2108 (<Bx,y>,y_bar>, x0*y_bar + y0*x_bar)`` where `B = 1 \times B22` is
2109 a symmetric positive-definite "bilinear form" matrix. Its
2110 dimension is the size of `B`, and it has rank two in dimensions
2111 larger than two. It reduces to the ``JordanSpinEJA`` when `B` is
2112 the identity matrix of order ``n``.
2114 We insist that the one-by-one upper-left identity block of `B` be
2115 passed in as well so that we can be passed a matrix of size zero
2116 to construct a trivial algebra.
2120 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
2121 ....: JordanSpinEJA)
2125 When no bilinear form is specified, the identity matrix is used,
2126 and the resulting algebra is the Jordan spin algebra::
2128 sage: B = matrix.identity(AA,3)
2129 sage: J0 = BilinearFormEJA(B)
2130 sage: J1 = JordanSpinEJA(3)
2131 sage: J0.multiplication_table() == J0.multiplication_table()
2134 An error is raised if the matrix `B` does not correspond to a
2135 positive-definite bilinear form::
2137 sage: B = matrix.random(QQ,2,3)
2138 sage: J = BilinearFormEJA(B)
2139 Traceback (most recent call last):
2141 ValueError: bilinear form is not positive-definite
2142 sage: B = matrix.zero(QQ,3)
2143 sage: J = BilinearFormEJA(B)
2144 Traceback (most recent call last):
2146 ValueError: bilinear form is not positive-definite
2150 We can create a zero-dimensional algebra::
2152 sage: B = matrix.identity(AA,0)
2153 sage: J = BilinearFormEJA(B)
2157 We can check the multiplication condition given in the Jordan, von
2158 Neumann, and Wigner paper (and also discussed on my "On the
2159 symmetry..." paper). Note that this relies heavily on the standard
2160 choice of basis, as does anything utilizing the bilinear form matrix::
2162 sage: set_random_seed()
2163 sage: n = ZZ.random_element(5)
2164 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2165 sage: B11 = matrix.identity(QQ,1)
2166 sage: B22 = M.transpose()*M
2167 sage: B = block_matrix(2,2,[ [B11,0 ],
2169 sage: J = BilinearFormEJA(B)
2170 sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
2171 sage: V = J.vector_space()
2172 sage: sis = [ J.from_vector(V([0] + (M.inverse()*ei).list()))
2173 ....: for ei in eis ]
2174 sage: actual = [ sis[i]*sis[j]
2175 ....: for i in range(n-1)
2176 ....: for j in range(n-1) ]
2177 sage: expected = [ J.one() if i == j else J.zero()
2178 ....: for i in range(n-1)
2179 ....: for j in range(n-1) ]
2180 sage: actual == expected
2183 def __init__(self
, B
, field
=AA
, **kwargs
):
2187 if not B
.is_positive_definite():
2188 raise ValueError("bilinear form is not positive-definite")
2190 V
= VectorSpace(field
, n
)
2191 mult_table
= [[V
.zero() for j
in range(n
)] for i
in range(n
)]
2200 z0
= (B
*x
).inner_product(y
)
2201 zbar
= y0
*xbar
+ x0
*ybar
2202 z
= V([z0
] + zbar
.list())
2203 mult_table
[i
][j
] = z
2205 # The rank of this algebra is two, unless we're in a
2206 # one-dimensional ambient space (because the rank is bounded
2207 # by the ambient dimension).
2208 super(BilinearFormEJA
, self
).__init
__(field
,
2212 self
.rank
.set_cache(min(n
,2))
2215 self
.one
.set_cache( self
.zero() )
2217 self
.one
.set_cache( self
.monomial(0) )
2220 def _max_random_instance_size():
2224 def random_instance(cls
, field
=AA
, **kwargs
):
2226 Return a random instance of this algebra.
2228 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2230 B
= matrix
.identity(field
, n
)
2231 return cls(B
, field
, **kwargs
)
2233 B11
= matrix
.identity(field
,1)
2234 M
= matrix
.random(field
, n
-1)
2235 I
= matrix
.identity(field
, n
-1)
2236 alpha
= field
.zero()
2237 while alpha
.is_zero():
2238 alpha
= field
.random_element().abs()
2239 B22
= M
.transpose()*M
+ alpha
*I
2241 from sage
.matrix
.special
import block_matrix
2242 B
= block_matrix(2,2, [ [B11
, ZZ(0) ],
2245 return cls(B
, field
, **kwargs
)
2247 def inner_product(self
, x
, y
):
2249 Half of the trace inner product.
2251 This is defined so that the special case of the Jordan spin
2252 algebra gets the usual inner product.
2256 sage: from mjo.eja.eja_algebra import BilinearFormEJA
2260 Ensure that this is one-half of the trace inner-product when
2261 the algebra isn't just the reals (when ``n`` isn't one). This
2262 is in Faraut and Koranyi, and also my "On the symmetry..."
2265 sage: set_random_seed()
2266 sage: J = BilinearFormEJA.random_instance()
2267 sage: n = J.dimension()
2268 sage: x = J.random_element()
2269 sage: y = J.random_element()
2270 sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
2274 return (self
._B
*x
.to_vector()).inner_product(y
.to_vector())
2277 class JordanSpinEJA(BilinearFormEJA
):
2279 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2280 with the usual inner product and jordan product ``x*y =
2281 (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2286 sage: from mjo.eja.eja_algebra import JordanSpinEJA
2290 This multiplication table can be verified by hand::
2292 sage: J = JordanSpinEJA(4)
2293 sage: e0,e1,e2,e3 = J.gens()
2309 We can change the generator prefix::
2311 sage: JordanSpinEJA(2, prefix='B').gens()
2316 Ensure that we have the usual inner product on `R^n`::
2318 sage: set_random_seed()
2319 sage: J = JordanSpinEJA.random_instance()
2320 sage: x,y = J.random_elements(2)
2321 sage: X = x.natural_representation()
2322 sage: Y = y.natural_representation()
2323 sage: x.inner_product(y) == J.natural_inner_product(X,Y)
2327 def __init__(self
, n
, field
=AA
, **kwargs
):
2328 # This is a special case of the BilinearFormEJA with the identity
2329 # matrix as its bilinear form.
2330 B
= matrix
.identity(field
, n
)
2331 super(JordanSpinEJA
, self
).__init
__(B
, field
, **kwargs
)
2334 def random_instance(cls
, field
=AA
, **kwargs
):
2336 Return a random instance of this type of algebra.
2338 Needed here to override the implementation for ``BilinearFormEJA``.
2340 n
= ZZ
.random_element(cls
._max
_random
_instance
_size
() + 1)
2341 return cls(n
, field
, **kwargs
)
2344 class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra
):
2346 The trivial Euclidean Jordan algebra consisting of only a zero element.
2350 sage: from mjo.eja.eja_algebra import TrivialEJA
2354 sage: J = TrivialEJA()
2361 sage: 7*J.one()*12*J.one()
2363 sage: J.one().inner_product(J.one())
2365 sage: J.one().norm()
2367 sage: J.one().subalgebra_generated_by()
2368 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2373 def __init__(self
, field
=AA
, **kwargs
):
2375 super(TrivialEJA
, self
).__init
__(field
,
2379 # The rank is zero using my definition, namely the dimension of the
2380 # largest subalgebra generated by any element.
2381 self
.rank
.set_cache(0)
2382 self
.one
.set_cache( self
.zero() )
2385 def random_instance(cls
, field
=AA
, **kwargs
):
2386 # We don't take a "size" argument so the superclass method is
2387 # inappropriate for us.
2388 return cls(field
, **kwargs
)
2390 class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra
):
2392 The external (orthogonal) direct sum of two other Euclidean Jordan
2393 algebras. Essentially the Cartesian product of its two factors.
2394 Every Euclidean Jordan algebra decomposes into an orthogonal
2395 direct sum of simple Euclidean Jordan algebras, so no generality
2396 is lost by providing only this construction.
2400 sage: from mjo.eja.eja_algebra import (HadamardEJA,
2401 ....: RealSymmetricEJA,
2406 sage: J1 = HadamardEJA(2)
2407 sage: J2 = RealSymmetricEJA(3)
2408 sage: J = DirectSumEJA(J1,J2)
2415 def __init__(self
, J1
, J2
, field
=AA
, **kwargs
):
2416 self
._factors
= (J1
, J2
)
2420 V
= VectorSpace(field
, n
)
2421 mult_table
= [ [ V
.zero() for j
in range(n
) ]
2425 p
= (J1
.monomial(i
)*J1
.monomial(j
)).to_vector()
2426 mult_table
[i
][j
] = V(p
.list() + [field
.zero()]*n2
)
2430 p
= (J2
.monomial(i
)*J2
.monomial(j
)).to_vector()
2431 mult_table
[n1
+i
][n1
+j
] = V([field
.zero()]*n1
+ p
.list())
2433 super(DirectSumEJA
, self
).__init
__(field
,
2437 self
.rank
.set_cache(J1
.rank() + J2
.rank())
2442 Return the pair of this algebra's factors.
2446 sage: from mjo.eja.eja_algebra import (HadamardEJA,
2447 ....: JordanSpinEJA,
2452 sage: J1 = HadamardEJA(2,QQ)
2453 sage: J2 = JordanSpinEJA(3,QQ)
2454 sage: J = DirectSumEJA(J1,J2)
2456 (Euclidean Jordan algebra of dimension 2 over Rational Field,
2457 Euclidean Jordan algebra of dimension 3 over Rational Field)
2460 return self
._factors
2462 def projections(self
):
2464 Return a pair of projections onto this algebra's factors.
2468 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
2469 ....: ComplexHermitianEJA,
2474 sage: J1 = JordanSpinEJA(2)
2475 sage: J2 = ComplexHermitianEJA(2)
2476 sage: J = DirectSumEJA(J1,J2)
2477 sage: (pi_left, pi_right) = J.projections()
2478 sage: J.one().to_vector()
2480 sage: pi_left(J.one()).to_vector()
2482 sage: pi_right(J.one()).to_vector()
2486 (J1
,J2
) = self
.factors()
2488 pi_left
= lambda x
: J1
.from_vector(x
.to_vector()[:n
])
2489 pi_right
= lambda x
: J2
.from_vector(x
.to_vector()[n
:])
2490 return (pi_left
, pi_right
)
2492 def inclusions(self
):
2494 Return the pair of inclusion maps from our factors into us.
2498 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
2499 ....: RealSymmetricEJA,
2504 sage: J1 = JordanSpinEJA(3)
2505 sage: J2 = RealSymmetricEJA(2)
2506 sage: J = DirectSumEJA(J1,J2)
2507 sage: (iota_left, iota_right) = J.inclusions()
2508 sage: iota_left(J1.zero()) == J.zero()
2510 sage: iota_right(J2.zero()) == J.zero()
2512 sage: J1.one().to_vector()
2514 sage: iota_left(J1.one()).to_vector()
2516 sage: J2.one().to_vector()
2518 sage: iota_right(J2.one()).to_vector()
2520 sage: J.one().to_vector()
2524 (J1
,J2
) = self
.factors()
2526 V_basis
= self
.vector_space().basis()
2527 I1
= matrix
.column(self
.base_ring(), V_basis
[:n
])
2528 I2
= matrix
.column(self
.base_ring(), V_basis
[n
:])
2529 iota_left
= lambda x
: self
.from_vector(I1
*x
.to_vector())
2530 iota_right
= lambda x
: self
.from_vector(I2
*+x
.to_vector())
2531 return (iota_left
, iota_right
)
2533 def inner_product(self
, x
, y
):
2535 The standard Cartesian inner-product.
2537 We project ``x`` and ``y`` onto our factors, and add up the
2538 inner-products from the subalgebras.
2543 sage: from mjo.eja.eja_algebra import (HadamardEJA,
2544 ....: QuaternionHermitianEJA,
2549 sage: J1 = HadamardEJA(3)
2550 sage: J2 = QuaternionHermitianEJA(2,QQ,normalize_basis=False)
2551 sage: J = DirectSumEJA(J1,J2)
2556 sage: x1.inner_product(x2)
2558 sage: y1.inner_product(y2)
2560 sage: J.one().inner_product(J.one())
2564 (pi_left
, pi_right
) = self
.projections()
2570 return (x1
.inner_product(y1
) + x2
.inner_product(y2
))