X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=mjo%2Feja%2Feja_algebra.py;h=44d83147bbeceb4acae52d573749545dec3c4b99;hb=9708043704809263d2bd543de38c46b458b873cb;hp=7c2eda211232ceac757a6d6caa40ae6722528349;hpb=33b5476e1422a36972979e768f1829cec1b421a5;p=sage.d.git diff --git a/mjo/eja/eja_algebra.py b/mjo/eja/eja_algebra.py index 7c2eda2..44d8314 100644 --- a/mjo/eja/eja_algebra.py +++ b/mjo/eja/eja_algebra.py @@ -3,9 +3,20 @@ 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 itertools import izip, repeat +from itertools import repeat from sage.algebras.quatalg.quaternion_algebra import QuaternionAlgebra from sage.categories.magmatic_algebras import MagmaticAlgebras @@ -13,38 +24,56 @@ 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.misc.table import table from sage.modules.free_module import FreeModule, VectorSpace -from sage.rings.integer_ring import ZZ -from sage.rings.number_field.number_field import NumberField, QuadraticField -from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing -from sage.rings.rational_field import QQ -from sage.rings.real_lazy import CLF, RLF - +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_operator import FiniteDimensionalEuclideanJordanAlgebraOperator from mjo.eja.eja_utils import _mat2vec class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): - # This is an ugly hack needed to prevent the category framework - # from implementing a coercion from our base ring (e.g. the - # rationals) into the algebra. First of all -- such a coercion is - # nonsense to begin with. But more importantly, it tries to do so - # in the category of rings, and since our algebras aren't - # associative they generally won't be rings. - _no_generic_basering_coercion = True + + 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, prefix='e', category=None, - natural_basis=None): + matrix_basis=None, + check_field=True, + check_axioms=True): """ SETUP:: - sage: from mjo.eja.eja_algebra import random_eja + sage: from mjo.eja.eja_algebra import ( + ....: FiniteDimensionalEuclideanJordanAlgebra, + ....: JordanSpinEJA, + ....: random_eja) EXAMPLES: @@ -56,9 +85,37 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): sage: x*y == y*x True + TESTS: + + The ``field`` we're given must be real with ``check_field=True``:: + + sage: JordanSpinEJA(2,QQbar) + Traceback (most recent call last): + ... + ValueError: scalar field is not real + + The multiplication table must be square with ``check_axioms=True``:: + + sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,((),())) + Traceback (most recent call last): + ... + ValueError: multiplication table is not square + """ - self._rank = rank - self._natural_basis = natural_basis + 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 table had better be square + n = len(mult_table) + if check_axioms: + if not all( len(l) == n for l in mult_table ): + raise ValueError("multiplication table is not square") + + self._matrix_basis = matrix_basis if category is None: category = MagmaticAlgebras(field).FiniteDimensional() @@ -66,7 +123,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): fda = super(FiniteDimensionalEuclideanJordanAlgebra, self) fda.__init__(field, - range(len(mult_table)), + range(n), prefix=prefix, category=category) self.print_options(bracket='') @@ -77,13 +134,27 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): # 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. - self._multiplication_table = [ map(lambda x: self.from_vector(x), ls) - for ls in mult_table ] - + self._multiplication_table = [ [ self.vector_space().zero() + for i in range(n) ] + for j in range(n) ] + # take advantage of symmetry + for i in range(n): + for j in range(i+1): + elt = self.from_vector(mult_table[i][j]) + self._multiplication_table[i][j] = elt + self._multiplication_table[j][i] = elt + + if check_axioms: + if not self._is_commutative(): + raise ValueError("algebra is not commutative") + 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): """ - Construct an element of this algebra from its natural + Construct an element of this algebra from its vector or matrix representation. This gets called only after the parent element _call_ method @@ -92,7 +163,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): SETUP:: sage: from mjo.eja.eja_algebra import (JordanSpinEJA, - ....: RealCartesianProductEJA, + ....: HadamardEJA, ....: RealSymmetricEJA) EXAMPLES: @@ -111,16 +182,16 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): sage: J(A) Traceback (most recent call last): ... - ArithmeticError: vector is not in free module + ValueError: not an element of this algebra TESTS: Ensure that we can convert any element of the two non-matrix - simple algebras (whose natural representations are their usual - vector representations) back and forth faithfully:: + simple algebras (whose matrix representations are columns) + back and forth faithfully:: sage: set_random_seed() - sage: J = RealCartesianProductEJA.random_instance() + sage: J = HadamardEJA.random_instance() sage: x = J.random_element() sage: J(x.to_vector().column()) == x True @@ -128,29 +199,37 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): sage: x = J.random_element() sage: J(x.to_vector().column()) == x True - """ + 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) - natural_basis = self.natural_basis() - basis_space = natural_basis[0].matrix_space() - if elt not in basis_space: - raise ValueError("not a naturally-represented algebra element") + 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 natural-basis coordinates + # 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. - V = VectorSpace(basis_space.base_ring(), elt.nrows()*elt.ncols()) - W = V.span_of_basis( _mat2vec(s) for s in natural_basis ) - coords = W.coordinate_vector(_mat2vec(elt)) - return self.from_vector(coords) + V = VectorSpace(self.base_ring(), elt.nrows()*elt.ncols()) + W = V.span_of_basis( _mat2vec(s) for s in self.matrix_basis() ) + try: + coords = W.coordinate_vector(_mat2vec(elt)) + except ArithmeticError: # vector is not in free module + raise ValueError(msg) + + return self.from_vector(coords) def _repr_(self): """ @@ -164,8 +243,8 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): Ensure that it says what we think it says:: - sage: JordanSpinEJA(2, field=QQ) - Euclidean Jordan algebra of dimension 2 over Rational Field + 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 @@ -176,165 +255,74 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): def product_on_basis(self, i, j): return self._multiplication_table[i][j] - def _a_regular_element(self): - """ - Guess a regular element. Needed to compute the basis for our - characteristic polynomial coefficients. - - SETUP:: - - sage: from mjo.eja.eja_algebra import random_eja - - TESTS: - - Ensure that this hacky method succeeds for every algebra that we - know how to construct:: - - sage: set_random_seed() - sage: J = random_eja() - sage: J._a_regular_element().is_regular() - True - - """ - 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 - - - @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). - """ - z = self._a_regular_element() - # Don't use the parent vector space directly here in case this - # happens to be a subalgebra. In that case, we would be e.g. - # two-dimensional but span_of_basis() would expect three - # coordinates. - V = VectorSpace(self.base_ring(), self.vector_space().dimension()) - basis = [ (z**k).to_vector() for k in range(self.rank()) ] - V1 = V.span_of_basis( basis ) - 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) - + def _is_commutative(self): + r""" + Whether or not this algebra's multiplication table 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. + """ + 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. + """ + 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. + """ + + # Used to check whether or not something is zero in an inexact + # ring. This number is sufficient to allow the construction of + # QuaternionHermitianEJA(2, RDF) with check_axioms=True. + epsilon = 1e-16 + + 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) + + if self.base_ring().is_exact(): + if diff != 0: + return False + else: + if diff.abs() > epsilon: + return False + + return True @cached_method - def _charpoly_matrix_system(self): - """ - 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. + def characteristic_polynomial_of(self): """ - r = self.rank() - n = self.dimension() - - # Turn my vector space into a module so that "vectors" can - # have multivatiate polynomial entries. - names = tuple('X' + str(i) for i in range(1,n+1)) - R = PolynomialRing(self.base_ring(), names) - - # Using change_ring() on the parent's vector space doesn't work - # here because, in a subalgebra, that vector space has a basis - # and change_ring() tries to bring the basis along with it. And - # that doesn't work unless the new ring is a PID, which it usually - # won't be. - V = FreeModule(R,n) - - # Now let x = (X1,X2,...,Xn) be the vector whose entries are - # indeterminates... - x = V(names) - - # And figure out the "left multiplication by x" matrix in - # that setting. - lmbx_cols = [] - monomial_matrices = [ self.monomial(i).operator().matrix() - for i in range(n) ] # don't recompute these! - for k in range(n): - ek = self.monomial(k).to_vector() - lmbx_cols.append( - sum( x[i]*(monomial_matrices[i]*ek) - for i in range(n) ) ) - Lx = matrix.column(R, lmbx_cols) - - # Now we can compute powers of x "symbolically" - x_powers = [self.one().to_vector(), x] - for d in range(2, r+1): - x_powers.append( Lx*(x_powers[-1]) ) - - idmat = matrix.identity(R, n) - - W = self._charpoly_basis_space() - W = W.change_ring(R.fraction_field()) - - # 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_powers = [ W.coordinate_vector(xp) for xp in x_powers ] - l2 = [idmat.column(k-1) for k in range(r+1, n+1)] - A_of_x = matrix.column(R, n, (x_powers[:r] + l2)) - return (A_of_x, x, x_powers[r], A_of_x.det()) - - - @cached_method - def characteristic_polynomial(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 @@ -346,7 +334,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): SETUP:: - sage: from mjo.eja.eja_algebra import JordanSpinEJA + sage: from mjo.eja.eja_algebra import JordanSpinEJA, TrivialEJA EXAMPLES: @@ -354,42 +342,64 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): 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().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 sum( a[k]*(t**k) for k in xrange(len(a)) ) + 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:: + + sage: J = HadamardEJA(2) + sage: J.coordinate_polynomial_ring() + Multivariate Polynomial Ring in X1, X2... + sage: J = RealSymmetricEJA(3,QQ) + 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): """ @@ -401,7 +411,9 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): SETUP:: - sage: from mjo.eja.eja_algebra import random_eja + sage: from mjo.eja.eja_algebra import (random_eja, + ....: HadamardEJA, + ....: BilinearFormEJA) EXAMPLES: @@ -414,10 +426,34 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): 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: (n == 1) or (x.inner_product(y) == (x*y).trace()/2) + True """ - X = x.natural_representation() - Y = y.natural_representation() - return self.natural_inner_product(X,Y) + B = self._inner_product_matrix + return (B*x.to_vector()).inner_product(y.to_vector()) def is_trivial(self): @@ -428,15 +464,19 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): SETUP:: - sage: from mjo.eja.eja_algebra import ComplexHermitianEJA + sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA, + ....: TrivialEJA) EXAMPLES:: sage: J = ComplexHermitianEJA(3) sage: J.is_trivial() False - sage: A = J.zero().subalgebra_generated_by() - sage: A.is_trivial() + + :: + + sage: J = TrivialEJA() + sage: J.is_trivial() True """ @@ -470,25 +510,40 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): """ M = list(self._multiplication_table) # copy - for i in xrange(len(M)): + for i in range(len(M)): # M had better be "square" M[i] = [self.monomial(i)] + M[i] M = [["*"] + list(self.gens())] + M return table(M, header_row=True, header_column=True, frame=True) - def natural_basis(self): + def matrix_basis(self): """ - Return a more-natural representation of this algebra's basis. + Return an (often more natural) representation of this algebras + basis as an ordered tuple of matrices. - 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.) + 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. - Note that this will always return a matrix. The standard basis - in `R^n` will be returned as `n`-by-`1` column matrices. + 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:: @@ -500,10 +555,10 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): sage: J = RealSymmetricEJA(2) sage: J.basis() Finite family {0: e0, 1: e1, 2: e2} - sage: J.natural_basis() + sage: J.matrix_basis() ( - [1 0] [ 0 1/2*sqrt2] [0 0] - [0 0], [1/2*sqrt2 0], [0 1] + [1 0] [ 0 0.7071067811865475?] [0 0] + [0 0], [0.7071067811865475? 0], [0 1] ) :: @@ -511,43 +566,38 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): sage: J = JordanSpinEJA(2) sage: J.basis() Finite family {0: e0, 1: e1} - sage: J.natural_basis() + sage: J.matrix_basis() ( [1] [0] [0], [1] ) - """ - if self._natural_basis is None: - M = self.natural_basis_space() + 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 natural_basis_space(self): - """ - Return the matrix space in which this algebra's natural basis - elements live. + def matrix_space(self): """ - if self._natural_basis is None or len(self._natural_basis) == 0: - return MatrixSpace(self.base_ring(), self.dimension(), 1) - else: - return self._natural_basis[0].matrix_space() - + 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). - @staticmethod - def natural_inner_product(X,Y): - """ - Compute the inner product of two naturally-represented elements. + 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. - For example in the real symmetric matrix EJA, this will compute - the trace inner-product of two n-by-n symmetric matrices. The - default should work for the real cartesian product EJA, the - Jordan spin EJA, and the real symmetric matrices. The others - will have to be overridden. + Matrix algebras override this with something more useful. """ - return (X.conjugate_transpose()*Y).trace() + 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() @cached_method @@ -557,12 +607,12 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): SETUP:: - sage: from mjo.eja.eja_algebra import (RealCartesianProductEJA, + sage: from mjo.eja.eja_algebra import (HadamardEJA, ....: random_eja) EXAMPLES:: - sage: J = RealCartesianProductEJA(5) + sage: J = HadamardEJA(5) sage: J.one() e0 + e1 + e2 + e3 + e4 @@ -585,6 +635,16 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): sage: actual == expected True + 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 + """ # We can brute-force compute the matrices of the operators # that correspond to the basis elements of this algebra. @@ -596,33 +656,224 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): # appeal to the "long vectors" isometry. oper_vecs = [ _mat2vec(g.operator().matrix()) for g in self.gens() ] - # Now we use basis linear algebra to find the coefficients, + # 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.column(self.base_ring(), oper_vecs) + A = matrix(self.base_ring(), oper_vecs) # 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()) ) - # Now if there's an identity element in the algebra, this should work. - coeffs = A.solve_right(b) - return self.linear_combination(zip(self.gens(), coeffs)) + # 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)) - def random_element(self): - # Temporary workaround for https://trac.sagemath.org/ticket/28327 - if self.is_trivial(): - return self.zero() - else: - s = super(FiniteDimensionalEuclideanJordanAlgebra, self) - return s.random_element() + def peirce_decomposition(self, c): + """ + 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. + + INPUT: + + - ``c`` -- an idempotent of this algebra. + + OUTPUT: + + A triple (J0, J5, J1) containing two subalgebras and one subspace + of this algebra, + + - ``J0`` -- the algebra on the eigenspace of ``c.operator()`` + corresponding to the eigenvalue zero. + + - ``J5`` -- the eigenspace (NOT a subalgebra) of ``c.operator()`` + corresponding to the eigenvalue one-half. + + - ``J1`` -- the algebra on the eigenspace of ``c.operator()`` + corresponding to the eigenvalue one. + + 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. + + SETUP:: + + sage: from mjo.eja.eja_algebra import random_eja, RealSymmetricEJA + + 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] + + TESTS: + + Every algebra decomposes trivially with respect to its identity + element:: + + 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 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:: + + 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 - def random_elements(self, count): + """ + 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) + + return (J0, J5, J1) + + + def random_element(self, thorough=False): + r""" + Return a random element of this algebra. + + 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. + + INPUT: + + - ``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 + + """ + # 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)) + + 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 + SETUP:: sage: from mjo.eja.eja_algebra import JordanSpinEJA @@ -637,23 +888,69 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): True """ - return tuple( self.random_element() for idx in xrange(count) ) + return tuple( self.random_element(thorough) + for idx in range(count) ) - def rank(self): + @cached_method + def _charpoly_coefficients(self): + r""" + The `r` polynomial coefficients of the "characteristic polynomial + of" function. """ - Return the rank of this EJA. + 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] - ALGORITHM: + @cached_method + def rank(self): + r""" + Return the rank of this EJA. - 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. + 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. SETUP:: - sage: from mjo.eja.eja_algebra import (JordanSpinEJA, + sage: from mjo.eja.eja_algebra import (HadamardEJA, + ....: JordanSpinEJA, ....: RealSymmetricEJA, ....: ComplexHermitianEJA, ....: QuaternionHermitianEJA, @@ -683,15 +980,29 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): TESTS: Ensure that every EJA that we know how to construct has a - positive integer rank:: + positive integer rank, unless the algebra is trivial in + which case its rank will be zero:: sage: set_random_seed() - sage: r = random_eja().rank() - sage: r in ZZ and r > 0 + sage: J = random_eja() + sage: r = J.rank() + sage: r in ZZ + True + sage: r > 0 or (r == 0 and J.is_trivial()) + True + + Ensure that computing the rank actually works, since the ranks + of all simple algebras are known and will be cached by default:: + + 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 """ - return self._rank + return len(self._charpoly_coefficients()) def vector_space(self): @@ -714,195 +1025,220 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): Element = FiniteDimensionalEuclideanJordanAlgebraElement +class RationalBasisEuclideanJordanAlgebraNg(FiniteDimensionalEuclideanJordanAlgebra): + def __init__(self, + field, + basis, + jordan_product, + inner_product, + orthonormalize=True, + prefix='e', + category=None, + check_field=True, + check_axioms=True): -class KnownRankEJA(object): - """ - A class for algebras that we actually know we can construct. The - main issue is that, for most of our methods to make sense, we need - to know the rank of our algebra. Thus we can't simply generate a - "random" algebra, or even check that a given basis and product - satisfy the axioms; because even if everything looks OK, we wouldn't - know the rank we need to actuallty build the thing. - - Not really a subclass of FDEJA because doing that causes method - resolution errors, e.g. + n = len(basis) + vector_basis = basis - TypeError: Error when calling the metaclass bases - Cannot create a consistent method resolution - order (MRO) for bases FiniteDimensionalEuclideanJordanAlgebra, - KnownRankEJA + from sage.matrix.matrix import is_Matrix + basis_is_matrices = False + degree = 0 + if n > 0: + if is_Matrix(basis[0]): + basis_is_matrices = True + vector_basis = tuple( map(_mat2vec,basis) ) + degree = basis[0].nrows()**2 + else: + degree = basis[0].degree() + + V = VectorSpace(field, degree) + + self._deorthonormalization_matrix = matrix.identity(field,n) + if orthonormalize: + A = matrix(field, vector_basis) + # uh oh, this is only the "usual" inner product + Q,R = A.gram_schmidt(orthonormal=True) + self._deorthonormalization_matrix = R.inverse().transpose() + vector_basis = Q.rows() + W = V.span_of_basis( vector_basis ) + if basis_is_matrices: + from mjo.eja.eja_utils import _vec2mat + basis = tuple( map(_vec2mat,vector_basis) ) + + mult_table = [ [0 for i in range(n)] for j in range(n) ] + ip_table = [ [0 for i in range(n)] for j in range(n) ] + + for i in range(n): + for j in range(i+1): + # do another mat2vec because the multiplication + # table is in terms of vectors + elt = _mat2vec(jordan_product(basis[i],basis[j])) + elt = W.coordinate_vector(elt) + mult_table[i][j] = elt + mult_table[j][i] = elt + ip = inner_product(basis[i],basis[j]) + ip_table[i][j] = ip + ip_table[j][i] = ip + + self._inner_product_matrix = matrix(field,ip_table) + + if basis_is_matrices: + for m in basis: + m.set_immutable() + + super().__init__(field, + mult_table, + prefix, + category, + basis, + check_field, + check_axioms) + +class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra): + r""" + Algebras whose basis consists of vectors with rational + entries. Equivalently, algebras whose multiplication tables + contain only rational coefficients. + + When an EJA has a basis that can be made rational, we can speed up + the computation of its characteristic polynomial by doing it over + ``QQ``. All of the named EJA constructors that we provide fall + into this category. """ - @staticmethod - def _max_test_case_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 quite vague -- 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. - - We default to five in this class, which is safe in `R^n`. The - matrix algebra subclasses (or any class where the "size" is - interpreted to be far less than the dimension) should override - with a smaller number. - """ - return 5 + @cached_method + def _charpoly_coefficients(self): + r""" + Override the parent method with something that tries to compute + over a faster (non-extension) field. - @classmethod - def random_instance(cls, field=QQ, **kwargs): - """ - Return a random instance of this type of algebra. + SETUP:: - Beware, this will crash for "most instances" because the - constructor below looks wrong. - """ - n = ZZ.random_element(cls._max_test_case_size()) + 1 - return cls(n, field, **kwargs) + sage: from mjo.eja.eja_algebra import JordanSpinEJA + EXAMPLES: -class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra, - KnownRankEJA): - """ - Return the Euclidean Jordan Algebra corresponding to the set - `R^n` under the Hadamard product. + 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):: - 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. + 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: + # There's no need to construct *another* algebra over the + # rationals if this one is already over the rationals. + superclass = super(RationalBasisEuclideanJordanAlgebra, self) + return superclass._charpoly_coefficients() + + mult_table = tuple( + tuple(map(lambda x: x.to_vector(), ls)) + for ls in self._multiplication_table + ) + + # Do the computation over the rationals. The answer will be + # the same, because our basis coordinates are (essentially) + # rational. + J = FiniteDimensionalEuclideanJordanAlgebra(QQ, + mult_table, + check_field=False, + check_axioms=False) + a = J._charpoly_coefficients() + return tuple(map(lambda x: x.change_ring(self.base_ring()), a)) + + +class ConcreteEuclideanJordanAlgebra: + 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 RealCartesianProductEJA - - EXAMPLES: - - This multiplication table can be verified by hand:: - - 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 + sage: from mjo.eja.eja_algebra import ConcreteEuclideanJordanAlgebra TESTS: - We can change the generator prefix:: + Our basis is normalized with respect to the algebra's inner + product, unless we specify otherwise:: - sage: RealCartesianProductEJA(3, prefix='r').gens() - (r0, r1, r2) + sage: set_random_seed() + sage: J = ConcreteEuclideanJordanAlgebra.random_instance() + sage: all( b.norm() == 1 for b in J.gens() ) + True - """ - def __init__(self, n, field=QQ, **kwargs): - V = VectorSpace(field, n) - mult_table = [ [ V.gen(i)*(i == j) for j in xrange(n) ] - for i in xrange(n) ] + 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:: - fdeja = super(RealCartesianProductEJA, self) - return fdeja.__init__(field, mult_table, rank=n, **kwargs) + sage: set_random_seed() + sage: J = ConcreteEuclideanJordanAlgebra.random_instance() + sage: x = J.random_element() + sage: x.operator().is_self_adjoint() + True + """ - def inner_product(self, x, y): + @staticmethod + def _max_random_instance_size(): """ - Faster to reimplement than to use natural representations. - - SETUP:: - - sage: from mjo.eja.eja_algebra import RealCartesianProductEJA - - TESTS: - - Ensure that this is the usual inner product for the algebras - over `R^n`:: - - sage: set_random_seed() - sage: J = RealCartesianProductEJA.random_instance() - sage: x,y = J.random_elements(2) - sage: X = x.natural_representation() - sage: Y = y.natural_representation() - sage: x.inner_product(y) == J.natural_inner_product(X,Y) - True + 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. """ - return x.to_vector().inner_product(y.to_vector()) - - -def random_eja(): - """ - Return a "random" finite-dimensional Euclidean Jordan Algebra. - - ALGORITHM: - - For now, we choose a random natural number ``n`` (greater than zero) - and then give you back one of the following: - - * The cartesian product of the rational numbers ``n`` times; this is - ``QQ^n`` with the Hadamard product. - - * The Jordan spin algebra on ``QQ^n``. - - * The ``n``-by-``n`` rational symmetric matrices with the symmetric - product. - - * The ``n``-by-``n`` complex-rational Hermitian matrices embedded - in the space of ``2n``-by-``2n`` real symmetric matrices. - - * 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. - - SETUP:: - - sage: from mjo.eja.eja_algebra import random_eja - - TESTS:: - - sage: random_eja() - Euclidean Jordan algebra of dimension... - - """ - classname = choice(KnownRankEJA.__subclasses__()) - return classname.random_instance() - - + raise NotImplementedError + @classmethod + def random_instance(cls, field=AA, **kwargs): + """ + Return a random instance of this type of algebra. + This method should be implemented in each subclass. + """ + from sage.misc.prandom import choice + eja_class = choice(cls.__subclasses__()) + return eja_class.random_instance(field) class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra): - @staticmethod - def _max_test_case_size(): - # Play it safe, since this will be squared and the underlying - # field can have dimension 4 (quaternions) too. - return 2 - def __init__(self, field, basis, rank, normalize_basis=True, **kwargs): + def __init__(self, field, basis, normalize_basis=True, **kwargs): """ Compared to the superclass constructor, we take a basis instead of a multiplication table because the latter can be computed in terms of the former when the product is known (like it is here). """ - # Used in this class's fast _charpoly_coeff() override. + # Used in this class's fast _charpoly_coefficients() override. self._basis_normalizers = None # We're going to loop through this a few times, so now's a good # time to ensure that it isn't a generator expression. basis = tuple(basis) - if rank > 1 and normalize_basis: + algebra_dim = len(basis) + degree = 0 # size of the matrices + if algebra_dim > 0: + degree = basis[0].nrows() + + if algebra_dim > 1 and normalize_basis: # We'll need sqrt(2) to normalize the basis, and this # winds up in the multiplication table, so the whole # algebra needs to be over the field extension. @@ -910,84 +1246,99 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra): z = R.gen() p = z**2 - 2 if p.is_irreducible(): - field = NumberField(p, 'sqrt2', embedding=RLF(2).sqrt()) + field = field.extension(p, 'sqrt2', embedding=RLF(2).sqrt()) basis = tuple( s.change_ring(field) for s in basis ) self._basis_normalizers = tuple( - ~(self.natural_inner_product(s,s).sqrt()) for s in basis ) - basis = tuple(s*c for (s,c) in izip(basis,self._basis_normalizers)) + ~(self.matrix_inner_product(s,s).sqrt()) for s in basis ) + basis = tuple(s*c for (s,c) in zip(basis,self._basis_normalizers)) - Qs = self.multiplication_table_from_matrix_basis(basis) + # Now compute the multiplication and inner product tables. + # We have to do this *after* normalizing the basis, because + # scaling affects the answers. + V = VectorSpace(field, degree**2) + W = V.span_of_basis( _mat2vec(s) for s in basis ) + mult_table = [[W.zero() for j in range(algebra_dim)] + for i in range(algebra_dim)] + ip_table = [[W.zero() for j in range(algebra_dim)] + for i in range(algebra_dim)] + for i in range(algebra_dim): + for j in range(algebra_dim): + mat_entry = (basis[i]*basis[j] + basis[j]*basis[i])/2 + mult_table[i][j] = W.coordinate_vector(_mat2vec(mat_entry)) - fdeja = super(MatrixEuclideanJordanAlgebra, self) - return fdeja.__init__(field, - Qs, - rank=rank, - natural_basis=basis, - **kwargs) + try: + # HACK: ignore the error here if we don't need the + # inner product (as is the case when we construct + # a dummy QQ-algebra for fast charpoly coefficients. + ip_table[i][j] = self.matrix_inner_product(basis[i], + basis[j]) + except: + pass + + try: + # HACK PART DEUX + self._inner_product_matrix = matrix(field,ip_table) + except: + pass + + super(MatrixEuclideanJordanAlgebra, self).__init__(field, + mult_table, + matrix_basis=basis, + **kwargs) + + if algebra_dim == 0: + self.one.set_cache(self.zero()) + else: + n = basis[0].nrows() + # The identity wrt (A,B) -> (AB + BA)/2 is independent of the + # details of this algebra. + self.one.set_cache(self(matrix.identity(field,n))) @cached_method - def _charpoly_coeff(self, i): - """ + def _charpoly_coefficients(self): + r""" Override the parent method with something that tries to compute over a faster (non-extension) field. """ - if self._basis_normalizers is None: - # We didn't normalize, so assume that the basis we started - # with had entries in a nice field. - return super(MatrixEuclideanJordanAlgebra, self)._charpoly_coeff(i) - else: - basis = ( (b/n) for (b,n) in izip(self.natural_basis(), - self._basis_normalizers) ) - - # Do this over the rationals and convert back at the end. - J = MatrixEuclideanJordanAlgebra(QQ, - basis, - self.rank(), - normalize_basis=False) - (_,x,_,_) = J._charpoly_matrix_system() - p = J._charpoly_coeff(i) - # p might be missing some vars, have to substitute "optionally" - pairs = izip(x.base_ring().gens(), self._basis_normalizers) - substitutions = { v: v*c for (v,c) in pairs } - result = p.subs(substitutions) - - # The result of "subs" can be either a coefficient-ring - # element or a polynomial. Gotta handle both cases. - if result in QQ: - return self.base_ring()(result) - else: - return result.change_ring(self.base_ring()) - - - @staticmethod - 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. - """ - # 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_of_basis( _mat2vec(s) for s in basis ) - n = len(basis) - mult_table = [[W.zero() for j in xrange(n)] for i in xrange(n)] - for i in xrange(n): - for j in xrange(n): - mat_entry = (basis[i]*basis[j] + basis[j]*basis[i])/2 - mult_table[i][j] = W.coordinate_vector(_mat2vec(mat_entry)) - - return mult_table + if self._basis_normalizers is None or self.base_ring() is QQ: + # We didn't normalize, or the basis we started with had + # entries in a nice field already. Just compute the thing. + return super(MatrixEuclideanJordanAlgebra, self)._charpoly_coefficients() + + basis = ( (b/n) for (b,n) in zip(self.matrix_basis(), + self._basis_normalizers) ) + + # Do this over the rationals and convert back at the end. + # Only works because we know the entries of the basis are + # integers. The argument ``check_axioms=False`` is required + # because the trace inner-product method for this + # class is a stub and can't actually be checked. + J = MatrixEuclideanJordanAlgebra(QQ, + basis, + normalize_basis=False, + check_field=False, + check_axioms=False) + a = J._charpoly_coefficients() + + # Unfortunately, changing the basis does change the + # coefficients of the characteristic polynomial, but since + # these are really the coefficients of the "characteristic + # polynomial of" function, everything is still nice and + # unevaluated. It's therefore "obvious" how scaling the + # basis affects the coordinate variables X1, X2, et + # cetera. Scaling the first basis vector up by "n" adds a + # factor of 1/n into every "X1" term, for example. So here + # we simply undo the basis_normalizer scaling that we + # performed earlier. + # + # The a[0] access here is safe because trivial algebras + # won't have any basis normalizers and therefore won't + # make it to this "else" branch. + XS = a[0].parent().gens() + subs_dict = { XS[i]: self._basis_normalizers[i]*XS[i] + for i in range(len(XS)) } + return tuple( a_i.subs(subs_dict) for a_i in a ) @staticmethod @@ -1014,22 +1365,17 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra): """ raise NotImplementedError - @classmethod - def natural_inner_product(cls,X,Y): + def matrix_inner_product(cls,X,Y): Xu = cls.real_unembed(X) Yu = cls.real_unembed(Y) tr = (Xu*Yu).trace() - if tr in RLF: - # It's real already. - return tr - # Otherwise, try the thing that works for complex numbers; and - # if that doesn't work, the thing that works for quaternions. try: - return tr.vector()[0] # real part, imag part is index 1 + # Works in QQ, AA, RDF, et cetera. + return tr.real() except AttributeError: - # A quaternions doesn't have a vector() method, but does + # 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] @@ -1053,7 +1399,8 @@ class RealMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): return M -class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra, KnownRankEJA): +class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra, + ConcreteEuclideanJordanAlgebra): """ The rank-n simple EJA consisting of real symmetric n-by-n matrices, the usual symmetric Jordan product, and the trace inner @@ -1074,12 +1421,20 @@ class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra, KnownRankEJA): sage: e2*e2 e2 + In theory, our "field" can be any subfield of the reals:: + + sage: RealSymmetricEJA(2, RDF) + Euclidean Jordan algebra of dimension 3 over Real Double Field + sage: RealSymmetricEJA(2, RR) + Euclidean Jordan algebra of dimension 3 over Real Field with + 53 bits of precision + TESTS: The dimension of this algebra is `(n^2 + n) / 2`:: sage: set_random_seed() - sage: n_max = RealSymmetricEJA._max_test_case_size() + sage: n_max = RealSymmetricEJA._max_random_instance_size() sage: n = ZZ.random_element(1, n_max) sage: J = RealSymmetricEJA(n) sage: J.dimension() == (n^2 + n)/2 @@ -1090,9 +1445,9 @@ class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra, KnownRankEJA): sage: set_random_seed() sage: J = RealSymmetricEJA.random_instance() sage: x,y = J.random_elements(2) - sage: actual = (x*y).natural_representation() - sage: X = x.natural_representation() - sage: Y = y.natural_representation() + 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 @@ -1104,24 +1459,10 @@ class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra, KnownRankEJA): sage: RealSymmetricEJA(3, prefix='q').gens() (q0, q1, q2, q3, q4, q5) - Our natural basis is normalized with respect to the natural inner - product unless we specify otherwise:: - - sage: set_random_seed() - sage: J = RealSymmetricEJA.random_instance() - sage: all( b.norm() == 1 for b in J.gens() ) - True - - Since our natural basis is normalized with respect to the natural - 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:: + We can construct the (trivial) algebra of rank zero:: - sage: set_random_seed() - sage: x = RealSymmetricEJA.random_instance().random_element() - sage: x.operator().matrix().is_symmetric() - True + sage: RealSymmetricEJA(0) + Euclidean Jordan algebra of dimension 0 over Algebraic Real Field """ @classmethod @@ -1145,8 +1486,8 @@ class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra, KnownRankEJA): # 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): + for i in range(n): + for j in range(i+1): Eij = matrix(field, n, lambda k,l: k==i and l==j) if i == j: Sij = Eij @@ -1157,13 +1498,24 @@ class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra, KnownRankEJA): @staticmethod - def _max_test_case_size(): + def _max_random_instance_size(): return 4 # Dimension 10 + @classmethod + def random_instance(cls, field=AA, **kwargs): + """ + Return a random instance of this type of algebra. + """ + n = ZZ.random_element(cls._max_random_instance_size() + 1) + return cls(n, field, **kwargs) - def __init__(self, n, field=QQ, **kwargs): + def __init__(self, n, field=AA, **kwargs): basis = self._denormalized_basis(n, field) - super(RealSymmetricEJA, self).__init__(field, basis, n, **kwargs) + super(RealSymmetricEJA, self).__init__(field, + basis, + check_axioms=False, + **kwargs) + self.rank.set_cache(n) class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): @@ -1181,7 +1533,7 @@ class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): EXAMPLES:: - sage: F = QuadraticField(-1, 'i') + sage: F = QuadraticField(-1, 'I') sage: x1 = F(4 - 2*i) sage: x2 = F(1 + 2*i) sage: x3 = F(-i) @@ -1199,9 +1551,8 @@ class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): Embedding is a homomorphism (isomorphism, in fact):: sage: set_random_seed() - sage: n_max = ComplexMatrixEuclideanJordanAlgebra._max_test_case_size() - sage: n = ZZ.random_element(n_max) - sage: F = QuadraticField(-1, 'i') + 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) @@ -1214,15 +1565,17 @@ class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): n = M.nrows() if M.ncols() != n: raise ValueError("the matrix 'M' must be square") - field = M.base_ring() + + # We don't need any adjoined elements... + field = M.base_ring().base_ring() + blocks = [] for z in M.list(): - a = z.vector()[0] # real part, I guess - b = z.vector()[1] # imag part, I guess + 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]])) - # We can drop the imaginaries here. - return matrix.block(field.base_ring(), n, blocks) + return matrix.block(field, n, blocks) @staticmethod @@ -1242,15 +1595,15 @@ class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): ....: [ 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] + [ 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: F = QuadraticField(-1, 'I') sage: M = random_matrix(F, 3) sage: Me = ComplexMatrixEuclideanJordanAlgebra.real_embed(M) sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(Me) == M @@ -1263,17 +1616,24 @@ class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): if not n.mod(2).is_zero(): raise ValueError("the matrix 'M' must be a complex embedding") - field = QQ + # 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() - F = NumberField(z**2 + 1,'i', embedding=CLF(-1).sqrt()) + 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 xrange(n/2): - for j in xrange(n/2): + 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') @@ -1286,9 +1646,9 @@ class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): @classmethod - def natural_inner_product(cls,X,Y): + def matrix_inner_product(cls,X,Y): """ - Compute a natural inner product in this algebra directly from + Compute a matrix inner product in this algebra directly from its real embedding. SETUP:: @@ -1303,20 +1663,21 @@ class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): sage: set_random_seed() sage: J = ComplexHermitianEJA.random_instance() sage: x,y = J.random_elements(2) - sage: Xe = x.natural_representation() - sage: Ye = y.natural_representation() + 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().vector()[0] - sage: actual = ComplexHermitianEJA.natural_inner_product(Xe,Ye) + sage: expected = (X*Y).trace().real() + sage: actual = ComplexHermitianEJA.matrix_inner_product(Xe,Ye) sage: actual == expected True """ - return RealMatrixEuclideanJordanAlgebra.natural_inner_product(X,Y)/2 + return RealMatrixEuclideanJordanAlgebra.matrix_inner_product(X,Y)/2 -class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra, KnownRankEJA): +class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra, + ConcreteEuclideanJordanAlgebra): """ The rank-n simple EJA consisting of complex Hermitian n-by-n matrices over the real numbers, the usual symmetric Jordan product, @@ -1327,12 +1688,22 @@ class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra, KnownRankEJA): sage: from mjo.eja.eja_algebra import ComplexHermitianEJA + EXAMPLES: + + In theory, our "field" can be any subfield of the reals:: + + sage: ComplexHermitianEJA(2, RDF) + Euclidean Jordan algebra of dimension 4 over Real Double Field + sage: ComplexHermitianEJA(2, RR) + Euclidean Jordan algebra of dimension 4 over Real Field with + 53 bits of precision + TESTS: The dimension of this algebra is `n^2`:: sage: set_random_seed() - sage: n_max = ComplexHermitianEJA._max_test_case_size() + sage: n_max = ComplexHermitianEJA._max_random_instance_size() sage: n = ZZ.random_element(1, n_max) sage: J = ComplexHermitianEJA(n) sage: J.dimension() == n^2 @@ -1343,9 +1714,9 @@ class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra, KnownRankEJA): sage: set_random_seed() sage: J = ComplexHermitianEJA.random_instance() sage: x,y = J.random_elements(2) - sage: actual = (x*y).natural_representation() - sage: X = x.natural_representation() - sage: Y = y.natural_representation() + 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 @@ -1357,24 +1728,10 @@ class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra, KnownRankEJA): sage: ComplexHermitianEJA(2, prefix='z').gens() (z0, z1, z2, z3) - Our natural basis is normalized with respect to the natural inner - product unless we specify otherwise:: - - sage: set_random_seed() - sage: J = ComplexHermitianEJA.random_instance() - sage: all( b.norm() == 1 for b in J.gens() ) - True - - Since our natural basis is normalized with respect to the natural - 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:: + We can construct the (trivial) algebra of rank zero:: - sage: set_random_seed() - sage: x = ComplexHermitianEJA.random_instance().random_element() - sage: x.operator().matrix().is_symmetric() - True + sage: ComplexHermitianEJA(0) + Euclidean Jordan algebra of dimension 0 over Algebraic Real Field """ @@ -1403,9 +1760,9 @@ class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra, KnownRankEJA): True """ - R = PolynomialRing(QQ, 'z') + R = PolynomialRing(field, 'z') z = R.gen() - F = NumberField(z**2 + 1, 'I', embedding=CLF(-1).sqrt()) + F = field.extension(z**2 + 1, 'I') I = F.gen() # This is like the symmetric case, but we need to be careful: @@ -1414,8 +1771,8 @@ class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra, KnownRankEJA): # * The diagonal will (as a result) be real. # S = [] - for i in xrange(n): - for j in xrange(i+1): + 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) @@ -1432,10 +1789,25 @@ class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra, KnownRankEJA): return ( s.change_ring(field) for s in S ) - def __init__(self, n, field=QQ, **kwargs): + def __init__(self, n, field=AA, **kwargs): basis = self._denormalized_basis(n,field) - super(ComplexHermitianEJA,self).__init__(field, basis, n, **kwargs) + super(ComplexHermitianEJA,self).__init__(field, + basis, + check_axioms=False, + **kwargs) + self.rank.set_cache(n) + @staticmethod + def _max_random_instance_size(): + return 3 # Dimension 9 + + @classmethod + def random_instance(cls, field=AA, **kwargs): + """ + Return a random instance of this type of algebra. + """ + n = ZZ.random_element(cls._max_random_instance_size() + 1) + return cls(n, field, **kwargs) class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): @staticmethod @@ -1467,8 +1839,7 @@ class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): Embedding is a homomorphism (isomorphism, in fact):: sage: set_random_seed() - sage: n_max = QuaternionMatrixEuclideanJordanAlgebra._max_test_case_size() - sage: n = ZZ.random_element(n_max) + sage: n = ZZ.random_element(2) sage: Q = QuaternionAlgebra(QQ,-1,-1) sage: X = random_matrix(Q, n) sage: Y = random_matrix(Q, n) @@ -1484,7 +1855,7 @@ class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): if M.ncols() != n: raise ValueError("the matrix 'M' must be square") - F = QuadraticField(-1, 'i') + F = QuadraticField(-1, 'I') i = F.gen() blocks = [] @@ -1552,27 +1923,27 @@ class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): # 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): + 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].vector()[0] # real part - z += submat[0,0].vector()[1]*i # imag part - z += submat[0,1].vector()[0]*j # real part - z += submat[0,1].vector()[1]*k # imag part + 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 natural_inner_product(cls,X,Y): + def matrix_inner_product(cls,X,Y): """ - Compute a natural inner product in this algebra directly from + Compute a matrix inner product in this algebra directly from its real embedding. SETUP:: @@ -1587,22 +1958,22 @@ class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): sage: set_random_seed() sage: J = QuaternionHermitianEJA.random_instance() sage: x,y = J.random_elements(2) - sage: Xe = x.natural_representation() - sage: Ye = y.natural_representation() + 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.natural_inner_product(Xe,Ye) + sage: actual = QuaternionHermitianEJA.matrix_inner_product(Xe,Ye) sage: actual == expected True """ - return RealMatrixEuclideanJordanAlgebra.natural_inner_product(X,Y)/4 + return RealMatrixEuclideanJordanAlgebra.matrix_inner_product(X,Y)/4 class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra, - KnownRankEJA): - """ + ConcreteEuclideanJordanAlgebra): + 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 @@ -1612,12 +1983,22 @@ class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra, sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA + EXAMPLES: + + In theory, our "field" can be any subfield of the reals:: + + sage: QuaternionHermitianEJA(2, RDF) + Euclidean Jordan algebra of dimension 6 over Real Double Field + sage: QuaternionHermitianEJA(2, RR) + Euclidean Jordan algebra of dimension 6 over Real Field with + 53 bits of precision + TESTS: The dimension of this algebra is `2*n^2 - n`:: sage: set_random_seed() - sage: n_max = QuaternionHermitianEJA._max_test_case_size() + sage: n_max = QuaternionHermitianEJA._max_random_instance_size() sage: n = ZZ.random_element(1, n_max) sage: J = QuaternionHermitianEJA(n) sage: J.dimension() == 2*(n^2) - n @@ -1628,9 +2009,9 @@ class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra, sage: set_random_seed() sage: J = QuaternionHermitianEJA.random_instance() sage: x,y = J.random_elements(2) - sage: actual = (x*y).natural_representation() - sage: X = x.natural_representation() - sage: Y = y.natural_representation() + 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 @@ -1642,24 +2023,10 @@ class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra, sage: QuaternionHermitianEJA(2, prefix='a').gens() (a0, a1, a2, a3, a4, a5) - Our natural basis is normalized with respect to the natural inner - product unless we specify otherwise:: - - sage: set_random_seed() - sage: J = QuaternionHermitianEJA.random_instance() - sage: all( b.norm() == 1 for b in J.gens() ) - True - - Since our natural basis is normalized with respect to the natural - 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:: + We can construct the (trivial) algebra of rank zero:: - sage: set_random_seed() - sage: x = QuaternionHermitianEJA.random_instance().random_element() - sage: x.operator().matrix().is_symmetric() - True + sage: QuaternionHermitianEJA(0) + Euclidean Jordan algebra of dimension 0 over Algebraic Real Field """ @classmethod @@ -1695,8 +2062,8 @@ class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra, # * The diagonal will (as a result) be real. # S = [] - for i in xrange(n): - for j in xrange(i+1): + 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) @@ -1718,16 +2085,276 @@ class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra, return ( s.change_ring(field) for s in S ) - def __init__(self, n, field=QQ, **kwargs): + def __init__(self, n, field=AA, **kwargs): basis = self._denormalized_basis(n,field) - super(QuaternionHermitianEJA,self).__init__(field, basis, n, **kwargs) + super(QuaternionHermitianEJA,self).__init__(field, + basis, + check_axioms=False, + **kwargs) + self.rank.set_cache(n) + + @staticmethod + def _max_random_instance_size(): + r""" + The maximum rank of a random QuaternionHermitianEJA. + """ + return 2 # Dimension 6 + + @classmethod + def random_instance(cls, field=AA, **kwargs): + """ + Return a random instance of this type of algebra. + """ + n = ZZ.random_element(cls._max_random_instance_size() + 1) + return cls(n, field, **kwargs) + + +class HadamardEJA(RationalBasisEuclideanJordanAlgebra, + 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, field=AA, **kwargs): + V = VectorSpace(field, n) + mult_table = [ [ V.gen(i)*(i == j) for j in range(n) ] + for i in range(n) ] + + # Inner products are real numbers and not algebra + # elements, so once we turn the algebra element + # into a vector in inner_product(), we never go + # back. As a result -- contrary to what we do with + # self._multiplication_table -- we store the inner + # product table as a plain old matrix and not as + # an algebra operator. + ip_table = matrix.identity(field,n) + self._inner_product_matrix = ip_table + + super(HadamardEJA, self).__init__(field, + mult_table, + 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, field=AA, **kwargs): + """ + Return a random instance of this type of algebra. + """ + n = ZZ.random_element(cls._max_random_instance_size() + 1) + return cls(n, field, **kwargs) + + +class BilinearFormEJA(RationalBasisEuclideanJordanAlgebra, + 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:: + + 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) + sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis() + sage: V = J.vector_space() + sage: sis = [ J.from_vector(V([0] + (M.inverse()*ei).list())) + ....: 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, field=AA, **kwargs): + n = B.nrows() + + if not B.is_positive_definite(): + raise ValueError("bilinear form is not positive-definite") + + V = VectorSpace(field, n) + mult_table = [[V.zero() for j in range(n)] for i in range(n)] + for i in range(n): + for j in range(n): + x = V.gen(i) + y = V.gen(j) + x0 = x[0] + xbar = x[1:] + y0 = y[0] + ybar = y[1:] + z0 = (B*x).inner_product(y) + zbar = y0*xbar + x0*ybar + z = V([z0] + zbar.list()) + mult_table[i][j] = z + + # Inner products are real numbers and not algebra + # elements, so once we turn the algebra element + # into a vector in inner_product(), we never go + # back. As a result -- contrary to what we do with + # self._multiplication_table -- we store the inner + # product table as a plain old matrix and not as + # an algebra operator. + ip_table = B + self._inner_product_matrix = ip_table + + super(BilinearFormEJA, self).__init__(field, + mult_table, + check_axioms=False, + **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, field=AA, **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(field, n) + return cls(B, field, **kwargs) + + B11 = matrix.identity(field,1) + M = matrix.random(field, n-1) + I = matrix.identity(field, n-1) + alpha = field.zero() + while alpha.is_zero(): + alpha = field.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, field, **kwargs) -class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra, KnownRankEJA): + +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:: @@ -1760,50 +2387,324 @@ class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra, KnownRankEJA): 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, field=AA, **kwargs): + # This is a special case of the BilinearFormEJA with the identity + # matrix as its bilinear form. + B = matrix.identity(field, n) + super(JordanSpinEJA, self).__init__(B, field, **kwargs) + + @staticmethod + def _max_random_instance_size(): + r""" + The maximum dimension of a random JordanSpinEJA. + """ + return 5 + + @classmethod + def random_instance(cls, field=AA, **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, field, **kwargs) + + +class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra, + 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, field=AA, **kwargs): + mult_table = [] + self._inner_product_matrix = matrix(field,0) + super(TrivialEJA, self).__init__(field, + mult_table, + check_axioms=False, + **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, field=AA, **kwargs): + # We don't take a "size" argument so the superclass method is + # inappropriate for us. + return cls(field, **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(AA) + sage: J2 = random_eja(QQ) + sage: J = DirectSumEJA(J1,J2) + Traceback (most recent call last): + ... + ValueError: algebras must share the same base field + """ - def __init__(self, n, field=QQ, **kwargs): + 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 xrange(n)] for i in xrange(n)] - for i in xrange(n): - for j in xrange(n): - x = V.gen(i) - y = V.gen(j) - x0 = x[0] - xbar = x[1:] - y0 = y[0] - ybar = y[1:] - # z = x*y - z0 = x.inner_product(y) - zbar = y0*xbar + x0*ybar - z = V([z0] + zbar.list()) - mult_table[i][j] = z + mult_table = [ [ V.zero() for j in range(n) ] + for i in range(n) ] + for i in range(n1): + for j in range(n1): + p = (J1.monomial(i)*J1.monomial(j)).to_vector() + mult_table[i][j] = V(p.list() + [field.zero()]*n2) - # 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, self) - return fdeja.__init__(field, mult_table, rank=min(n,2), **kwargs) + for i in range(n2): + for j in range(n2): + p = (J2.monomial(i)*J2.monomial(j)).to_vector() + mult_table[n1+i][n1+j] = V([field.zero()]*n1 + p.list()) + + super(DirectSumEJA, self).__init__(field, + mult_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,QQ) + sage: J2 = JordanSpinEJA(3,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) - def inner_product(self, x, y): """ - Faster to reimplement than to use natural representations. + 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 + 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: - Ensure that this is the usual inner product for the algebras - over `R^n`:: + 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: J = JordanSpinEJA.random_instance() - sage: x,y = J.random_elements(2) - sage: X = x.natural_representation() - sage: Y = y.natural_representation() - sage: x.inner_product(y) == J.natural_inner_product(X,Y) + 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 """ - return x.to_vector().inner_product(y.to_vector()) + (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): + 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,QQ) + sage: J2 = QuaternionHermitianEJA(2,QQ,normalize_basis=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