X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=mjo%2Feja%2Feja_algebra.py;h=359b7404a6fe5d003151e1e82d6c78090b0ea9f6;hb=3c7644ecfe369b6f83aa707b87d7a1f9aa246e27;hp=fb840edb93ed564c7372eacac2e3b9913a4a2350;hpb=f9e1e977040c2c3b1019a0bcc5776384a13e96c4;p=sage.d.git diff --git a/mjo/eja/eja_algebra.py b/mjo/eja/eja_algebra.py index fb840ed..e32bc24 100644 --- a/mjo/eja/eja_algebra.py +++ b/mjo/eja/eja_algebra.py @@ -3,73 +3,89 @@ Euclidean Jordan Algebras. These are formally-real Jordan Algebras; specifically those where u^2 + v^2 = 0 implies that u = v = 0. They are used in optimization, and have some additional nice methods beyond what can be supported in a general Jordan Algebra. + + +SETUP:: + + sage: from mjo.eja.eja_algebra import random_eja + +EXAMPLES:: + + sage: random_eja() + Euclidean Jordan algebra of dimension... + """ -from sage.algebras.finite_dimensional_algebras.finite_dimensional_algebra import FiniteDimensionalAlgebra +from itertools import repeat + from sage.algebras.quatalg.quaternion_algebra import QuaternionAlgebra -from sage.categories.finite_dimensional_algebras_with_basis import FiniteDimensionalAlgebrasWithBasis +from sage.categories.magmatic_algebras import MagmaticAlgebras +from sage.combinat.free_module import CombinatorialFreeModule from sage.matrix.constructor import matrix +from sage.matrix.matrix_space import MatrixSpace from sage.misc.cachefunc import cached_method -from sage.misc.prandom import choice -from sage.modules.free_module import VectorSpace -from sage.rings.integer_ring import ZZ -from sage.rings.number_field.number_field import QuadraticField -from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing -from sage.rings.rational_field import QQ -from sage.structure.element import is_Matrix -from sage.structure.category_object import normalize_names - +from sage.misc.table import table +from sage.modules.free_module import FreeModule, VectorSpace +from sage.rings.all import (ZZ, QQ, AA, QQbar, RR, RLF, CLF, + PolynomialRing, + QuadraticField) from mjo.eja.eja_element import FiniteDimensionalEuclideanJordanAlgebraElement -from mjo.eja.eja_utils import _vec2mat, _mat2vec +from mjo.eja.eja_operator import FiniteDimensionalEuclideanJordanAlgebraOperator +from mjo.eja.eja_utils import _mat2vec -class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra): - @staticmethod - def __classcall_private__(cls, - field, - mult_table, - rank, - names='e', - assume_associative=False, - category=None, - natural_basis=None): - n = len(mult_table) - mult_table = [b.base_extend(field) for b in mult_table] - for b in mult_table: - b.set_immutable() - if not (is_Matrix(b) and b.dimensions() == (n, n)): - raise ValueError("input is not a multiplication table") - mult_table = tuple(mult_table) - - cat = FiniteDimensionalAlgebrasWithBasis(field) - cat.or_subcategory(category) - if assume_associative: - cat = cat.Associative() - - names = normalize_names(n, names) - - fda = super(FiniteDimensionalEuclideanJordanAlgebra, cls) - return fda.__classcall__(cls, - field, - mult_table, - rank, - assume_associative=assume_associative, - names=names, - category=cat, - natural_basis=natural_basis) +class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): + r""" + The lowest-level class for representing a Euclidean Jordan algebra. + """ + def _coerce_map_from_base_ring(self): + """ + Disable the map from the base ring into the algebra. + + Performing a nonsense conversion like this automatically + is counterpedagogical. The fallback is to try the usual + element constructor, which should also fail. + + SETUP:: + + sage: from mjo.eja.eja_algebra import random_eja + TESTS:: + + sage: set_random_seed() + sage: J = random_eja() + sage: J(1) + Traceback (most recent call last): + ... + ValueError: not an element of this algebra + + """ + return None def __init__(self, field, - mult_table, - rank, - names='e', - assume_associative=False, + multiplication_table, + inner_product_table, + prefix='e', category=None, - natural_basis=None): + matrix_basis=None, + check_field=True, + check_axioms=True): """ + INPUT: + + * field -- the scalar field for this algebra (must be real) + + * multiplication_table -- the multiplication table for this + algebra's implicit basis. Only the lower-triangular portion + of the table is used, since the multiplication is assumed + to be commutative. + SETUP:: - sage: from mjo.eja.eja_algebra import random_eja + sage: from mjo.eja.eja_algebra import ( + ....: FiniteDimensionalEuclideanJordanAlgebra, + ....: JordanSpinEJA, + ....: random_eja) EXAMPLES: @@ -77,173 +93,340 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra): sage: set_random_seed() sage: J = random_eja() - sage: x = J.random_element() - sage: y = J.random_element() + sage: x,y = J.random_elements(2) sage: x*y == y*x True + An error is raised if the Jordan product is not commutative:: + + sage: JP = ((1,2),(0,0)) + sage: IP = ((1,0),(0,1)) + sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,JP,IP) + Traceback (most recent call last): + ... + ValueError: Jordan product is not commutative + + An error is raised if the inner-product is not commutative:: + + sage: JP = ((1,0),(0,1)) + sage: IP = ((1,2),(0,0)) + sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,JP,IP) + Traceback (most recent call last): + ... + ValueError: inner-product is not commutative + + TESTS: + + The ``field`` we're given must be real with ``check_field=True``:: + + sage: JordanSpinEJA(2, field=QQbar) + Traceback (most recent call last): + ... + ValueError: scalar field is not real + sage: JordanSpinEJA(2, field=QQbar, check_field=False) + Euclidean Jordan algebra of dimension 2 over Algebraic Field + + The multiplication table must be square with ``check_axioms=True``:: + + sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,((),()),((1,),)) + Traceback (most recent call last): + ... + ValueError: multiplication table is not square + + The multiplication and inner-product tables must be the same + size (and in particular, the inner-product table must also be + square) with ``check_axioms=True``:: + + sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,((1,),),(())) + Traceback (most recent call last): + ... + ValueError: multiplication and inner-product tables are + different sizes + sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,((1,),),((1,2),)) + Traceback (most recent call last): + ... + ValueError: multiplication and inner-product tables are + different sizes + """ - self._rank = rank - self._natural_basis = natural_basis - self._multiplication_table = mult_table + if check_field: + if not field.is_subring(RR): + # Note: this does return true for the real algebraic + # field, the rationals, and any quadratic field where + # we've specified a real embedding. + raise ValueError("scalar field is not real") + + + # The multiplication and inner-product tables should be square + # if the user wants us to verify them. And we verify them as + # soon as possible, because we want to exploit their symmetry. + n = len(multiplication_table) + if check_axioms: + if not all( len(l) == n for l in multiplication_table ): + raise ValueError("multiplication table is not square") + + # If the multiplication table is square, we can check if + # the inner-product table is square by comparing it to the + # multiplication table's dimensions. + msg = "multiplication and inner-product tables are different sizes" + if not len(inner_product_table) == n: + raise ValueError(msg) + + if not all( len(l) == n for l in inner_product_table ): + raise ValueError(msg) + + # Check commutativity of the Jordan product (symmetry of + # the multiplication table) and the commutativity of the + # inner-product (symmetry of the inner-product table) + # first if we're going to check them at all.. This has to + # be done before we define product_on_basis(), because + # that method assumes that self._multiplication_table is + # symmetric. And it has to be done before we build + # self._inner_product_matrix, because the process used to + # construct it assumes symmetry as well. + if not all( multiplication_table[j][i] + == multiplication_table[i][j] + for i in range(n) + for j in range(i+1) ): + raise ValueError("Jordan product is not commutative") + + if not all( inner_product_table[j][i] + == inner_product_table[i][j] + for i in range(n) + for j in range(i+1) ): + raise ValueError("inner-product is not commutative") + + self._matrix_basis = matrix_basis + + if category is None: + category = MagmaticAlgebras(field).FiniteDimensional() + category = category.WithBasis().Unital() + fda = super(FiniteDimensionalEuclideanJordanAlgebra, self) fda.__init__(field, - mult_table, - names=names, + range(n), + prefix=prefix, category=category) - - - def _repr_(self): + self.print_options(bracket='') + + # The multiplication table we're given is necessarily in terms + # of vectors, because we don't have an algebra yet for + # anything to be an element of. However, it's faster in the + # long run to have the multiplication table be in terms of + # algebra elements. We do this after calling the superclass + # constructor so that from_vector() knows what to do. + # + # Note: we take advantage of symmetry here, and only store + # the lower-triangular portion of the table. + self._multiplication_table = [ [ self.vector_space().zero() + for j in range(i+1) ] + for i in range(n) ] + + for i in range(n): + for j in range(i+1): + elt = self.from_vector(multiplication_table[i][j]) + self._multiplication_table[i][j] = elt + + self._multiplication_table = tuple(map(tuple, self._multiplication_table)) + + # Save our inner product as a matrix, since the efficiency of + # matrix multiplication will usually outweigh the fact that we + # have to store a redundant upper- or lower-triangular part. + # Pre-cache the fact that these are Hermitian (real symmetric, + # in fact) in case some e.g. matrix multiplication routine can + # take advantage of it. + ip_matrix_constructor = lambda i,j: inner_product_table[i][j] if j <= i else inner_product_table[j][i] + self._inner_product_matrix = matrix(field, n, ip_matrix_constructor) + self._inner_product_matrix._cache = {'hermitian': True} + self._inner_product_matrix.set_immutable() + + if check_axioms: + if not self._is_jordanian(): + raise ValueError("Jordan identity does not hold") + if not self._inner_product_is_associative(): + raise ValueError("inner product is not associative") + + def _element_constructor_(self, elt): """ - Return a string representation of ``self``. + Construct an element of this algebra from its vector or matrix + representation. + + This gets called only after the parent element _call_ method + fails to find a coercion for the argument. SETUP:: - sage: from mjo.eja.eja_algebra import JordanSpinEJA + sage: from mjo.eja.eja_algebra import (JordanSpinEJA, + ....: HadamardEJA, + ....: RealSymmetricEJA) + + EXAMPLES: + + The identity in `S^n` is converted to the identity in the EJA:: + + sage: J = RealSymmetricEJA(3) + sage: I = matrix.identity(QQ,3) + sage: J(I) == J.one() + True + + This skew-symmetric matrix can't be represented in the EJA:: + + sage: J = RealSymmetricEJA(3) + sage: A = matrix(QQ,3, lambda i,j: i-j) + sage: J(A) + Traceback (most recent call last): + ... + ValueError: not an element of this algebra TESTS: - Ensure that it says what we think it says:: + Ensure that we can convert any element of the two non-matrix + simple algebras (whose matrix representations are columns) + back and forth faithfully:: - sage: JordanSpinEJA(2, field=QQ) - Euclidean Jordan algebra of degree 2 over Rational Field - sage: JordanSpinEJA(3, field=RDF) - Euclidean Jordan algebra of degree 3 over Real Double Field + sage: set_random_seed() + sage: J = HadamardEJA.random_instance() + sage: x = J.random_element() + sage: J(x.to_vector().column()) == x + True + sage: J = JordanSpinEJA.random_instance() + sage: x = J.random_element() + sage: J(x.to_vector().column()) == x + True """ - fmt = "Euclidean Jordan algebra of degree {} over {}" - return fmt.format(self.degree(), self.base_ring()) + msg = "not an element of this algebra" + if elt == 0: + # The superclass implementation of random_element() + # needs to be able to coerce "0" into the algebra. + return self.zero() + elif elt in self.base_ring(): + # Ensure that no base ring -> algebra coercion is performed + # by this method. There's some stupidity in sage that would + # otherwise propagate to this method; for example, sage thinks + # that the integer 3 belongs to the space of 2-by-2 matrices. + raise ValueError(msg) + + if elt not in self.matrix_space(): + raise ValueError(msg) + + # Thanks for nothing! Matrix spaces aren't vector spaces in + # Sage, so we have to figure out its matrix-basis coordinates + # ourselves. We use the basis space's ring instead of the + # element's ring because the basis space might be an algebraic + # closure whereas the base ring of the 3-by-3 identity matrix + # could be QQ instead of QQbar. + # + # We pass check=False because the matrix basis is "guaranteed" + # to be linearly independent... right? Ha ha. + V = VectorSpace(self.base_ring(), elt.nrows()*elt.ncols()) + W = V.span_of_basis( (_mat2vec(s) for s in self.matrix_basis()), + check=False) + + try: + coords = W.coordinate_vector(_mat2vec(elt)) + except ArithmeticError: # vector is not in free module + raise ValueError(msg) + return self.from_vector(coords) - def _a_regular_element(self): + def _repr_(self): """ - Guess a regular element. Needed to compute the basis for our - characteristic polynomial coefficients. + Return a string representation of ``self``. SETUP:: - sage: from mjo.eja.eja_algebra import random_eja + sage: from mjo.eja.eja_algebra import JordanSpinEJA TESTS: - Ensure that this hacky method succeeds for every algebra that we - know how to construct:: + Ensure that it says what we think it says:: - sage: set_random_seed() - sage: J = random_eja() - sage: J._a_regular_element().is_regular() - True + sage: JordanSpinEJA(2, field=AA) + Euclidean Jordan algebra of dimension 2 over Algebraic Real Field + sage: JordanSpinEJA(3, field=RDF) + Euclidean Jordan algebra of dimension 3 over Real Double Field """ - gs = self.gens() - z = self.sum( (i+1)*gs[i] for i in range(len(gs)) ) - if not z.is_regular(): - raise ValueError("don't know a regular element") - return z + fmt = "Euclidean Jordan algebra of dimension {} over {}" + return fmt.format(self.dimension(), self.base_ring()) + + def product_on_basis(self, i, j): + # We only stored the lower-triangular portion of the + # multiplication table. + if j <= i: + return self._multiplication_table[i][j] + else: + return self._multiplication_table[j][i] + def _is_commutative(self): + r""" + Whether or not this algebra's multiplication table is commutative. - @cached_method - def _charpoly_basis_space(self): - """ - Return the vector space spanned by the basis used in our - characteristic polynomial coefficients. This is used not only to - compute those coefficients, but also any time we need to - evaluate the coefficients (like when we compute the trace or - determinant). + This method should of course always return ``True``, unless + this algebra was constructed with ``check_axioms=False`` and + passed an invalid multiplication table. """ - z = self._a_regular_element() - V = self.vector_space() - V1 = V.span_of_basis( (z**k).vector() for k in range(self.rank()) ) - b = (V1.basis() + V1.complement().basis()) - return V.span_of_basis(b) - - - @cached_method - def _charpoly_coeff(self, i): - """ - Return the coefficient polynomial "a_{i}" of this algebra's - general characteristic polynomial. - - Having this be a separate cached method lets us compute and - store the trace/determinant (a_{r-1} and a_{0} respectively) - separate from the entire characteristic polynomial. - """ - (A_of_x, x, xr, detA) = self._charpoly_matrix_system() - R = A_of_x.base_ring() - if i >= self.rank(): - # Guaranteed by theory - return R.zero() - - # Danger: the in-place modification is done for performance - # reasons (reconstructing a matrix with huge polynomial - # entries is slow), but I don't know how cached_method works, - # so it's highly possible that we're modifying some global - # list variable by reference, here. In other words, you - # probably shouldn't call this method twice on the same - # algebra, at the same time, in two threads - Ai_orig = A_of_x.column(i) - A_of_x.set_column(i,xr) - numerator = A_of_x.det() - A_of_x.set_column(i,Ai_orig) - - # We're relying on the theory here to ensure that each a_i is - # indeed back in R, and the added negative signs are to make - # the whole charpoly expression sum to zero. - return R(-numerator/detA) - - - @cached_method - def _charpoly_matrix_system(self): + return all( self.product_on_basis(i,j) == self.product_on_basis(i,j) + for i in range(self.dimension()) + for j in range(self.dimension()) ) + + def _is_jordanian(self): + r""" + Whether or not this algebra's multiplication table respects the + Jordan identity `(x^{2})(xy) = x(x^{2}y)`. + + We only check one arrangement of `x` and `y`, so for a + ``True`` result to be truly true, you should also check + :meth:`_is_commutative`. This method should of course always + return ``True``, unless this algebra was constructed with + ``check_axioms=False`` and passed an invalid multiplication table. """ - Compute the matrix whose entries A_ij are polynomials in - X1,...,XN, the vector ``x`` of variables X1,...,XN, the vector - corresponding to `x^r` and the determinent of the matrix A = - [A_ij]. In other words, all of the fixed (cachable) data needed - to compute the coefficients of the characteristic polynomial. + return all( (self.monomial(i)**2)*(self.monomial(i)*self.monomial(j)) + == + (self.monomial(i))*((self.monomial(i)**2)*self.monomial(j)) + for i in range(self.dimension()) + for j in range(self.dimension()) ) + + def _inner_product_is_associative(self): + r""" + Return whether or not this algebra's inner product `B` is + associative; that is, whether or not `B(xy,z) = B(x,yz)`. + + This method should of course always return ``True``, unless + this algebra was constructed with ``check_axioms=False`` and + passed an invalid multiplication table. """ - r = self.rank() - n = self.dimension() - - # Construct a new algebra over a multivariate polynomial ring... - names = ['X' + str(i) for i in range(1,n+1)] - R = PolynomialRing(self.base_ring(), names) - J = FiniteDimensionalEuclideanJordanAlgebra(R, - self._multiplication_table, - r) - idmat = matrix.identity(J.base_ring(), n) + # Used to check whether or not something is zero in an inexact + # ring. This number is sufficient to allow the construction of + # QuaternionHermitianEJA(2, field=RDF) with check_axioms=True. + epsilon = 1e-16 - W = self._charpoly_basis_space() - W = W.change_ring(R.fraction_field()) + for i in range(self.dimension()): + for j in range(self.dimension()): + for k in range(self.dimension()): + x = self.monomial(i) + y = self.monomial(j) + z = self.monomial(k) + diff = (x*y).inner_product(z) - x.inner_product(y*z) - # Starting with the standard coordinates x = (X1,X2,...,Xn) - # and then converting the entries to W-coordinates allows us - # to pass in the standard coordinates to the charpoly and get - # back the right answer. Specifically, with x = (X1,X2,...,Xn), - # we have - # - # W.coordinates(x^2) eval'd at (standard z-coords) - # = - # W-coords of (z^2) - # = - # W-coords of (standard coords of x^2 eval'd at std-coords of z) - # - # We want the middle equivalent thing in our matrix, but use - # the first equivalent thing instead so that we can pass in - # standard coordinates. - x = J(W(R.gens())) - l1 = [matrix.column(W.coordinates((x**k).vector())) for k in range(r)] - l2 = [idmat.column(k-1).column() for k in range(r+1, n+1)] - A_of_x = matrix.block(R, 1, n, (l1 + l2)) - xr = W.coordinates((x**r).vector()) - return (A_of_x, x, xr, A_of_x.det()) + if self.base_ring().is_exact(): + if diff != 0: + return False + else: + if diff.abs() > epsilon: + return False + return True @cached_method - def characteristic_polynomial(self): + def characteristic_polynomial_of(self): """ - Return a characteristic polynomial that works for all elements - of this algebra. + Return the algebra's "characteristic polynomial of" function, + which is itself a multivariate polynomial that, when evaluated + at the coordinates of some algebra element, returns that + element's characteristic polynomial. The resulting polynomial has `n+1` variables, where `n` is the dimension of this algebra. The first `n` variables correspond to @@ -255,7 +438,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra): SETUP:: - sage: from mjo.eja.eja_algebra import JordanSpinEJA + sage: from mjo.eja.eja_algebra import JordanSpinEJA, TrivialEJA EXAMPLES: @@ -263,42 +446,64 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra): Alizadeh, Example 11.11:: sage: J = JordanSpinEJA(3) - sage: p = J.characteristic_polynomial(); p + sage: p = J.characteristic_polynomial_of(); p X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2 - sage: xvec = J.one().vector() + sage: xvec = J.one().to_vector() sage: p(*xvec) t^2 - 2*t + 1 + By definition, the characteristic polynomial is a monic + degree-zero polynomial in a rank-zero algebra. Note that + Cayley-Hamilton is indeed satisfied since the polynomial + ``1`` evaluates to the identity element of the algebra on + any argument:: + + sage: J = TrivialEJA() + sage: J.characteristic_polynomial_of() + 1 + """ r = self.rank() n = self.dimension() - # The list of coefficient polynomials a_1, a_2, ..., a_n. - a = [ self._charpoly_coeff(i) for i in range(n) ] + # The list of coefficient polynomials a_0, a_1, a_2, ..., a_(r-1). + a = self._charpoly_coefficients() # We go to a bit of trouble here to reorder the # indeterminates, so that it's easier to evaluate the # characteristic polynomial at x's coordinates and get back # something in terms of t, which is what we want. - R = a[0].parent() S = PolynomialRing(self.base_ring(),'t') t = S.gen(0) - S = PolynomialRing(S, R.variable_names()) - t = S(t) - - # Note: all entries past the rth should be zero. The - # coefficient of the highest power (x^r) is 1, but it doesn't - # appear in the solution vector which contains coefficients - # for the other powers (to make them sum to x^r). - if (r < n): - a[r] = 1 # corresponds to x^r - else: - # When the rank is equal to the dimension, trying to - # assign a[r] goes out-of-bounds. - a.append(1) # corresponds to x^r + if r > 0: + R = a[0].parent() + S = PolynomialRing(S, R.variable_names()) + t = S(t) + + return (t**r + sum( a[k]*(t**k) for k in range(r) )) + + def coordinate_polynomial_ring(self): + r""" + The multivariate polynomial ring in which this algebra's + :meth:`characteristic_polynomial_of` lives. + + SETUP:: + + sage: from mjo.eja.eja_algebra import (HadamardEJA, + ....: RealSymmetricEJA) + + EXAMPLES:: - return sum( a[k]*(t**k) for k in range(len(a)) ) + sage: J = HadamardEJA(2) + sage: J.coordinate_polynomial_ring() + Multivariate Polynomial Ring in X1, X2... + sage: J = RealSymmetricEJA(3,field=QQ,orthonormalize=False) + sage: J.coordinate_polynomial_ring() + Multivariate Polynomial Ring in X1, X2, X3, X4, X5, X6... + """ + var_names = tuple( "X%d" % z for z in range(1, self.dimension()+1) ) + return PolynomialRing(self.base_ring(), var_names) def inner_product(self, x, y): """ @@ -310,39 +515,149 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra): SETUP:: - sage: from mjo.eja.eja_algebra import random_eja + sage: from mjo.eja.eja_algebra import (random_eja, + ....: HadamardEJA, + ....: BilinearFormEJA) EXAMPLES: - The inner product must satisfy its axiom for this algebra to truly - be a Euclidean Jordan Algebra:: + Our inner product is "associative," which means the following for + a symmetric bilinear form:: sage: set_random_seed() sage: J = random_eja() + sage: x,y,z = J.random_elements(3) + sage: (x*y).inner_product(z) == y.inner_product(x*z) + True + + TESTS: + + Ensure that this is the usual inner product for the algebras + over `R^n`:: + + sage: set_random_seed() + sage: J = HadamardEJA.random_instance() + sage: x,y = J.random_elements(2) + sage: actual = x.inner_product(y) + sage: expected = x.to_vector().inner_product(y.to_vector()) + sage: actual == expected + True + + Ensure that this is one-half of the trace inner-product in a + BilinearFormEJA that isn't just the reals (when ``n`` isn't + one). This is in Faraut and Koranyi, and also my "On the + symmetry..." paper:: + + sage: set_random_seed() + sage: J = BilinearFormEJA.random_instance() + sage: n = J.dimension() sage: x = J.random_element() sage: y = J.random_element() - sage: z = J.random_element() - sage: (x*y).inner_product(z) == y.inner_product(x*z) + sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2) + True + """ + B = self._inner_product_matrix + return (B*x.to_vector()).inner_product(y.to_vector()) + + + def is_trivial(self): + """ + Return whether or not this algebra is trivial. + + A trivial algebra contains only the zero element. + + SETUP:: + + sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA, + ....: TrivialEJA) + + EXAMPLES:: + + sage: J = ComplexHermitianEJA(3) + sage: J.is_trivial() + False + + :: + + sage: J = TrivialEJA() + sage: J.is_trivial() True """ - if (not x in self) or (not y in self): - raise TypeError("arguments must live in this algebra") - return x.trace_inner_product(y) + return self.dimension() == 0 + + + def multiplication_table(self): + """ + Return a visual representation of this algebra's multiplication + table (on basis elements). + SETUP:: + + sage: from mjo.eja.eja_algebra import JordanSpinEJA + + EXAMPLES:: + + sage: J = JordanSpinEJA(4) + sage: J.multiplication_table() + +----++----+----+----+----+ + | * || e0 | e1 | e2 | e3 | + +====++====+====+====+====+ + | e0 || e0 | e1 | e2 | e3 | + +----++----+----+----+----+ + | e1 || e1 | e0 | 0 | 0 | + +----++----+----+----+----+ + | e2 || e2 | 0 | e0 | 0 | + +----++----+----+----+----+ + | e3 || e3 | 0 | 0 | e0 | + +----++----+----+----+----+ - def natural_basis(self): """ - Return a more-natural representation of this algebra's basis. + n = self.dimension() + M = [ [ self.zero() for j in range(n) ] + for i in range(n) ] + for i in range(n): + for j in range(i+1): + M[i][j] = self._multiplication_table[i][j] + M[j][i] = M[i][j] + + for i in range(n): + # Prepend the left "header" column entry Can't do this in + # the loop because it messes up the symmetry. + M[i] = [self.monomial(i)] + M[i] - Every finite-dimensional Euclidean Jordan Algebra is a direct - sum of five simple algebras, four of which comprise Hermitian - matrices. This method returns the original "natural" basis - for our underlying vector space. (Typically, the natural basis - is used to construct the multiplication table in the first place.) + # Prepend the header row. + M = [["*"] + list(self.gens())] + M + return table(M, header_row=True, header_column=True, frame=True) - Note that this will always return a matrix. The standard basis - in `R^n` will be returned as `n`-by-`1` column matrices. + + def matrix_basis(self): + """ + Return an (often more natural) representation of this algebras + basis as an ordered tuple of matrices. + + Every finite-dimensional Euclidean Jordan Algebra is a, up to + Jordan isomorphism, a direct sum of five simple + algebras---four of which comprise Hermitian matrices. And the + last type of algebra can of course be thought of as `n`-by-`1` + column matrices (ambiguusly called column vectors) to avoid + special cases. As a result, matrices (and column vectors) are + a natural representation format for Euclidean Jordan algebra + elements. + + But, when we construct an algebra from a basis of matrices, + those matrix representations are lost in favor of coordinate + vectors *with respect to* that basis. We could eventually + convert back if we tried hard enough, but having the original + representations handy is valuable enough that we simply store + them and return them from this method. + + Why implement this for non-matrix algebras? Avoiding special + cases for the :class:`BilinearFormEJA` pays with simplicity in + its own right. But mainly, we would like to be able to assume + that elements of a :class:`DirectSumEJA` can be displayed + nicely, without having to have special classes for direct sums + one of whose components was a matrix algebra. SETUP:: @@ -353,607 +668,860 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra): sage: J = RealSymmetricEJA(2) sage: J.basis() - Family (e0, e1, e2) - sage: J.natural_basis() + Finite family {0: e0, 1: e1, 2: e2} + sage: J.matrix_basis() ( - [1 0] [0 1] [0 0] - [0 0], [1 0], [0 1] + [1 0] [ 0 0.7071067811865475?] [0 0] + [0 0], [0.7071067811865475? 0], [0 1] ) :: sage: J = JordanSpinEJA(2) sage: J.basis() - Family (e0, e1) - sage: J.natural_basis() + Finite family {0: e0, 1: e1} + sage: J.matrix_basis() ( [1] [0] [0], [1] ) - """ - if self._natural_basis is None: - return tuple( b.vector().column() for b in self.basis() ) + if self._matrix_basis is None: + M = self.matrix_space() + return tuple( M(b.to_vector()) for b in self.basis() ) else: - return self._natural_basis + return self._matrix_basis - def rank(self): + def matrix_space(self): """ - Return the rank of this EJA. + Return the matrix space in which this algebra's elements live, if + we think of them as matrices (including column vectors of the + appropriate size). + + Generally this will be an `n`-by-`1` column-vector space, + except when the algebra is trivial. There it's `n`-by-`n` + (where `n` is zero), to ensure that two elements of the matrix + space (empty matrices) can be multiplied. + + Matrix algebras override this with something more useful. + """ + if self.is_trivial(): + return MatrixSpace(self.base_ring(), 0) + elif self._matrix_basis is None or len(self._matrix_basis) == 0: + return MatrixSpace(self.base_ring(), self.dimension(), 1) + else: + return self._matrix_basis[0].matrix_space() - ALGORITHM: - The author knows of no algorithm to compute the rank of an EJA - where only the multiplication table is known. In lieu of one, we - require the rank to be specified when the algebra is created, - and simply pass along that number here. + @cached_method + def one(self): + """ + Return the unit element of this algebra. SETUP:: - sage: from mjo.eja.eja_algebra import (JordanSpinEJA, - ....: RealSymmetricEJA, - ....: ComplexHermitianEJA, - ....: QuaternionHermitianEJA, + sage: from mjo.eja.eja_algebra import (HadamardEJA, ....: random_eja) - EXAMPLES: - - The rank of the Jordan spin algebra is always two:: + EXAMPLES:: - sage: JordanSpinEJA(2).rank() - 2 - sage: JordanSpinEJA(3).rank() - 2 - sage: JordanSpinEJA(4).rank() - 2 + sage: J = HadamardEJA(5) + sage: J.one() + e0 + e1 + e2 + e3 + e4 - The rank of the `n`-by-`n` Hermitian real, complex, or - quaternion matrices is `n`:: + TESTS: - sage: RealSymmetricEJA(2).rank() - 2 - sage: ComplexHermitianEJA(2).rank() - 2 - sage: QuaternionHermitianEJA(2).rank() - 2 - sage: RealSymmetricEJA(5).rank() - 5 - sage: ComplexHermitianEJA(5).rank() - 5 - sage: QuaternionHermitianEJA(5).rank() - 5 + The identity element acts like the identity:: - TESTS: + sage: set_random_seed() + sage: J = random_eja() + sage: x = J.random_element() + sage: J.one()*x == x and x*J.one() == x + True - Ensure that every EJA that we know how to construct has a - positive integer rank:: + The matrix of the unit element's operator is the identity:: sage: set_random_seed() - sage: r = random_eja().rank() - sage: r in ZZ and r > 0 + sage: J = random_eja() + sage: actual = J.one().operator().matrix() + sage: expected = matrix.identity(J.base_ring(), J.dimension()) + sage: actual == expected True - """ - return self._rank + Ensure that the cached unit element (often precomputed by + hand) agrees with the computed one:: + sage: set_random_seed() + sage: J = random_eja() + sage: cached = J.one() + sage: J.one.clear_cache() + sage: J.one() == cached + True - def vector_space(self): """ - Return the vector space that underlies this algebra. + # We can brute-force compute the matrices of the operators + # that correspond to the basis elements of this algebra. + # If some linear combination of those basis elements is the + # algebra identity, then the same linear combination of + # their matrices has to be the identity matrix. + # + # Of course, matrices aren't vectors in sage, so we have to + # appeal to the "long vectors" isometry. + oper_vecs = [ _mat2vec(g.operator().matrix()) for g in self.gens() ] - SETUP:: + # Now we use basic linear algebra to find the coefficients, + # of the matrices-as-vectors-linear-combination, which should + # work for the original algebra basis too. + A = matrix(self.base_ring(), oper_vecs) - sage: from mjo.eja.eja_algebra import RealSymmetricEJA + # We used the isometry on the left-hand side already, but we + # still need to do it for the right-hand side. Recall that we + # wanted something that summed to the identity matrix. + b = _mat2vec( matrix.identity(self.base_ring(), self.dimension()) ) - EXAMPLES:: + # Now if there's an identity element in the algebra, this + # should work. We solve on the left to avoid having to + # transpose the matrix "A". + return self.from_vector(A.solve_left(b)) - sage: J = RealSymmetricEJA(2) - sage: J.vector_space() - Vector space of dimension 3 over Rational Field + def peirce_decomposition(self, c): """ - return self.zero().vector().parent().ambient_vector_space() + The Peirce decomposition of this algebra relative to the + idempotent ``c``. + In the future, this can be extended to a complete system of + orthogonal idempotents. - Element = FiniteDimensionalEuclideanJordanAlgebraElement + INPUT: + - ``c`` -- an idempotent of this algebra. -class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra): - """ - Return the Euclidean Jordan Algebra corresponding to the set - `R^n` under the Hadamard product. + OUTPUT: - Note: this is nothing more than the Cartesian product of ``n`` - copies of the spin algebra. Once Cartesian product algebras - are implemented, this can go. + A triple (J0, J5, J1) containing two subalgebras and one subspace + of this algebra, - SETUP:: + - ``J0`` -- the algebra on the eigenspace of ``c.operator()`` + corresponding to the eigenvalue zero. - sage: from mjo.eja.eja_algebra import RealCartesianProductEJA + - ``J5`` -- the eigenspace (NOT a subalgebra) of ``c.operator()`` + corresponding to the eigenvalue one-half. - EXAMPLES: + - ``J1`` -- the algebra on the eigenspace of ``c.operator()`` + corresponding to the eigenvalue one. - This multiplication table can be verified by hand:: + These are the only possible eigenspaces for that operator, and this + algebra is a direct sum of them. The spaces ``J0`` and ``J1`` are + orthogonal, and are subalgebras of this algebra with the appropriate + restrictions. - sage: J = RealCartesianProductEJA(3) - sage: e0,e1,e2 = J.gens() - sage: e0*e0 - e0 - sage: e0*e1 - 0 - sage: e0*e2 - 0 - sage: e1*e1 - e1 - sage: e1*e2 - 0 - sage: e2*e2 - e2 + SETUP:: - """ - @staticmethod - def __classcall_private__(cls, n, field=QQ): - # The FiniteDimensionalAlgebra constructor takes a list of - # matrices, the ith representing right multiplication by the ith - # basis element in the vector space. So if e_1 = (1,0,0), then - # right (Hadamard) multiplication of x by e_1 picks out the first - # component of x; and likewise for the ith basis element e_i. - Qs = [ matrix(field, n, n, lambda k,j: 1*(k == j == i)) - for i in xrange(n) ] - - fdeja = super(RealCartesianProductEJA, cls) - return fdeja.__classcall_private__(cls, field, Qs, rank=n) + sage: from mjo.eja.eja_algebra import random_eja, RealSymmetricEJA - def inner_product(self, x, y): - return _usual_ip(x,y) + EXAMPLES: + The canonical example comes from the symmetric matrices, which + decompose into diagonal and off-diagonal parts:: + + sage: J = RealSymmetricEJA(3) + sage: C = matrix(QQ, [ [1,0,0], + ....: [0,1,0], + ....: [0,0,0] ]) + sage: c = J(C) + sage: J0,J5,J1 = J.peirce_decomposition(c) + sage: J0 + Euclidean Jordan algebra of dimension 1... + sage: J5 + Vector space of degree 6 and dimension 2... + sage: J1 + Euclidean Jordan algebra of dimension 3... + sage: J0.one().to_matrix() + [0 0 0] + [0 0 0] + [0 0 1] + sage: orig_df = AA.options.display_format + sage: AA.options.display_format = 'radical' + sage: J.from_vector(J5.basis()[0]).to_matrix() + [ 0 0 1/2*sqrt(2)] + [ 0 0 0] + [1/2*sqrt(2) 0 0] + sage: J.from_vector(J5.basis()[1]).to_matrix() + [ 0 0 0] + [ 0 0 1/2*sqrt(2)] + [ 0 1/2*sqrt(2) 0] + sage: AA.options.display_format = orig_df + sage: J1.one().to_matrix() + [1 0 0] + [0 1 0] + [0 0 0] -def random_eja(): - """ - Return a "random" finite-dimensional Euclidean Jordan Algebra. + TESTS: - ALGORITHM: + Every algebra decomposes trivially with respect to its identity + element:: - For now, we choose a random natural number ``n`` (greater than zero) - and then give you back one of the following: + sage: set_random_seed() + sage: J = random_eja() + sage: J0,J5,J1 = J.peirce_decomposition(J.one()) + sage: J0.dimension() == 0 and J5.dimension() == 0 + True + sage: J1.superalgebra() == J and J1.dimension() == J.dimension() + True - * The cartesian product of the rational numbers ``n`` times; this is - ``QQ^n`` with the Hadamard product. + The decomposition is into eigenspaces, and its components are + therefore necessarily orthogonal. Moreover, the identity + elements in the two subalgebras are the projections onto their + respective subspaces of the superalgebra's identity element:: - * The Jordan spin algebra on ``QQ^n``. + sage: set_random_seed() + sage: J = random_eja() + sage: x = J.random_element() + sage: if not J.is_trivial(): + ....: while x.is_nilpotent(): + ....: x = J.random_element() + sage: c = x.subalgebra_idempotent() + sage: J0,J5,J1 = J.peirce_decomposition(c) + sage: ipsum = 0 + sage: for (w,y,z) in zip(J0.basis(), J5.basis(), J1.basis()): + ....: w = w.superalgebra_element() + ....: y = J.from_vector(y) + ....: z = z.superalgebra_element() + ....: ipsum += w.inner_product(y).abs() + ....: ipsum += w.inner_product(z).abs() + ....: ipsum += y.inner_product(z).abs() + sage: ipsum + 0 + sage: J1(c) == J1.one() + True + sage: J0(J.one() - c) == J0.one() + True - * The ``n``-by-``n`` rational symmetric matrices with the symmetric - product. + """ + if not c.is_idempotent(): + raise ValueError("element is not idempotent: %s" % c) + + from mjo.eja.eja_subalgebra import FiniteDimensionalEuclideanJordanSubalgebra + + # Default these to what they should be if they turn out to be + # trivial, because eigenspaces_left() won't return eigenvalues + # corresponding to trivial spaces (e.g. it returns only the + # eigenspace corresponding to lambda=1 if you take the + # decomposition relative to the identity element). + trivial = FiniteDimensionalEuclideanJordanSubalgebra(self, ()) + J0 = trivial # eigenvalue zero + J5 = VectorSpace(self.base_ring(), 0) # eigenvalue one-half + J1 = trivial # eigenvalue one + + for (eigval, eigspace) in c.operator().matrix().right_eigenspaces(): + if eigval == ~(self.base_ring()(2)): + J5 = eigspace + else: + gens = tuple( self.from_vector(b) for b in eigspace.basis() ) + subalg = FiniteDimensionalEuclideanJordanSubalgebra(self, + gens, + check_axioms=False) + if eigval == 0: + J0 = subalg + elif eigval == 1: + J1 = subalg + else: + raise ValueError("unexpected eigenvalue: %s" % eigval) - * The ``n``-by-``n`` complex-rational Hermitian matrices embedded - in the space of ``2n``-by-``2n`` real symmetric matrices. + return (J0, J5, J1) - * The ``n``-by-``n`` quaternion-rational Hermitian matrices embedded - in the space of ``4n``-by-``4n`` real symmetric matrices. - Later this might be extended to return Cartesian products of the - EJAs above. + def random_element(self, thorough=False): + r""" + Return a random element of this algebra. - SETUP:: + Our algebra superclass method only returns a linear + combination of at most two basis elements. We instead + want the vector space "random element" method that + returns a more diverse selection. - sage: from mjo.eja.eja_algebra import random_eja + INPUT: - TESTS:: + - ``thorough`` -- (boolean; default False) whether or not we + should generate irrational coefficients for the random + element when our base ring is irrational; this slows the + algebra operations to a crawl, but any truly random method + should include them - sage: random_eja() - Euclidean Jordan algebra of degree... + """ + # For a general base ring... maybe we can trust this to do the + # right thing? Unlikely, but. + V = self.vector_space() + v = V.random_element() + + if self.base_ring() is AA: + # The "random element" method of the algebraic reals is + # stupid at the moment, and only returns integers between + # -2 and 2, inclusive: + # + # https://trac.sagemath.org/ticket/30875 + # + # Instead, we implement our own "random vector" method, + # and then coerce that into the algebra. We use the vector + # space degree here instead of the dimension because a + # subalgebra could (for example) be spanned by only two + # vectors, each with five coordinates. We need to + # generate all five coordinates. + if thorough: + v *= QQbar.random_element().real() + else: + v *= QQ.random_element() - """ + return self.from_vector(V.coordinate_vector(v)) - # The max_n component lets us choose different upper bounds on the - # value "n" that gets passed to the constructor. This is needed - # because e.g. R^{10} is reasonable to test, while the Hermitian - # 10-by-10 quaternion matrices are not. - (constructor, max_n) = choice([(RealCartesianProductEJA, 6), - (JordanSpinEJA, 6), - (RealSymmetricEJA, 5), - (ComplexHermitianEJA, 4), - (QuaternionHermitianEJA, 3)]) - n = ZZ.random_element(1, max_n) - return constructor(n, field=QQ) + def random_elements(self, count, thorough=False): + """ + Return ``count`` random elements as a tuple. + INPUT: + - ``thorough`` -- (boolean; default False) whether or not we + should generate irrational coefficients for the random + elements when our base ring is irrational; this slows the + algebra operations to a crawl, but any truly random method + should include them -def _real_symmetric_basis(n, field=QQ): - """ - Return a basis for the space of real symmetric n-by-n matrices. - """ - # The basis of symmetric matrices, as matrices, in their R^(n-by-n) - # coordinates. - S = [] - for i in xrange(n): - for j in xrange(i+1): - Eij = matrix(field, n, lambda k,l: k==i and l==j) - if i == j: - Sij = Eij - else: - # Beware, orthogonal but not normalized! - Sij = Eij + Eij.transpose() - S.append(Sij) - return tuple(S) + SETUP:: + sage: from mjo.eja.eja_algebra import JordanSpinEJA -def _complex_hermitian_basis(n, field=QQ): - """ - Returns a basis for the space of complex Hermitian n-by-n matrices. + EXAMPLES:: - SETUP:: + sage: J = JordanSpinEJA(3) + sage: x,y,z = J.random_elements(3) + sage: all( [ x in J, y in J, z in J ]) + True + sage: len( J.random_elements(10) ) == 10 + True - sage: from mjo.eja.eja_algebra import _complex_hermitian_basis + """ + return tuple( self.random_element(thorough) + for idx in range(count) ) - TESTS:: - sage: set_random_seed() - sage: n = ZZ.random_element(1,5) - sage: all( M.is_symmetric() for M in _complex_hermitian_basis(n) ) - True + @cached_method + def _charpoly_coefficients(self): + r""" + The `r` polynomial coefficients of the "characteristic polynomial + of" function. + """ + n = self.dimension() + R = self.coordinate_polynomial_ring() + vars = R.gens() + F = R.fraction_field() + + def L_x_i_j(i,j): + # From a result in my book, these are the entries of the + # basis representation of L_x. + return sum( vars[k]*self.monomial(k).operator().matrix()[i,j] + for k in range(n) ) + + L_x = matrix(F, n, n, L_x_i_j) + + r = None + if self.rank.is_in_cache(): + r = self.rank() + # There's no need to pad the system with redundant + # columns if we *know* they'll be redundant. + n = r + + # Compute an extra power in case the rank is equal to + # the dimension (otherwise, we would stop at x^(r-1)). + x_powers = [ (L_x**k)*self.one().to_vector() + for k in range(n+1) ] + A = matrix.column(F, x_powers[:n]) + AE = A.extended_echelon_form() + E = AE[:,n:] + A_rref = AE[:,:n] + if r is None: + r = A_rref.rank() + b = x_powers[r] + + # The theory says that only the first "r" coefficients are + # nonzero, and they actually live in the original polynomial + # ring and not the fraction field. We negate them because + # in the actual characteristic polynomial, they get moved + # to the other side where x^r lives. + return -A_rref.solve_right(E*b).change_ring(R)[:r] - """ - F = QuadraticField(-1, 'I') - I = F.gen() - - # This is like the symmetric case, but we need to be careful: - # - # * We want conjugate-symmetry, not just symmetry. - # * The diagonal will (as a result) be real. - # - S = [] - for i in xrange(n): - for j in xrange(i+1): - Eij = matrix(field, n, lambda k,l: k==i and l==j) - if i == j: - Sij = _embed_complex_matrix(Eij) - S.append(Sij) - else: - # Beware, orthogonal but not normalized! The second one - # has a minus because it's conjugated. - Sij_real = _embed_complex_matrix(Eij + Eij.transpose()) - S.append(Sij_real) - Sij_imag = _embed_complex_matrix(I*Eij - I*Eij.transpose()) - S.append(Sij_imag) - return tuple(S) + @cached_method + def rank(self): + r""" + Return the rank of this EJA. + This is a cached method because we know the rank a priori for + all of the algebras we can construct. Thus we can avoid the + expensive ``_charpoly_coefficients()`` call unless we truly + need to compute the whole characteristic polynomial. -def _quaternion_hermitian_basis(n, field=QQ): - """ - Returns a basis for the space of quaternion Hermitian n-by-n matrices. + SETUP:: - SETUP:: + sage: from mjo.eja.eja_algebra import (HadamardEJA, + ....: JordanSpinEJA, + ....: RealSymmetricEJA, + ....: ComplexHermitianEJA, + ....: QuaternionHermitianEJA, + ....: random_eja) - sage: from mjo.eja.eja_algebra import _quaternion_hermitian_basis + EXAMPLES: - TESTS:: + The rank of the Jordan spin algebra is always two:: - sage: set_random_seed() - sage: n = ZZ.random_element(1,5) - sage: all( M.is_symmetric() for M in _quaternion_hermitian_basis(n) ) - True + sage: JordanSpinEJA(2).rank() + 2 + sage: JordanSpinEJA(3).rank() + 2 + sage: JordanSpinEJA(4).rank() + 2 - """ - Q = QuaternionAlgebra(QQ,-1,-1) - I,J,K = Q.gens() - - # This is like the symmetric case, but we need to be careful: - # - # * We want conjugate-symmetry, not just symmetry. - # * The diagonal will (as a result) be real. - # - S = [] - for i in xrange(n): - for j in xrange(i+1): - Eij = matrix(Q, n, lambda k,l: k==i and l==j) - if i == j: - Sij = _embed_quaternion_matrix(Eij) - S.append(Sij) - else: - # Beware, orthogonal but not normalized! The second, - # third, and fourth ones have a minus because they're - # conjugated. - Sij_real = _embed_quaternion_matrix(Eij + Eij.transpose()) - S.append(Sij_real) - Sij_I = _embed_quaternion_matrix(I*Eij - I*Eij.transpose()) - S.append(Sij_I) - Sij_J = _embed_quaternion_matrix(J*Eij - J*Eij.transpose()) - S.append(Sij_J) - Sij_K = _embed_quaternion_matrix(K*Eij - K*Eij.transpose()) - S.append(Sij_K) - return tuple(S) - - - -def _multiplication_table_from_matrix_basis(basis): - """ - At least three of the five simple Euclidean Jordan algebras have the - symmetric multiplication (A,B) |-> (AB + BA)/2, where the - multiplication on the right is matrix multiplication. Given a basis - for the underlying matrix space, this function returns a - multiplication table (obtained by looping through the basis - elements) for an algebra of those matrices. A reordered copy - of the basis is also returned to work around the fact that - the ``span()`` in this function will change the order of the basis - from what we think it is, to... something else. - """ - # In S^2, for example, we nominally have four coordinates even - # though the space is of dimension three only. The vector space V - # is supposed to hold the entire long vector, and the subspace W - # of V will be spanned by the vectors that arise from symmetric - # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3. - field = basis[0].base_ring() - dimension = basis[0].nrows() - - V = VectorSpace(field, dimension**2) - W = V.span( _mat2vec(s) for s in basis ) - - # Taking the span above reorders our basis (thanks, jerk!) so we - # need to put our "matrix basis" in the same order as the - # (reordered) vector basis. - S = tuple( _vec2mat(b) for b in W.basis() ) - - Qs = [] - for s in S: - # Brute force the multiplication-by-s matrix by looping - # through all elements of the basis and doing the computation - # to find out what the corresponding row should be. BEWARE: - # these multiplication tables won't be symmetric! It therefore - # becomes REALLY IMPORTANT that the underlying algebra - # constructor uses ROW vectors and not COLUMN vectors. That's - # why we're computing rows here and not columns. - Q_rows = [] - for t in S: - this_row = _mat2vec((s*t + t*s)/2) - Q_rows.append(W.coordinates(this_row)) - Q = matrix(field, W.dimension(), Q_rows) - Qs.append(Q) - - return (Qs, S) - - -def _embed_complex_matrix(M): - """ - Embed the n-by-n complex matrix ``M`` into the space of real - matrices of size 2n-by-2n via the map the sends each entry `z = a + - bi` to the block matrix ``[[a,b],[-b,a]]``. + The rank of the `n`-by-`n` Hermitian real, complex, or + quaternion matrices is `n`:: - SETUP:: + sage: RealSymmetricEJA(4).rank() + 4 + sage: ComplexHermitianEJA(3).rank() + 3 + sage: QuaternionHermitianEJA(2).rank() + 2 - sage: from mjo.eja.eja_algebra import _embed_complex_matrix + TESTS: - EXAMPLES:: + Ensure that every EJA that we know how to construct has a + positive integer rank, unless the algebra is trivial in + which case its rank will be zero:: - sage: F = QuadraticField(-1,'i') - sage: x1 = F(4 - 2*i) - sage: x2 = F(1 + 2*i) - sage: x3 = F(-i) - sage: x4 = F(6) - sage: M = matrix(F,2,[[x1,x2],[x3,x4]]) - sage: _embed_complex_matrix(M) - [ 4 -2| 1 2] - [ 2 4|-2 1] - [-----+-----] - [ 0 -1| 6 0] - [ 1 0| 0 6] + sage: set_random_seed() + sage: J = random_eja() + sage: r = J.rank() + sage: r in ZZ + True + sage: r > 0 or (r == 0 and J.is_trivial()) + True - TESTS: + Ensure that computing the rank actually works, since the ranks + of all simple algebras are known and will be cached by default:: - Embedding is a homomorphism (isomorphism, in fact):: + sage: set_random_seed() # long time + sage: J = random_eja() # long time + sage: caches = J.rank() # long time + sage: J.rank.clear_cache() # long time + sage: J.rank() == cached # long time + True - sage: set_random_seed() - sage: n = ZZ.random_element(5) - sage: F = QuadraticField(-1, 'i') - sage: X = random_matrix(F, n) - sage: Y = random_matrix(F, n) - sage: actual = _embed_complex_matrix(X) * _embed_complex_matrix(Y) - sage: expected = _embed_complex_matrix(X*Y) - sage: actual == expected - True + """ + return len(self._charpoly_coefficients()) - """ - n = M.nrows() - if M.ncols() != n: - raise ValueError("the matrix 'M' must be square") - field = M.base_ring() - blocks = [] - for z in M.list(): - a = z.real() - b = z.imag() - blocks.append(matrix(field, 2, [[a,b],[-b,a]])) - # We can drop the imaginaries here. - return matrix.block(field.base_ring(), n, blocks) + def vector_space(self): + """ + Return the vector space that underlies this algebra. + + SETUP:: + sage: from mjo.eja.eja_algebra import RealSymmetricEJA + + EXAMPLES:: + + sage: J = RealSymmetricEJA(2) + sage: J.vector_space() + Vector space of dimension 3 over... + + """ + return self.zero().to_vector().parent().ambient_vector_space() -def _unembed_complex_matrix(M): - """ - The inverse of _embed_complex_matrix(). + + Element = FiniteDimensionalEuclideanJordanAlgebraElement + +class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra): + r""" + New class for algebras whose supplied basis elements have all rational entries. SETUP:: - sage: from mjo.eja.eja_algebra import (_embed_complex_matrix, - ....: _unembed_complex_matrix) + sage: from mjo.eja.eja_algebra import BilinearFormEJA - EXAMPLES:: + EXAMPLES: + + The supplied basis is orthonormalized by default:: + + sage: B = matrix(QQ, [[1, 0, 0], [0, 25, -32], [0, -32, 41]]) + sage: J = BilinearFormEJA(B) + sage: J.matrix_basis() + ( + [1] [ 0] [ 0] + [0] [1/5] [32/5] + [0], [ 0], [ 5] + ) + + """ + def __init__(self, + basis, + jordan_product, + inner_product, + field=AA, + orthonormalize=True, + prefix='e', + category=None, + check_field=True, + check_axioms=True): + + if check_field: + # Abuse the check_field parameter to check that the entries of + # out basis (in ambient coordinates) are in the field QQ. + if not all( all(b_i in QQ for b_i in b.list()) for b in basis ): + raise TypeError("basis not rational") + + # Temporary(?) hack to ensure that the matrix and vector bases + # are over the same ring. + basis = tuple( b.change_ring(field) for b in basis ) + + n = len(basis) + vector_basis = basis + + from sage.structure.element import is_Matrix + basis_is_matrices = False + + degree = 0 + if n > 0: + if is_Matrix(basis[0]): + basis_is_matrices = True + from mjo.eja.eja_utils import _vec2mat + vector_basis = tuple( map(_mat2vec,basis) ) + degree = basis[0].nrows()**2 + else: + degree = basis[0].degree() + + V = VectorSpace(field, degree) + + # If we were asked to orthonormalize, and if the orthonormal + # basis is different from the given one, then we also want to + # compute multiplication and inner-product tables for the + # deorthonormalized basis. These can be used later to + # construct a deorthonormalized copy of this algebra over QQ + # in which several operations are much faster. + self._rational_algebra = None + + if orthonormalize: + if self.base_ring() is not QQ: + # There's no point in constructing the extra algebra if this + # one is already rational. If the original basis is rational + # but normalization would make it irrational, then this whole + # constructor will just fail anyway as it tries to stick an + # irrational number into a rational algebra. + # + # Note: the same Jordan and inner-products work here, + # because they are necessarily defined with respect to + # ambient coordinates and not any particular basis. + self._rational_algebra = RationalBasisEuclideanJordanAlgebra( + basis, + jordan_product, + inner_product, + field=QQ, + orthonormalize=False, + prefix=prefix, + category=category, + check_field=False, + check_axioms=False) + + # Compute the deorthonormalized tables before we orthonormalize + # the given basis. The "check" parameter here guarantees that + # the basis is linearly-independent. + W = V.span_of_basis( vector_basis, check=check_axioms) + + # Note: the Jordan and inner-products are defined in terms + # of the ambient basis. It's important that their arguments + # are in ambient coordinates as well. + for i in range(n): + for j in range(i+1): + # given basis w.r.t. ambient coords + q_i = vector_basis[i] + q_j = vector_basis[j] + + if basis_is_matrices: + q_i = _vec2mat(q_i) + q_j = _vec2mat(q_j) + + elt = jordan_product(q_i, q_j) + ip = inner_product(q_i, q_j) + + if basis_is_matrices: + # do another mat2vec because the multiplication + # table is in terms of vectors + elt = _mat2vec(elt) + + # We overwrite the name "vector_basis" in a second, but never modify it + # in place, to this effectively makes a copy of it. + deortho_vector_basis = vector_basis + self._deortho_matrix = None + + if orthonormalize: + from mjo.eja.eja_utils import gram_schmidt + if basis_is_matrices: + vector_ip = lambda x,y: inner_product(_vec2mat(x), _vec2mat(y)) + vector_basis = gram_schmidt(vector_basis, vector_ip) + else: + vector_basis = gram_schmidt(vector_basis, inner_product) + + # Normalize the "matrix" basis, too! + basis = vector_basis + + if basis_is_matrices: + basis = tuple( map(_vec2mat,basis) ) + + W = V.span_of_basis( vector_basis, check=check_axioms) + + # Now "W" is the vector space of our algebra coordinates. The + # variables "X1", "X2",... refer to the entries of vectors in + # W. Thus to convert back and forth between the orthonormal + # coordinates and the given ones, we need to stick the original + # basis in W. + U = V.span_of_basis( deortho_vector_basis, check=check_axioms) + self._deortho_matrix = matrix( U.coordinate_vector(q) + for q in vector_basis ) + + # If the superclass constructor is going to verify the + # symmetry of this table, it has better at least be + # square... + if check_axioms: + mult_table = [ [0 for j in range(n)] for i in range(n) ] + ip_table = [ [0 for j in range(n)] for i in range(n) ] + else: + mult_table = [ [0 for j in range(i+1)] for i in range(n) ] + ip_table = [ [0 for j in range(i+1)] for i in range(n) ] + + # Note: the Jordan and inner-products are defined in terms + # of the ambient basis. It's important that their arguments + # are in ambient coordinates as well. + for i in range(n): + for j in range(i+1): + # ortho basis w.r.t. ambient coords + q_i = vector_basis[i] + q_j = vector_basis[j] + + if basis_is_matrices: + q_i = _vec2mat(q_i) + q_j = _vec2mat(q_j) + + elt = jordan_product(q_i, q_j) + ip = inner_product(q_i, q_j) + + if basis_is_matrices: + # do another mat2vec because the multiplication + # table is in terms of vectors + elt = _mat2vec(elt) + + elt = W.coordinate_vector(elt) + mult_table[i][j] = elt + ip_table[i][j] = ip + if check_axioms: + # The tables are square if we're verifying that they + # are commutative. + mult_table[j][i] = elt + ip_table[j][i] = ip + + if basis_is_matrices: + for m in basis: + m.set_immutable() + else: + basis = tuple( x.column() for x in basis ) - sage: A = matrix(QQ,[ [ 1, 2, 3, 4], - ....: [-2, 1, -4, 3], - ....: [ 9, 10, 11, 12], - ....: [-10, 9, -12, 11] ]) - sage: _unembed_complex_matrix(A) - [ 2*i + 1 4*i + 3] - [ 10*i + 9 12*i + 11] + super().__init__(field, + mult_table, + ip_table, + prefix, + category, + basis, # matrix basis + check_field, + check_axioms) + + @cached_method + def _charpoly_coefficients(self): + r""" + SETUP:: + + sage: from mjo.eja.eja_algebra import (BilinearFormEJA, + ....: JordanSpinEJA) + + EXAMPLES: + + The base ring of the resulting polynomial coefficients is what + it should be, and not the rationals (unless the algebra was + already over the rationals):: + + sage: J = JordanSpinEJA(3) + sage: J._charpoly_coefficients() + (X1^2 - X2^2 - X3^2, -2*X1) + sage: a0 = J._charpoly_coefficients()[0] + sage: J.base_ring() + Algebraic Real Field + sage: a0.base_ring() + Algebraic Real Field + + """ + if self.base_ring() is QQ or self._rational_algebra is None: + # There's no need to construct *another* algebra over the + # rationals if this one is already over the + # rationals. Likewise, if we never orthonormalized our + # basis, we might as well just use the given one. + superclass = super(RationalBasisEuclideanJordanAlgebra, self) + return superclass._charpoly_coefficients() + + # Do the computation over the rationals. The answer will be + # the same, because all we've done is a change of basis. + # Then, change back from QQ to our real base ring + a = ( a_i.change_ring(self.base_ring()) + for a_i in self._rational_algebra._charpoly_coefficients() ) + + # Now convert the coordinate variables back to the + # deorthonormalized ones. + R = self.coordinate_polynomial_ring() + from sage.modules.free_module_element import vector + X = vector(R, R.gens()) + BX = self._deortho_matrix*X + + subs_dict = { X[i]: BX[i] for i in range(len(X)) } + return tuple( a_i.subs(subs_dict) for a_i in a ) + +class ConcreteEuclideanJordanAlgebra(RationalBasisEuclideanJordanAlgebra): + r""" + A class for the Euclidean Jordan algebras that we know by name. + + These are the Jordan algebras whose basis, multiplication table, + rank, and so on are known a priori. More to the point, they are + the Euclidean Jordan algebras for which we are able to conjure up + a "random instance." + + SETUP:: + + sage: from mjo.eja.eja_algebra import ConcreteEuclideanJordanAlgebra TESTS: - Unembedding is the inverse of embedding:: + Our basis is normalized with respect to the algebra's inner + product, unless we specify otherwise:: sage: set_random_seed() - sage: F = QuadraticField(-1, 'i') - sage: M = random_matrix(F, 3) - sage: _unembed_complex_matrix(_embed_complex_matrix(M)) == M + sage: J = ConcreteEuclideanJordanAlgebra.random_instance() + sage: all( b.norm() == 1 for b in J.gens() ) True + Since our basis is orthonormal with respect to the algebra's inner + product, and since we know that this algebra is an EJA, any + left-multiplication operator's matrix will be symmetric because + natural->EJA basis representation is an isometry and within the + EJA the operator is self-adjoint by the Jordan axiom:: + + sage: set_random_seed() + sage: J = ConcreteEuclideanJordanAlgebra.random_instance() + sage: x = J.random_element() + sage: x.operator().is_self_adjoint() + True """ - n = ZZ(M.nrows()) - if M.ncols() != n: - raise ValueError("the matrix 'M' must be square") - if not n.mod(2).is_zero(): - raise ValueError("the matrix 'M' must be a complex embedding") - - F = QuadraticField(-1, 'i') - i = F.gen() - - # Go top-left to bottom-right (reading order), converting every - # 2-by-2 block we see to a single complex element. - elements = [] - for k in xrange(n/2): - for j in xrange(n/2): - submat = M[2*k:2*k+2,2*j:2*j+2] - if submat[0,0] != submat[1,1]: - raise ValueError('bad on-diagonal submatrix') - if submat[0,1] != -submat[1,0]: - raise ValueError('bad off-diagonal submatrix') - z = submat[0,0] + submat[0,1]*i - elements.append(z) - - return matrix(F, n/2, elements) - - -def _embed_quaternion_matrix(M): - """ - Embed the n-by-n quaternion matrix ``M`` into the space of real - matrices of size 4n-by-4n by first sending each quaternion entry - `z = a + bi + cj + dk` to the block-complex matrix - ``[[a + bi, c+di],[-c + di, a-bi]]`, and then embedding those into - a real matrix. - SETUP:: + @staticmethod + def _max_random_instance_size(): + """ + Return an integer "size" that is an upper bound on the size of + this algebra when it is used in a random test + case. Unfortunately, the term "size" is ambiguous -- when + dealing with `R^n` under either the Hadamard or Jordan spin + product, the "size" refers to the dimension `n`. When dealing + with a matrix algebra (real symmetric or complex/quaternion + Hermitian), it refers to the size of the matrix, which is far + less than the dimension of the underlying vector space. + + This method must be implemented in each subclass. + """ + raise NotImplementedError - sage: from mjo.eja.eja_algebra import _embed_quaternion_matrix + @classmethod + def random_instance(cls, *args, **kwargs): + """ + Return a random instance of this type of algebra. - EXAMPLES:: + This method should be implemented in each subclass. + """ + from sage.misc.prandom import choice + eja_class = choice(cls.__subclasses__()) - sage: Q = QuaternionAlgebra(QQ,-1,-1) - sage: i,j,k = Q.gens() - sage: x = 1 + 2*i + 3*j + 4*k - sage: M = matrix(Q, 1, [[x]]) - sage: _embed_quaternion_matrix(M) - [ 1 2 3 4] - [-2 1 -4 3] - [-3 4 1 -2] - [-4 -3 2 1] + # These all bubble up to the RationalBasisEuclideanJordanAlgebra + # superclass constructor, so any (kw)args valid there are also + # valid here. + return eja_class.random_instance(*args, **kwargs) - Embedding is a homomorphism (isomorphism, in fact):: - sage: set_random_seed() - sage: n = ZZ.random_element(5) - sage: Q = QuaternionAlgebra(QQ,-1,-1) - sage: X = random_matrix(Q, n) - sage: Y = random_matrix(Q, n) - sage: actual = _embed_quaternion_matrix(X)*_embed_quaternion_matrix(Y) - sage: expected = _embed_quaternion_matrix(X*Y) - sage: actual == expected - True +class MatrixEuclideanJordanAlgebra: + @staticmethod + def real_embed(M): + """ + Embed the matrix ``M`` into a space of real matrices. - """ - quaternions = M.base_ring() - n = M.nrows() - if M.ncols() != n: - raise ValueError("the matrix 'M' must be square") - - F = QuadraticField(-1, 'i') - i = F.gen() - - blocks = [] - for z in M.list(): - t = z.coefficient_tuple() - a = t[0] - b = t[1] - c = t[2] - d = t[3] - cplx_matrix = matrix(F, 2, [[ a + b*i, c + d*i], - [-c + d*i, a - b*i]]) - blocks.append(_embed_complex_matrix(cplx_matrix)) - - # We should have real entries by now, so use the realest field - # we've got for the return value. - return matrix.block(quaternions.base_ring(), n, blocks) - - -def _unembed_quaternion_matrix(M): - """ - The inverse of _embed_quaternion_matrix(). + The matrix ``M`` can have entries in any field at the moment: + the real numbers, complex numbers, or quaternions. And although + they are not a field, we can probably support octonions at some + point, too. This function returns a real matrix that "acts like" + the original with respect to matrix multiplication; i.e. - SETUP:: + real_embed(M*N) = real_embed(M)*real_embed(N) - sage: from mjo.eja.eja_algebra import (_embed_quaternion_matrix, - ....: _unembed_quaternion_matrix) + """ + raise NotImplementedError - EXAMPLES:: - sage: M = matrix(QQ, [[ 1, 2, 3, 4], - ....: [-2, 1, -4, 3], - ....: [-3, 4, 1, -2], - ....: [-4, -3, 2, 1]]) - sage: _unembed_quaternion_matrix(M) - [1 + 2*i + 3*j + 4*k] + @staticmethod + def real_unembed(M): + """ + The inverse of :meth:`real_embed`. + """ + raise NotImplementedError - TESTS: + @staticmethod + def jordan_product(X,Y): + return (X*Y + Y*X)/2 + + @classmethod + def trace_inner_product(cls,X,Y): + Xu = cls.real_unembed(X) + Yu = cls.real_unembed(Y) + tr = (Xu*Yu).trace() + + try: + # Works in QQ, AA, RDF, et cetera. + return tr.real() + except AttributeError: + # A quaternion doesn't have a real() method, but does + # have coefficient_tuple() method that returns the + # coefficients of 1, i, j, and k -- in that order. + return tr.coefficient_tuple()[0] + + +class RealMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): + @staticmethod + def real_embed(M): + """ + The identity function, for embedding real matrices into real + matrices. + """ + return M - Unembedding is the inverse of embedding:: + @staticmethod + def real_unembed(M): + """ + The identity function, for unembedding real matrices from real + matrices. + """ + return M - sage: set_random_seed() - sage: Q = QuaternionAlgebra(QQ, -1, -1) - sage: M = random_matrix(Q, 3) - sage: _unembed_quaternion_matrix(_embed_quaternion_matrix(M)) == M - True - """ - n = ZZ(M.nrows()) - if M.ncols() != n: - raise ValueError("the matrix 'M' must be square") - if not n.mod(4).is_zero(): - raise ValueError("the matrix 'M' must be a complex embedding") - - Q = QuaternionAlgebra(QQ,-1,-1) - i,j,k = Q.gens() - - # Go top-left to bottom-right (reading order), converting every - # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1 - # quaternion block. - elements = [] - for l in xrange(n/4): - for m in xrange(n/4): - submat = _unembed_complex_matrix(M[4*l:4*l+4,4*m:4*m+4]) - if submat[0,0] != submat[1,1].conjugate(): - raise ValueError('bad on-diagonal submatrix') - if submat[0,1] != -submat[1,0].conjugate(): - raise ValueError('bad off-diagonal submatrix') - z = submat[0,0].real() + submat[0,0].imag()*i - z += submat[0,1].real()*j + submat[0,1].imag()*k - elements.append(z) - - return matrix(Q, n/4, elements) - - -# The usual inner product on R^n. -def _usual_ip(x,y): - return x.vector().inner_product(y.vector()) - -# The inner product used for the real symmetric simple EJA. -# We keep it as a separate function because e.g. the complex -# algebra uses the same inner product, except divided by 2. -def _matrix_ip(X,Y): - X_mat = X.natural_representation() - Y_mat = Y.natural_representation() - return (X_mat*Y_mat).trace() - - -class RealSymmetricEJA(FiniteDimensionalEuclideanJordanAlgebra): +class RealSymmetricEJA(ConcreteEuclideanJordanAlgebra, + RealMatrixEuclideanJordanAlgebra): """ The rank-n simple EJA consisting of real symmetric n-by-n matrices, the usual symmetric Jordan product, and the trace inner @@ -970,54 +1538,267 @@ class RealSymmetricEJA(FiniteDimensionalEuclideanJordanAlgebra): sage: e0*e0 e0 sage: e1*e1 - e0 + e2 + 1/2*e0 + 1/2*e2 sage: e2*e2 e2 + In theory, our "field" can be any subfield of the reals:: + + sage: RealSymmetricEJA(2, field=RDF) + Euclidean Jordan algebra of dimension 3 over Real Double Field + sage: RealSymmetricEJA(2, field=RR) + Euclidean Jordan algebra of dimension 3 over Real Field with + 53 bits of precision + TESTS: - The degree of this algebra is `(n^2 + n) / 2`:: + The dimension of this algebra is `(n^2 + n) / 2`:: sage: set_random_seed() - sage: n = ZZ.random_element(1,5) + sage: n_max = RealSymmetricEJA._max_random_instance_size() + sage: n = ZZ.random_element(1, n_max) sage: J = RealSymmetricEJA(n) - sage: J.degree() == (n^2 + n)/2 + sage: J.dimension() == (n^2 + n)/2 True The Jordan multiplication is what we think it is:: sage: set_random_seed() - sage: n = ZZ.random_element(1,5) - sage: J = RealSymmetricEJA(n) - sage: x = J.random_element() - sage: y = J.random_element() - sage: actual = (x*y).natural_representation() - sage: X = x.natural_representation() - sage: Y = y.natural_representation() + sage: J = RealSymmetricEJA.random_instance() + sage: x,y = J.random_elements(2) + sage: actual = (x*y).to_matrix() + sage: X = x.to_matrix() + sage: Y = y.to_matrix() sage: expected = (X*Y + Y*X)/2 sage: actual == expected True sage: J(expected) == x*y True + We can change the generator prefix:: + + sage: RealSymmetricEJA(3, prefix='q').gens() + (q0, q1, q2, q3, q4, q5) + + We can construct the (trivial) algebra of rank zero:: + + sage: RealSymmetricEJA(0) + Euclidean Jordan algebra of dimension 0 over Algebraic Real Field + """ + @classmethod + def _denormalized_basis(cls, n): + """ + Return a basis for the space of real symmetric n-by-n matrices. + + SETUP:: + + sage: from mjo.eja.eja_algebra import RealSymmetricEJA + + TESTS:: + + sage: set_random_seed() + sage: n = ZZ.random_element(1,5) + sage: B = RealSymmetricEJA._denormalized_basis(n) + sage: all( M.is_symmetric() for M in B) + True + + """ + # The basis of symmetric matrices, as matrices, in their R^(n-by-n) + # coordinates. + S = [] + for i in range(n): + for j in range(i+1): + Eij = matrix(ZZ, n, lambda k,l: k==i and l==j) + if i == j: + Sij = Eij + else: + Sij = Eij + Eij.transpose() + S.append(Sij) + return tuple(S) + + @staticmethod - def __classcall_private__(cls, n, field=QQ): - S = _real_symmetric_basis(n, field=field) - (Qs, T) = _multiplication_table_from_matrix_basis(S) + def _max_random_instance_size(): + return 4 # Dimension 10 - fdeja = super(RealSymmetricEJA, cls) - return fdeja.__classcall_private__(cls, - field, - Qs, - rank=n, - natural_basis=T) + @classmethod + def random_instance(cls, **kwargs): + """ + Return a random instance of this type of algebra. + """ + n = ZZ.random_element(cls._max_random_instance_size() + 1) + return cls(n, **kwargs) - def inner_product(self, x, y): - return _matrix_ip(x,y) + def __init__(self, n, **kwargs): + super(RealSymmetricEJA, self).__init__(self._denormalized_basis(n), + self.jordan_product, + self.trace_inner_product, + **kwargs) + self.rank.set_cache(n) + self.one.set_cache(self(matrix.identity(ZZ,n))) -class ComplexHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra): +class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): + @staticmethod + def real_embed(M): + """ + Embed the n-by-n complex matrix ``M`` into the space of real + matrices of size 2n-by-2n via the map the sends each entry `z = a + + bi` to the block matrix ``[[a,b],[-b,a]]``. + + SETUP:: + + sage: from mjo.eja.eja_algebra import \ + ....: ComplexMatrixEuclideanJordanAlgebra + + EXAMPLES:: + + sage: F = QuadraticField(-1, 'I') + sage: x1 = F(4 - 2*i) + sage: x2 = F(1 + 2*i) + sage: x3 = F(-i) + sage: x4 = F(6) + sage: M = matrix(F,2,[[x1,x2],[x3,x4]]) + sage: ComplexMatrixEuclideanJordanAlgebra.real_embed(M) + [ 4 -2| 1 2] + [ 2 4|-2 1] + [-----+-----] + [ 0 -1| 6 0] + [ 1 0| 0 6] + + TESTS: + + Embedding is a homomorphism (isomorphism, in fact):: + + sage: set_random_seed() + sage: n = ZZ.random_element(3) + sage: F = QuadraticField(-1, 'I') + sage: X = random_matrix(F, n) + sage: Y = random_matrix(F, n) + sage: Xe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X) + sage: Ye = ComplexMatrixEuclideanJordanAlgebra.real_embed(Y) + sage: XYe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X*Y) + sage: Xe*Ye == XYe + True + + """ + n = M.nrows() + if M.ncols() != n: + raise ValueError("the matrix 'M' must be square") + + # We don't need any adjoined elements... + field = M.base_ring().base_ring() + + blocks = [] + for z in M.list(): + a = z.list()[0] # real part, I guess + b = z.list()[1] # imag part, I guess + blocks.append(matrix(field, 2, [[a,b],[-b,a]])) + + return matrix.block(field, n, blocks) + + + @staticmethod + def real_unembed(M): + """ + The inverse of _embed_complex_matrix(). + + SETUP:: + + sage: from mjo.eja.eja_algebra import \ + ....: ComplexMatrixEuclideanJordanAlgebra + + EXAMPLES:: + + sage: A = matrix(QQ,[ [ 1, 2, 3, 4], + ....: [-2, 1, -4, 3], + ....: [ 9, 10, 11, 12], + ....: [-10, 9, -12, 11] ]) + sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(A) + [ 2*I + 1 4*I + 3] + [ 10*I + 9 12*I + 11] + + TESTS: + + Unembedding is the inverse of embedding:: + + sage: set_random_seed() + sage: F = QuadraticField(-1, 'I') + sage: M = random_matrix(F, 3) + sage: Me = ComplexMatrixEuclideanJordanAlgebra.real_embed(M) + sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(Me) == M + True + + """ + n = ZZ(M.nrows()) + if M.ncols() != n: + raise ValueError("the matrix 'M' must be square") + if not n.mod(2).is_zero(): + raise ValueError("the matrix 'M' must be a complex embedding") + + # If "M" was normalized, its base ring might have roots + # adjoined and they can stick around after unembedding. + field = M.base_ring() + R = PolynomialRing(field, 'z') + z = R.gen() + if field is AA: + # Sage doesn't know how to embed AA into QQbar, i.e. how + # to adjoin sqrt(-1) to AA. + F = QQbar + else: + F = field.extension(z**2 + 1, 'I', embedding=CLF(-1).sqrt()) + i = F.gen() + + # Go top-left to bottom-right (reading order), converting every + # 2-by-2 block we see to a single complex element. + elements = [] + for k in range(n/2): + for j in range(n/2): + submat = M[2*k:2*k+2,2*j:2*j+2] + if submat[0,0] != submat[1,1]: + raise ValueError('bad on-diagonal submatrix') + if submat[0,1] != -submat[1,0]: + raise ValueError('bad off-diagonal submatrix') + z = submat[0,0] + submat[0,1]*i + elements.append(z) + + return matrix(F, n/2, elements) + + + @classmethod + def trace_inner_product(cls,X,Y): + """ + Compute a matrix inner product in this algebra directly from + its real embedding. + + SETUP:: + + sage: from mjo.eja.eja_algebra import ComplexHermitianEJA + + TESTS: + + This gives the same answer as the slow, default method implemented + in :class:`MatrixEuclideanJordanAlgebra`:: + + sage: set_random_seed() + sage: J = ComplexHermitianEJA.random_instance() + sage: x,y = J.random_elements(2) + sage: Xe = x.to_matrix() + sage: Ye = y.to_matrix() + sage: X = ComplexHermitianEJA.real_unembed(Xe) + sage: Y = ComplexHermitianEJA.real_unembed(Ye) + sage: expected = (X*Y).trace().real() + sage: actual = ComplexHermitianEJA.trace_inner_product(Xe,Ye) + sage: actual == expected + True + + """ + return RealMatrixEuclideanJordanAlgebra.trace_inner_product(X,Y)/2 + + +class ComplexHermitianEJA(ConcreteEuclideanJordanAlgebra, + ComplexMatrixEuclideanJordanAlgebra): """ The rank-n simple EJA consisting of complex Hermitian n-by-n matrices over the real numbers, the usual symmetric Jordan product, @@ -1028,58 +1809,293 @@ class ComplexHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra): sage: from mjo.eja.eja_algebra import ComplexHermitianEJA + EXAMPLES: + + In theory, our "field" can be any subfield of the reals:: + + sage: ComplexHermitianEJA(2, field=RDF) + Euclidean Jordan algebra of dimension 4 over Real Double Field + sage: ComplexHermitianEJA(2, field=RR) + Euclidean Jordan algebra of dimension 4 over Real Field with + 53 bits of precision + TESTS: - The degree of this algebra is `n^2`:: + The dimension of this algebra is `n^2`:: sage: set_random_seed() - sage: n = ZZ.random_element(1,5) + sage: n_max = ComplexHermitianEJA._max_random_instance_size() + sage: n = ZZ.random_element(1, n_max) sage: J = ComplexHermitianEJA(n) - sage: J.degree() == n^2 + sage: J.dimension() == n^2 True The Jordan multiplication is what we think it is:: sage: set_random_seed() - sage: n = ZZ.random_element(1,5) - sage: J = ComplexHermitianEJA(n) - sage: x = J.random_element() - sage: y = J.random_element() - sage: actual = (x*y).natural_representation() - sage: X = x.natural_representation() - sage: Y = y.natural_representation() + sage: J = ComplexHermitianEJA.random_instance() + sage: x,y = J.random_elements(2) + sage: actual = (x*y).to_matrix() + sage: X = x.to_matrix() + sage: Y = y.to_matrix() sage: expected = (X*Y + Y*X)/2 sage: actual == expected True sage: J(expected) == x*y True + We can change the generator prefix:: + + sage: ComplexHermitianEJA(2, prefix='z').gens() + (z0, z1, z2, z3) + + We can construct the (trivial) algebra of rank zero:: + + sage: ComplexHermitianEJA(0) + Euclidean Jordan algebra of dimension 0 over Algebraic Real Field + """ - @staticmethod - def __classcall_private__(cls, n, field=QQ): - S = _complex_hermitian_basis(n) - (Qs, T) = _multiplication_table_from_matrix_basis(S) - fdeja = super(ComplexHermitianEJA, cls) - return fdeja.__classcall_private__(cls, - field, - Qs, - rank=n, - natural_basis=T) + @classmethod + def _denormalized_basis(cls, n): + """ + Returns a basis for the space of complex Hermitian n-by-n matrices. - def inner_product(self, x, y): - # Since a+bi on the diagonal is represented as + Why do we embed these? Basically, because all of numerical linear + algebra assumes that you're working with vectors consisting of `n` + entries from a field and scalars from the same field. There's no way + to tell SageMath that (for example) the vectors contain complex + numbers, while the scalar field is real. + + SETUP:: + + sage: from mjo.eja.eja_algebra import ComplexHermitianEJA + + TESTS:: + + sage: set_random_seed() + sage: n = ZZ.random_element(1,5) + sage: field = QuadraticField(2, 'sqrt2') + sage: B = ComplexHermitianEJA._denormalized_basis(n) + sage: all( M.is_symmetric() for M in B) + True + + """ + field = ZZ + R = PolynomialRing(field, 'z') + z = R.gen() + F = field.extension(z**2 + 1, 'I') + I = F.gen(1) + + # This is like the symmetric case, but we need to be careful: # - # a + bi = [ a b ] - # [ -b a ], + # * We want conjugate-symmetry, not just symmetry. + # * The diagonal will (as a result) be real. # - # we'll double-count the "a" entries if we take the trace of - # the embedding. - return _matrix_ip(x,y)/2 + S = [] + for i in range(n): + for j in range(i+1): + Eij = matrix(F, n, lambda k,l: k==i and l==j) + if i == j: + Sij = cls.real_embed(Eij) + S.append(Sij) + else: + # The second one has a minus because it's conjugated. + Sij_real = cls.real_embed(Eij + Eij.transpose()) + S.append(Sij_real) + Sij_imag = cls.real_embed(I*Eij - I*Eij.transpose()) + S.append(Sij_imag) + + # Since we embedded these, we can drop back to the "field" that we + # started with instead of the complex extension "F". + return tuple( s.change_ring(field) for s in S ) + + + def __init__(self, n, **kwargs): + super(ComplexHermitianEJA, self).__init__(self._denormalized_basis(n), + self.jordan_product, + self.trace_inner_product, + **kwargs) + self.rank.set_cache(n) + # TODO: pre-cache the identity! + + @staticmethod + def _max_random_instance_size(): + return 3 # Dimension 9 + @classmethod + def random_instance(cls, **kwargs): + """ + Return a random instance of this type of algebra. + """ + n = ZZ.random_element(cls._max_random_instance_size() + 1) + return cls(n, **kwargs) -class QuaternionHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra): - """ +class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): + @staticmethod + def real_embed(M): + """ + Embed the n-by-n quaternion matrix ``M`` into the space of real + matrices of size 4n-by-4n by first sending each quaternion entry `z + = a + bi + cj + dk` to the block-complex matrix ``[[a + bi, + c+di],[-c + di, a-bi]]`, and then embedding those into a real + matrix. + + SETUP:: + + sage: from mjo.eja.eja_algebra import \ + ....: QuaternionMatrixEuclideanJordanAlgebra + + EXAMPLES:: + + sage: Q = QuaternionAlgebra(QQ,-1,-1) + sage: i,j,k = Q.gens() + sage: x = 1 + 2*i + 3*j + 4*k + sage: M = matrix(Q, 1, [[x]]) + sage: QuaternionMatrixEuclideanJordanAlgebra.real_embed(M) + [ 1 2 3 4] + [-2 1 -4 3] + [-3 4 1 -2] + [-4 -3 2 1] + + Embedding is a homomorphism (isomorphism, in fact):: + + sage: set_random_seed() + sage: n = ZZ.random_element(2) + sage: Q = QuaternionAlgebra(QQ,-1,-1) + sage: X = random_matrix(Q, n) + sage: Y = random_matrix(Q, n) + sage: Xe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X) + sage: Ye = QuaternionMatrixEuclideanJordanAlgebra.real_embed(Y) + sage: XYe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X*Y) + sage: Xe*Ye == XYe + True + + """ + quaternions = M.base_ring() + n = M.nrows() + if M.ncols() != n: + raise ValueError("the matrix 'M' must be square") + + F = QuadraticField(-1, 'I') + i = F.gen() + + blocks = [] + for z in M.list(): + t = z.coefficient_tuple() + a = t[0] + b = t[1] + c = t[2] + d = t[3] + cplxM = matrix(F, 2, [[ a + b*i, c + d*i], + [-c + d*i, a - b*i]]) + realM = ComplexMatrixEuclideanJordanAlgebra.real_embed(cplxM) + blocks.append(realM) + + # We should have real entries by now, so use the realest field + # we've got for the return value. + return matrix.block(quaternions.base_ring(), n, blocks) + + + + @staticmethod + def real_unembed(M): + """ + The inverse of _embed_quaternion_matrix(). + + SETUP:: + + sage: from mjo.eja.eja_algebra import \ + ....: QuaternionMatrixEuclideanJordanAlgebra + + EXAMPLES:: + + sage: M = matrix(QQ, [[ 1, 2, 3, 4], + ....: [-2, 1, -4, 3], + ....: [-3, 4, 1, -2], + ....: [-4, -3, 2, 1]]) + sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(M) + [1 + 2*i + 3*j + 4*k] + + TESTS: + + Unembedding is the inverse of embedding:: + + sage: set_random_seed() + sage: Q = QuaternionAlgebra(QQ, -1, -1) + sage: M = random_matrix(Q, 3) + sage: Me = QuaternionMatrixEuclideanJordanAlgebra.real_embed(M) + sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(Me) == M + True + + """ + n = ZZ(M.nrows()) + if M.ncols() != n: + raise ValueError("the matrix 'M' must be square") + if not n.mod(4).is_zero(): + raise ValueError("the matrix 'M' must be a quaternion embedding") + + # Use the base ring of the matrix to ensure that its entries can be + # multiplied by elements of the quaternion algebra. + field = M.base_ring() + Q = QuaternionAlgebra(field,-1,-1) + i,j,k = Q.gens() + + # Go top-left to bottom-right (reading order), converting every + # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1 + # quaternion block. + elements = [] + for l in range(n/4): + for m in range(n/4): + submat = ComplexMatrixEuclideanJordanAlgebra.real_unembed( + M[4*l:4*l+4,4*m:4*m+4] ) + if submat[0,0] != submat[1,1].conjugate(): + raise ValueError('bad on-diagonal submatrix') + if submat[0,1] != -submat[1,0].conjugate(): + raise ValueError('bad off-diagonal submatrix') + z = submat[0,0].real() + z += submat[0,0].imag()*i + z += submat[0,1].real()*j + z += submat[0,1].imag()*k + elements.append(z) + + return matrix(Q, n/4, elements) + + + @classmethod + def trace_inner_product(cls,X,Y): + """ + Compute a matrix inner product in this algebra directly from + its real embedding. + + SETUP:: + + sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA + + TESTS: + + This gives the same answer as the slow, default method implemented + in :class:`MatrixEuclideanJordanAlgebra`:: + + sage: set_random_seed() + sage: J = QuaternionHermitianEJA.random_instance() + sage: x,y = J.random_elements(2) + sage: Xe = x.to_matrix() + sage: Ye = y.to_matrix() + sage: X = QuaternionHermitianEJA.real_unembed(Xe) + sage: Y = QuaternionHermitianEJA.real_unembed(Ye) + sage: expected = (X*Y).trace().coefficient_tuple()[0] + sage: actual = QuaternionHermitianEJA.trace_inner_product(Xe,Ye) + sage: actual == expected + True + + """ + return RealMatrixEuclideanJordanAlgebra.trace_inner_product(X,Y)/4 + + +class QuaternionHermitianEJA(ConcreteEuclideanJordanAlgebra, + QuaternionMatrixEuclideanJordanAlgebra): + r""" The rank-n simple EJA consisting of self-adjoint n-by-n quaternion matrices, the usual symmetric Jordan product, and the real-part-of-trace inner product. It has dimension `2n^2 - n` over @@ -1089,63 +2105,368 @@ class QuaternionHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra): sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA + EXAMPLES: + + In theory, our "field" can be any subfield of the reals:: + + sage: QuaternionHermitianEJA(2, field=RDF) + Euclidean Jordan algebra of dimension 6 over Real Double Field + sage: QuaternionHermitianEJA(2, field=RR) + Euclidean Jordan algebra of dimension 6 over Real Field with + 53 bits of precision + TESTS: - The degree of this algebra is `n^2`:: + The dimension of this algebra is `2*n^2 - n`:: sage: set_random_seed() - sage: n = ZZ.random_element(1,5) + sage: n_max = QuaternionHermitianEJA._max_random_instance_size() + sage: n = ZZ.random_element(1, n_max) sage: J = QuaternionHermitianEJA(n) - sage: J.degree() == 2*(n^2) - n + sage: J.dimension() == 2*(n^2) - n True The Jordan multiplication is what we think it is:: sage: set_random_seed() - sage: n = ZZ.random_element(1,5) - sage: J = QuaternionHermitianEJA(n) - sage: x = J.random_element() - sage: y = J.random_element() - sage: actual = (x*y).natural_representation() - sage: X = x.natural_representation() - sage: Y = y.natural_representation() + sage: J = QuaternionHermitianEJA.random_instance() + sage: x,y = J.random_elements(2) + sage: actual = (x*y).to_matrix() + sage: X = x.to_matrix() + sage: Y = y.to_matrix() sage: expected = (X*Y + Y*X)/2 sage: actual == expected True sage: J(expected) == x*y True + We can change the generator prefix:: + + sage: QuaternionHermitianEJA(2, prefix='a').gens() + (a0, a1, a2, a3, a4, a5) + + We can construct the (trivial) algebra of rank zero:: + + sage: QuaternionHermitianEJA(0) + Euclidean Jordan algebra of dimension 0 over Algebraic Real Field + """ - @staticmethod - def __classcall_private__(cls, n, field=QQ): - S = _quaternion_hermitian_basis(n) - (Qs, T) = _multiplication_table_from_matrix_basis(S) + @classmethod + def _denormalized_basis(cls, n): + """ + Returns a basis for the space of quaternion Hermitian n-by-n matrices. - fdeja = super(QuaternionHermitianEJA, cls) - return fdeja.__classcall_private__(cls, - field, - Qs, - rank=n, - natural_basis=T) + Why do we embed these? Basically, because all of numerical + linear algebra assumes that you're working with vectors consisting + of `n` entries from a field and scalars from the same field. There's + no way to tell SageMath that (for example) the vectors contain + complex numbers, while the scalar field is real. - def inner_product(self, x, y): - # Since a+bi+cj+dk on the diagonal is represented as + SETUP:: + + sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA + + TESTS:: + + sage: set_random_seed() + sage: n = ZZ.random_element(1,5) + sage: B = QuaternionHermitianEJA._denormalized_basis(n) + sage: all( M.is_symmetric() for M in B ) + True + + """ + field = ZZ + Q = QuaternionAlgebra(QQ,-1,-1) + I,J,K = Q.gens() + + # This is like the symmetric case, but we need to be careful: # - # a + bi +cj + dk = [ a b c d] - # [ -b a -d c] - # [ -c d a -b] - # [ -d -c b a], + # * We want conjugate-symmetry, not just symmetry. + # * The diagonal will (as a result) be real. # - # we'll quadruple-count the "a" entries if we take the trace of - # the embedding. - return _matrix_ip(x,y)/4 + S = [] + for i in range(n): + for j in range(i+1): + Eij = matrix(Q, n, lambda k,l: k==i and l==j) + if i == j: + Sij = cls.real_embed(Eij) + S.append(Sij) + else: + # The second, third, and fourth ones have a minus + # because they're conjugated. + Sij_real = cls.real_embed(Eij + Eij.transpose()) + S.append(Sij_real) + Sij_I = cls.real_embed(I*Eij - I*Eij.transpose()) + S.append(Sij_I) + Sij_J = cls.real_embed(J*Eij - J*Eij.transpose()) + S.append(Sij_J) + Sij_K = cls.real_embed(K*Eij - K*Eij.transpose()) + S.append(Sij_K) + + # Since we embedded these, we can drop back to the "field" that we + # started with instead of the quaternion algebra "Q". + return tuple( s.change_ring(field) for s in S ) + + + def __init__(self, n, **kwargs): + super(QuaternionHermitianEJA, self).__init__(self._denormalized_basis(n), + self.jordan_product, + self.trace_inner_product, + **kwargs) + self.rank.set_cache(n) + # TODO: cache one()! + + @staticmethod + def _max_random_instance_size(): + r""" + The maximum rank of a random QuaternionHermitianEJA. + """ + return 2 # Dimension 6 + + @classmethod + def random_instance(cls, **kwargs): + """ + Return a random instance of this type of algebra. + """ + n = ZZ.random_element(cls._max_random_instance_size() + 1) + return cls(n, **kwargs) -class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra): +class HadamardEJA(ConcreteEuclideanJordanAlgebra): + """ + Return the Euclidean Jordan Algebra corresponding to the set + `R^n` under the Hadamard product. + + Note: this is nothing more than the Cartesian product of ``n`` + copies of the spin algebra. Once Cartesian product algebras + are implemented, this can go. + + SETUP:: + + sage: from mjo.eja.eja_algebra import HadamardEJA + + EXAMPLES: + + This multiplication table can be verified by hand:: + + sage: J = HadamardEJA(3) + sage: e0,e1,e2 = J.gens() + sage: e0*e0 + e0 + sage: e0*e1 + 0 + sage: e0*e2 + 0 + sage: e1*e1 + e1 + sage: e1*e2 + 0 + sage: e2*e2 + e2 + + TESTS: + + We can change the generator prefix:: + + sage: HadamardEJA(3, prefix='r').gens() + (r0, r1, r2) + + """ + def __init__(self, n, **kwargs): + def jordan_product(x,y): + P = x.parent() + return P(tuple( xi*yi for (xi,yi) in zip(x,y) )) + def inner_product(x,y): + return x.inner_product(y) + + # Don't orthonormalize because our basis is already + # orthonormal with respect to our inner-product. + if not 'orthonormalize' in kwargs: + kwargs['orthonormalize'] = False + + # But also don't pass check_field=False here, because the user + # can pass in a field! + standard_basis = FreeModule(ZZ, n).basis() + super(HadamardEJA, self).__init__(standard_basis, + jordan_product, + inner_product, + check_axioms=False, + **kwargs) + self.rank.set_cache(n) + + if n == 0: + self.one.set_cache( self.zero() ) + else: + self.one.set_cache( sum(self.gens()) ) + + @staticmethod + def _max_random_instance_size(): + r""" + The maximum dimension of a random HadamardEJA. + """ + return 5 + + @classmethod + def random_instance(cls, **kwargs): + """ + Return a random instance of this type of algebra. + """ + n = ZZ.random_element(cls._max_random_instance_size() + 1) + return cls(n, **kwargs) + + +class BilinearFormEJA(ConcreteEuclideanJordanAlgebra): + r""" + The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)`` + with the half-trace inner product and jordan product ``x*y = + (,y_bar>, x0*y_bar + y0*x_bar)`` where `B = 1 \times B22` is + a symmetric positive-definite "bilinear form" matrix. Its + dimension is the size of `B`, and it has rank two in dimensions + larger than two. It reduces to the ``JordanSpinEJA`` when `B` is + the identity matrix of order ``n``. + + We insist that the one-by-one upper-left identity block of `B` be + passed in as well so that we can be passed a matrix of size zero + to construct a trivial algebra. + + SETUP:: + + sage: from mjo.eja.eja_algebra import (BilinearFormEJA, + ....: JordanSpinEJA) + + EXAMPLES: + + When no bilinear form is specified, the identity matrix is used, + and the resulting algebra is the Jordan spin algebra:: + + sage: B = matrix.identity(AA,3) + sage: J0 = BilinearFormEJA(B) + sage: J1 = JordanSpinEJA(3) + sage: J0.multiplication_table() == J0.multiplication_table() + True + + An error is raised if the matrix `B` does not correspond to a + positive-definite bilinear form:: + + sage: B = matrix.random(QQ,2,3) + sage: J = BilinearFormEJA(B) + Traceback (most recent call last): + ... + ValueError: bilinear form is not positive-definite + sage: B = matrix.zero(QQ,3) + sage: J = BilinearFormEJA(B) + Traceback (most recent call last): + ... + ValueError: bilinear form is not positive-definite + + TESTS: + + We can create a zero-dimensional algebra:: + + sage: B = matrix.identity(AA,0) + sage: J = BilinearFormEJA(B) + sage: J.basis() + Finite family {} + + We can check the multiplication condition given in the Jordan, von + Neumann, and Wigner paper (and also discussed on my "On the + symmetry..." paper). Note that this relies heavily on the standard + choice of basis, as does anything utilizing the bilinear form + matrix. We opt not to orthonormalize the basis, because if we + did, we would have to normalize the `s_{i}` in a similar manner:: + + sage: set_random_seed() + sage: n = ZZ.random_element(5) + sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular') + sage: B11 = matrix.identity(QQ,1) + sage: B22 = M.transpose()*M + sage: B = block_matrix(2,2,[ [B11,0 ], + ....: [0, B22 ] ]) + sage: J = BilinearFormEJA(B, orthonormalize=False) + sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis() + sage: V = J.vector_space() + sage: sis = [ J( V([0] + (M.inverse()*ei).list()).column() ) + ....: for ei in eis ] + sage: actual = [ sis[i]*sis[j] + ....: for i in range(n-1) + ....: for j in range(n-1) ] + sage: expected = [ J.one() if i == j else J.zero() + ....: for i in range(n-1) + ....: for j in range(n-1) ] + sage: actual == expected + True + """ + def __init__(self, B, **kwargs): + if not B.is_positive_definite(): + raise ValueError("bilinear form is not positive-definite") + + def inner_product(x,y): + return (B*x).inner_product(y) + + def jordan_product(x,y): + P = x.parent() + x0 = x[0] + xbar = x[1:] + y0 = y[0] + ybar = y[1:] + z0 = inner_product(x,y) + zbar = y0*xbar + x0*ybar + return P((z0,) + tuple(zbar)) + + n = B.nrows() + standard_basis = FreeModule(ZZ, n).basis() + super(BilinearFormEJA, self).__init__(standard_basis, + jordan_product, + inner_product, + **kwargs) + + # The rank of this algebra is two, unless we're in a + # one-dimensional ambient space (because the rank is bounded + # by the ambient dimension). + self.rank.set_cache(min(n,2)) + + if n == 0: + self.one.set_cache( self.zero() ) + else: + self.one.set_cache( self.monomial(0) ) + + @staticmethod + def _max_random_instance_size(): + r""" + The maximum dimension of a random BilinearFormEJA. + """ + return 5 + + @classmethod + def random_instance(cls, **kwargs): + """ + Return a random instance of this algebra. + """ + n = ZZ.random_element(cls._max_random_instance_size() + 1) + if n.is_zero(): + B = matrix.identity(ZZ, n) + return cls(B, **kwargs) + + B11 = matrix.identity(ZZ, 1) + M = matrix.random(ZZ, n-1) + I = matrix.identity(ZZ, n-1) + alpha = ZZ.zero() + while alpha.is_zero(): + alpha = ZZ.random_element().abs() + B22 = M.transpose()*M + alpha*I + + from sage.matrix.special import block_matrix + B = block_matrix(2,2, [ [B11, ZZ(0) ], + [ZZ(0), B22 ] ]) + + return cls(B, **kwargs) + + +class JordanSpinEJA(BilinearFormEJA): """ The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)`` with the usual inner product and jordan product ``x*y = - (, x0*y_bar + y0*x_bar)``. It has dimension `n` over + (, x0*y_bar + y0*x_bar)``. It has dimension `n` over the reals. SETUP:: @@ -1173,27 +2494,344 @@ class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra): sage: e2*e3 0 + We can change the generator prefix:: + + sage: JordanSpinEJA(2, prefix='B').gens() + (B0, B1) + + TESTS: + + Ensure that we have the usual inner product on `R^n`:: + + sage: set_random_seed() + sage: J = JordanSpinEJA.random_instance() + sage: x,y = J.random_elements(2) + sage: actual = x.inner_product(y) + sage: expected = x.to_vector().inner_product(y.to_vector()) + sage: actual == expected + True + """ + def __init__(self, n, **kwargs): + # This is a special case of the BilinearFormEJA with the + # identity matrix as its bilinear form. + B = matrix.identity(ZZ, n) + + # Don't orthonormalize because our basis is already + # orthonormal with respect to our inner-product. + if not 'orthonormalize' in kwargs: + kwargs['orthonormalize'] = False + + # But also don't pass check_field=False here, because the user + # can pass in a field! + super(JordanSpinEJA, self).__init__(B, + check_axioms=False, + **kwargs) + @staticmethod - def __classcall_private__(cls, n, field=QQ): - Qs = [] - id_matrix = matrix.identity(field, n) - for i in xrange(n): - ei = id_matrix.column(i) - Qi = matrix.zero(field, n) - Qi.set_row(0, ei) - Qi.set_column(0, ei) - Qi += matrix.diagonal(n, [ei[0]]*n) - # The addition of the diagonal matrix adds an extra ei[0] in the - # upper-left corner of the matrix. - Qi[0,0] = Qi[0,0] * ~field(2) - Qs.append(Qi) - - # The rank of the spin algebra is two, unless we're in a - # one-dimensional ambient space (because the rank is bounded by - # the ambient dimension). - fdeja = super(JordanSpinEJA, cls) - return fdeja.__classcall_private__(cls, field, Qs, rank=min(n,2)) + def _max_random_instance_size(): + r""" + The maximum dimension of a random JordanSpinEJA. + """ + return 5 + + @classmethod + def random_instance(cls, **kwargs): + """ + Return a random instance of this type of algebra. + + Needed here to override the implementation for ``BilinearFormEJA``. + """ + n = ZZ.random_element(cls._max_random_instance_size() + 1) + return cls(n, **kwargs) + + +class TrivialEJA(ConcreteEuclideanJordanAlgebra): + """ + The trivial Euclidean Jordan algebra consisting of only a zero element. + + SETUP:: + + sage: from mjo.eja.eja_algebra import TrivialEJA + + EXAMPLES:: + + sage: J = TrivialEJA() + sage: J.dimension() + 0 + sage: J.zero() + 0 + sage: J.one() + 0 + sage: 7*J.one()*12*J.one() + 0 + sage: J.one().inner_product(J.one()) + 0 + sage: J.one().norm() + 0 + sage: J.one().subalgebra_generated_by() + Euclidean Jordan algebra of dimension 0 over Algebraic Real Field + sage: J.rank() + 0 + + """ + def __init__(self, **kwargs): + jordan_product = lambda x,y: x + inner_product = lambda x,y: 0 + basis = () + super(TrivialEJA, self).__init__(basis, + jordan_product, + inner_product, + **kwargs) + # The rank is zero using my definition, namely the dimension of the + # largest subalgebra generated by any element. + self.rank.set_cache(0) + self.one.set_cache( self.zero() ) + + @classmethod + def random_instance(cls, **kwargs): + # We don't take a "size" argument so the superclass method is + # inappropriate for us. + return cls(**kwargs) + +class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra): + r""" + The external (orthogonal) direct sum of two other Euclidean Jordan + algebras. Essentially the Cartesian product of its two factors. + Every Euclidean Jordan algebra decomposes into an orthogonal + direct sum of simple Euclidean Jordan algebras, so no generality + is lost by providing only this construction. + + SETUP:: + + sage: from mjo.eja.eja_algebra import (random_eja, + ....: HadamardEJA, + ....: RealSymmetricEJA, + ....: DirectSumEJA) + + EXAMPLES:: + + sage: J1 = HadamardEJA(2) + sage: J2 = RealSymmetricEJA(3) + sage: J = DirectSumEJA(J1,J2) + sage: J.dimension() + 8 + sage: J.rank() + 5 + + TESTS: + + The external direct sum construction is only valid when the two factors + have the same base ring; an error is raised otherwise:: + + sage: set_random_seed() + sage: J1 = random_eja(field=AA) + sage: J2 = random_eja(field=QQ,orthonormalize=False) + sage: J = DirectSumEJA(J1,J2) + Traceback (most recent call last): + ... + ValueError: algebras must share the same base field + + """ + def __init__(self, J1, J2, **kwargs): + if J1.base_ring() != J2.base_ring(): + raise ValueError("algebras must share the same base field") + field = J1.base_ring() + + self._factors = (J1, J2) + n1 = J1.dimension() + n2 = J2.dimension() + n = n1+n2 + V = VectorSpace(field, n) + mult_table = [ [ V.zero() for j in range(i+1) ] + for i in range(n) ] + for i in range(n1): + for j in range(i+1): + p = (J1.monomial(i)*J1.monomial(j)).to_vector() + mult_table[i][j] = V(p.list() + [field.zero()]*n2) + + for i in range(n2): + for j in range(i+1): + p = (J2.monomial(i)*J2.monomial(j)).to_vector() + mult_table[n1+i][n1+j] = V([field.zero()]*n1 + p.list()) + + # TODO: build the IP table here from the two constituent IP + # matrices (it'll be block diagonal, I think). + ip_table = [ [ field.zero() for j in range(i+1) ] + for i in range(n) ] + super(DirectSumEJA, self).__init__(field, + mult_table, + ip_table, + check_axioms=False, + **kwargs) + self.rank.set_cache(J1.rank() + J2.rank()) + + + def factors(self): + r""" + Return the pair of this algebra's factors. + + SETUP:: + + sage: from mjo.eja.eja_algebra import (HadamardEJA, + ....: JordanSpinEJA, + ....: DirectSumEJA) + + EXAMPLES:: + + sage: J1 = HadamardEJA(2, field=QQ) + sage: J2 = JordanSpinEJA(3, field=QQ) + sage: J = DirectSumEJA(J1,J2) + sage: J.factors() + (Euclidean Jordan algebra of dimension 2 over Rational Field, + Euclidean Jordan algebra of dimension 3 over Rational Field) + + """ + return self._factors + + def projections(self): + r""" + Return a pair of projections onto this algebra's factors. + + SETUP:: + + sage: from mjo.eja.eja_algebra import (JordanSpinEJA, + ....: ComplexHermitianEJA, + ....: DirectSumEJA) + + EXAMPLES:: + + sage: J1 = JordanSpinEJA(2) + sage: J2 = ComplexHermitianEJA(2) + sage: J = DirectSumEJA(J1,J2) + sage: (pi_left, pi_right) = J.projections() + sage: J.one().to_vector() + (1, 0, 1, 0, 0, 1) + sage: pi_left(J.one()).to_vector() + (1, 0) + sage: pi_right(J.one()).to_vector() + (1, 0, 0, 1) + + """ + (J1,J2) = self.factors() + m = J1.dimension() + n = J2.dimension() + V_basis = self.vector_space().basis() + # Need to specify the dimensions explicitly so that we don't + # wind up with a zero-by-zero matrix when we want e.g. a + # zero-by-two matrix (important for composing things). + P1 = matrix(self.base_ring(), m, m+n, V_basis[:m]) + P2 = matrix(self.base_ring(), n, m+n, V_basis[m:]) + pi_left = FiniteDimensionalEuclideanJordanAlgebraOperator(self,J1,P1) + pi_right = FiniteDimensionalEuclideanJordanAlgebraOperator(self,J2,P2) + return (pi_left, pi_right) + + def inclusions(self): + r""" + Return the pair of inclusion maps from our factors into us. + + SETUP:: + + sage: from mjo.eja.eja_algebra import (random_eja, + ....: JordanSpinEJA, + ....: RealSymmetricEJA, + ....: DirectSumEJA) + + EXAMPLES:: + + sage: J1 = JordanSpinEJA(3) + sage: J2 = RealSymmetricEJA(2) + sage: J = DirectSumEJA(J1,J2) + sage: (iota_left, iota_right) = J.inclusions() + sage: iota_left(J1.zero()) == J.zero() + True + sage: iota_right(J2.zero()) == J.zero() + True + sage: J1.one().to_vector() + (1, 0, 0) + sage: iota_left(J1.one()).to_vector() + (1, 0, 0, 0, 0, 0) + sage: J2.one().to_vector() + (1, 0, 1) + sage: iota_right(J2.one()).to_vector() + (0, 0, 0, 1, 0, 1) + sage: J.one().to_vector() + (1, 0, 0, 1, 0, 1) + + TESTS: + + Composing a projection with the corresponding inclusion should + produce the identity map, and mismatching them should produce + the zero map:: + + sage: set_random_seed() + sage: J1 = random_eja() + sage: J2 = random_eja() + sage: J = DirectSumEJA(J1,J2) + sage: (iota_left, iota_right) = J.inclusions() + sage: (pi_left, pi_right) = J.projections() + sage: pi_left*iota_left == J1.one().operator() + True + sage: pi_right*iota_right == J2.one().operator() + True + sage: (pi_left*iota_right).is_zero() + True + sage: (pi_right*iota_left).is_zero() + True + + """ + (J1,J2) = self.factors() + m = J1.dimension() + n = J2.dimension() + V_basis = self.vector_space().basis() + # Need to specify the dimensions explicitly so that we don't + # wind up with a zero-by-zero matrix when we want e.g. a + # two-by-zero matrix (important for composing things). + I1 = matrix.column(self.base_ring(), m, m+n, V_basis[:m]) + I2 = matrix.column(self.base_ring(), n, m+n, V_basis[m:]) + iota_left = FiniteDimensionalEuclideanJordanAlgebraOperator(J1,self,I1) + iota_right = FiniteDimensionalEuclideanJordanAlgebraOperator(J2,self,I2) + return (iota_left, iota_right) def inner_product(self, x, y): - return _usual_ip(x,y) + r""" + The standard Cartesian inner-product. + + We project ``x`` and ``y`` onto our factors, and add up the + inner-products from the subalgebras. + + SETUP:: + + + sage: from mjo.eja.eja_algebra import (HadamardEJA, + ....: QuaternionHermitianEJA, + ....: DirectSumEJA) + + EXAMPLE:: + + sage: J1 = HadamardEJA(3,field=QQ) + sage: J2 = QuaternionHermitianEJA(2,field=QQ,orthonormalize=False) + sage: J = DirectSumEJA(J1,J2) + sage: x1 = J1.one() + sage: x2 = x1 + sage: y1 = J2.one() + sage: y2 = y1 + sage: x1.inner_product(x2) + 3 + sage: y1.inner_product(y2) + 2 + sage: J.one().inner_product(J.one()) + 5 + + """ + (pi_left, pi_right) = self.projections() + x1 = pi_left(x) + x2 = pi_right(x) + y1 = pi_left(y) + y2 = pi_right(y) + + return (x1.inner_product(y1) + x2.inner_product(y2)) + + + +random_eja = ConcreteEuclideanJordanAlgebra.random_instance