X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=mjo%2Feja%2Feja_algebra.py;h=c35bf256453cfa9a3758c034a47a043745eb0c2b;hb=7dc61e259969f86d4f31706fcfe411e1677389f2;hp=081abd303c32d35c9376ad99ecf2732a14ea3e9b;hpb=95ae8e7b0ddca840da9631603a2f37cca888468b;p=sage.d.git diff --git a/mjo/eja/eja_algebra.py b/mjo/eja/eja_algebra.py index 081abd3..c35bf25 100644 --- a/mjo/eja/eja_algebra.py +++ b/mjo/eja/eja_algebra.py @@ -133,7 +133,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): deortho_vector_basis = tuple( V(b.list()) for b in basis ) from mjo.eja.eja_utils import gram_schmidt - basis = gram_schmidt(basis, inner_product) + basis = tuple(gram_schmidt(basis, inner_product)) # Save the (possibly orthonormalized) matrix basis for # later... @@ -158,7 +158,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): # Now we actually compute the multiplication and inner-product # tables/matrices using the possibly-orthonormalized basis. - self._inner_product_matrix = matrix.zero(field, n) + self._inner_product_matrix = matrix.identity(field, n) self._multiplication_table = [ [0 for j in range(i+1)] for i in range(n) ] @@ -171,15 +171,20 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): q_i = basis[i] q_j = basis[j] - elt = jordan_product(q_i, q_j) - ip = inner_product(q_i, q_j) - # The jordan product returns a matrixy answer, so we # have to convert it to the algebra coordinates. + elt = jordan_product(q_i, q_j) elt = W.coordinate_vector(V(elt.list())) self._multiplication_table[i][j] = self.from_vector(elt) - self._inner_product_matrix[i,j] = ip - self._inner_product_matrix[j,i] = ip + + if not orthonormalize: + # If we're orthonormalizing the basis with respect + # to an inner-product, then the inner-product + # matrix with respect to the resulting basis is + # just going to be the identity. + ip = inner_product(q_i, q_j) + self._inner_product_matrix[i,j] = ip + self._inner_product_matrix[j,i] = ip self._inner_product_matrix._cache = {'hermitian': True} self._inner_product_matrix.set_immutable() @@ -315,7 +320,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): This method should of course always return ``True``, unless this algebra was constructed with ``check_axioms=False`` and - passed an invalid multiplication table. + passed an invalid Jordan or inner-product. """ # Used to check whether or not something is zero in an inexact @@ -732,7 +737,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): if self.is_trivial(): return MatrixSpace(self.base_ring(), 0) else: - return self._matrix_basis[0].matrix_space() + return self.matrix_basis()[0].parent() @cached_method @@ -745,23 +750,57 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): sage: from mjo.eja.eja_algebra import (HadamardEJA, ....: random_eja) - EXAMPLES:: + EXAMPLES: + + We can compute unit element in the Hadamard EJA:: sage: J = HadamardEJA(5) sage: J.one() e0 + e1 + e2 + e3 + e4 + The unit element in the Hadamard EJA is inherited in the + subalgebras generated by its elements:: + + sage: J = HadamardEJA(5) + sage: J.one() + e0 + e1 + e2 + e3 + e4 + sage: x = sum(J.gens()) + sage: A = x.subalgebra_generated_by(orthonormalize=False) + sage: A.one() + f0 + sage: A.one().superalgebra_element() + e0 + e1 + e2 + e3 + e4 + TESTS: - The identity element acts like the identity:: + The identity element acts like the identity, regardless of + whether or not we orthonormalize:: sage: set_random_seed() sage: J = random_eja() sage: x = J.random_element() sage: J.one()*x == x and x*J.one() == x True + sage: A = x.subalgebra_generated_by() + sage: y = A.random_element() + sage: A.one()*y == y and y*A.one() == y + True + + :: + + sage: set_random_seed() + sage: J = random_eja(field=QQ, orthonormalize=False) + sage: x = J.random_element() + sage: J.one()*x == x and x*J.one() == x + True + sage: A = x.subalgebra_generated_by(orthonormalize=False) + sage: y = A.random_element() + sage: A.one()*y == y and y*A.one() == y + True - The matrix of the unit element's operator is the identity:: + The matrix of the unit element's operator is the identity, + regardless of the base field and whether or not we + orthonormalize:: sage: set_random_seed() sage: J = random_eja() @@ -769,6 +808,27 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): sage: expected = matrix.identity(J.base_ring(), J.dimension()) sage: actual == expected True + sage: x = J.random_element() + sage: A = x.subalgebra_generated_by() + sage: actual = A.one().operator().matrix() + sage: expected = matrix.identity(A.base_ring(), A.dimension()) + sage: actual == expected + True + + :: + + sage: set_random_seed() + sage: J = random_eja(field=QQ, orthonormalize=False) + sage: actual = J.one().operator().matrix() + sage: expected = matrix.identity(J.base_ring(), J.dimension()) + sage: actual == expected + True + sage: x = J.random_element() + sage: A = x.subalgebra_generated_by(orthonormalize=False) + sage: actual = A.one().operator().matrix() + sage: expected = matrix.identity(A.base_ring(), A.dimension()) + sage: actual == expected + True Ensure that the cached unit element (often precomputed by hand) agrees with the computed one:: @@ -780,6 +840,15 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): sage: J.one() == cached True + :: + + sage: set_random_seed() + sage: J = random_eja(field=QQ, orthonormalize=False) + 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. @@ -1032,6 +1101,21 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): r""" The `r` polynomial coefficients of the "characteristic polynomial of" function. + + SETUP:: + + sage: from mjo.eja.eja_algebra import random_eja + + TESTS: + + The theory shows that these are all homogeneous polynomials of + a known degree:: + + sage: set_random_seed() + sage: J = random_eja() + sage: all(p.is_homogeneous() for p in J._charpoly_coefficients()) + True + """ n = self.dimension() R = self.coordinate_polynomial_ring() @@ -1067,10 +1151,17 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): # 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] + # 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. We don't bother to trim A_rref + # down to a square matrix and solve the resulting system, + # because the upper-left r-by-r portion of A_rref is + # guaranteed to be the identity matrix, so e.g. + # + # A_rref.solve_right(Y) + # + # would just be returning Y. + return (-E*b)[:r].change_ring(R) @cached_method def rank(self): @@ -1131,7 +1222,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): sage: set_random_seed() # long time sage: J = random_eja() # long time - sage: caches = J.rank() # long time + sage: cached = J.rank() # long time sage: J.rank.clear_cache() # long time sage: J.rank() == cached # long time True @@ -1196,6 +1287,7 @@ class RationalBasisEJA(FiniteDimensionalEJA): if not all( all(b_i in QQ for b_i in b.list()) for b in basis ): raise TypeError("basis not rational") + self._rational_algebra = None if field is not QQ: # There's no point in constructing the extra algebra if this # one is already rational. @@ -1256,7 +1348,14 @@ class RationalBasisEJA(FiniteDimensionalEJA): 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 + if self._deortho_matrix is None: + # This can happen if our base ring was, say, AA and we + # chose not to (or didn't need to) orthonormalize. It's + # still faster to do the computations over QQ even if + # the numbers in the boxes stay the same. + return tuple(a) + + # Otherwise, convert the coordinate variables back to the # deorthonormalized ones. R = self.coordinate_polynomial_ring() from sage.modules.free_module_element import vector @@ -1590,6 +1689,38 @@ class RealSymmetricEJA(ConcreteEJA, RealMatrixEJA): class ComplexMatrixEJA(MatrixEJA): + # A manual dictionary-cache for the complex_extension() method, + # since apparently @classmethods can't also be @cached_methods. + _complex_extension = {} + + @classmethod + def complex_extension(cls,field): + r""" + The complex field that we embed/unembed, as an extension + of the given ``field``. + """ + if field in cls._complex_extension: + return cls._complex_extension[field] + + # Sage doesn't know how to adjoin the complex "i" (the root of + # x^2 + 1) to a field in a general way. Here, we just enumerate + # all of the cases that I have cared to support so far. + 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 + elif not field.is_exact(): + # RDF or RR + F = field.complex_field() + else: + # Works for QQ and... maybe some other fields. + R = PolynomialRing(field, 'z') + z = R.gen() + F = field.extension(z**2 + 1, 'I', embedding=CLF(-1).sqrt()) + + cls._complex_extension[field] = F + return F + @staticmethod def dimension_over_reals(): return 2 @@ -1644,9 +1775,10 @@ class ComplexMatrixEJA(MatrixEJA): 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]])) + a = z.real() + b = z.imag() + blocks.append(matrix(field, 2, [ [ a, b], + [-b, a] ])) return matrix.block(field, n, blocks) @@ -1685,26 +1817,7 @@ class ComplexMatrixEJA(MatrixEJA): super(ComplexMatrixEJA,cls).real_unembed(M) n = ZZ(M.nrows()) d = cls.dimension_over_reals() - - # 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() - - # Sage doesn't know how to adjoin the complex "i" (the root of - # x^2 + 1) to a field in a general way. Here, we just enumerate - # all of the cases that I have cared to support so far. - 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 - elif not field.is_exact(): - # RDF or RR - F = field.complex_field() - else: - # Works for QQ and... maybe some other fields. - F = field.extension(z**2 + 1, 'I', embedding=CLF(-1).sqrt()) + F = cls.complex_extension(M.base_ring()) i = F.gen() # Go top-left to bottom-right (reading order), converting every @@ -1800,7 +1913,6 @@ class ComplexHermitianEJA(ConcreteEJA, ComplexMatrixEJA): 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 @@ -1818,18 +1930,27 @@ class ComplexHermitianEJA(ConcreteEJA, ComplexMatrixEJA): # * The diagonal will (as a result) be real. # S = [] + Eij = matrix.zero(F,n) for i in range(n): for j in range(i+1): - Eij = matrix(F, n, lambda k,l: k==i and l==j) + # "build" E_ij + Eij[i,j] = 1 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()) + Eij[j,i] = 1 # Eij = Eij + Eij.transpose() + Sij_real = cls.real_embed(Eij) S.append(Sij_real) - Sij_imag = cls.real_embed(I*Eij - I*Eij.transpose()) + # Eij = I*Eij - I*Eij.transpose() + Eij[i,j] = I + Eij[j,i] = -I + Sij_imag = cls.real_embed(Eij) S.append(Sij_imag) + Eij[j,i] = 0 + # "erase" E_ij + Eij[i,j] = 0 # Since we embedded these, we can drop back to the "field" that we # started with instead of the complex extension "F". @@ -1865,6 +1986,25 @@ class ComplexHermitianEJA(ConcreteEJA, ComplexMatrixEJA): return cls(n, **kwargs) class QuaternionMatrixEJA(MatrixEJA): + + # A manual dictionary-cache for the quaternion_extension() method, + # since apparently @classmethods can't also be @cached_methods. + _quaternion_extension = {} + + @classmethod + def quaternion_extension(cls,field): + r""" + The quaternion field that we embed/unembed, as an extension + of the given ``field``. + """ + if field in cls._quaternion_extension: + return cls._quaternion_extension[field] + + Q = QuaternionAlgebra(field,-1,-1) + + cls._quaternion_extension[field] = Q + return Q + @staticmethod def dimension_over_reals(): return 4 @@ -1969,8 +2109,7 @@ class QuaternionMatrixEJA(MatrixEJA): # 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) + Q = cls.quaternion_extension(M.base_ring()) i,j,k = Q.gens() # Go top-left to bottom-right (reading order), converting every @@ -2085,23 +2224,39 @@ class QuaternionHermitianEJA(ConcreteEJA, QuaternionMatrixEJA): # * The diagonal will (as a result) be real. # S = [] + Eij = matrix.zero(Q,n) for i in range(n): for j in range(i+1): - Eij = matrix(Q, n, lambda k,l: k==i and l==j) + # "build" E_ij + Eij[i,j] = 1 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()) + # Eij = Eij + Eij.transpose() + Eij[j,i] = 1 + Sij_real = cls.real_embed(Eij) S.append(Sij_real) - Sij_I = cls.real_embed(I*Eij - I*Eij.transpose()) + # Eij = I*(Eij - Eij.transpose()) + Eij[i,j] = I + Eij[j,i] = -I + Sij_I = cls.real_embed(Eij) S.append(Sij_I) - Sij_J = cls.real_embed(J*Eij - J*Eij.transpose()) + # Eij = J*(Eij - Eij.transpose()) + Eij[i,j] = J + Eij[j,i] = -J + Sij_J = cls.real_embed(Eij) S.append(Sij_J) - Sij_K = cls.real_embed(K*Eij - K*Eij.transpose()) + # Eij = K*(Eij - Eij.transpose()) + Eij[i,j] = K + Eij[j,i] = -K + Sij_K = cls.real_embed(Eij) S.append(Sij_K) + Eij[j,i] = 0 + # "erase" E_ij + Eij[i,j] = 0 # Since we embedded these, we can drop back to the "field" that we # started with instead of the quaternion algebra "Q". @@ -2307,10 +2462,21 @@ class BilinearFormEJA(ConcreteEJA): ....: 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") + # The matrix "B" is supplied by the user in most cases, + # so it makes sense to check whether or not its positive- + # definite unless we are specifically asked not to... + if ("check_axioms" not in kwargs) or kwargs["check_axioms"]: + if not B.is_positive_definite(): + raise ValueError("bilinear form is not positive-definite") + + # However, all of the other data for this EJA is computed + # by us in manner that guarantees the axioms are + # satisfied. So, again, unless we are specifically asked to + # verify things, we'll skip the rest of the checks. + if "check_axioms" not in kwargs: kwargs["check_axioms"] = False def inner_product(x,y): return (y.T*B*x)[0,0] @@ -2321,13 +2487,9 @@ class BilinearFormEJA(ConcreteEJA): xbar = x[1:,0] y0 = y[0,0] ybar = y[1:,0] - z0 = (y.T*x)[0,0] + z0 = inner_product(y,x) zbar = y0*xbar + x0*ybar - return P([0] + zbar.list()) - - # We know this is a valid EJA, but will double-check - # if the user passes check_axioms=True. - if "check_axioms" not in kwargs: kwargs["check_axioms"] = False + return P([z0] + zbar.list()) n = B.nrows() column_basis = tuple( b.column() for b in FreeModule(ZZ, n).basis() ) @@ -2512,100 +2674,119 @@ class TrivialEJA(ConcreteEJA): # inappropriate for us. return cls(**kwargs) -# class DirectSumEJA(ConcreteEJA): -# 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:: +class DirectSumEJA(FiniteDimensionalEJA): + 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. -# sage: from mjo.eja.eja_algebra import (HadamardEJA, -# ....: JordanSpinEJA, -# ....: DirectSumEJA) + SETUP:: -# EXAMPLES:: + sage: from mjo.eja.eja_algebra import (random_eja, + ....: HadamardEJA, + ....: RealSymmetricEJA, + ....: DirectSumEJA) -# 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) + EXAMPLES:: + + sage: J1 = HadamardEJA(2) + sage: J2 = RealSymmetricEJA(3) + sage: J = DirectSumEJA(J1,J2) + sage: J.dimension() + 8 + sage: J.rank() + 5 + sage: J.matrix_space() + The Cartesian product of (Full MatrixSpace of 2 by 1 dense matrices + over Algebraic Real Field, Full MatrixSpace of 3 by 3 dense matrices + over Algebraic Real Field) + + 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() + + M = J1.matrix_space().cartesian_product(J2.matrix_space()) + self._cartprod_algebra = J1.cartesian_product(J2) + + self._matrix_basis = tuple( [M((a,0)) for a in J1.matrix_basis()] + + [M((0,b)) for b in J2.matrix_basis()] ) + + n = len(self._matrix_basis) + self._sets = None + CombinatorialFreeModule.__init__( + self, + field, + range(n), + category=self._cartprod_algebra.category(), + bracket=False, + **kwargs) + self.rank.set_cache(J1.rank() + J2.rank()) + + + + def product(self,x,y): + r""" + SETUP:: + + sage: from mjo.eja.eja_algebra import (JordanSpinEJA, + ....: ComplexHermitianEJA, + ....: DirectSumEJA) + + TESTS:: + + sage: set_random_seed() + sage: J1 = JordanSpinEJA(3, field=QQ) + sage: J2 = ComplexHermitianEJA(2, field=QQ, orthonormalize=False) + sage: J = DirectSumEJA(J1,J2) + sage: J.random_element()*J.random_element() in J + True + + """ + xv = self._cartprod_algebra.from_vector(x.to_vector()) + yv = self._cartprod_algebra.from_vector(y.to_vector()) + return self.from_vector((xv*yv).to_vector()) + + + def cartesian_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.cartesian_factors() + (Euclidean Jordan algebra of dimension 2 over Rational Field, + Euclidean Jordan algebra of dimension 3 over Rational Field) + + """ + return self._cartprod_algebra.cartesian_factors() -# """ -# return self._factors # def projections(self): # r"""