X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=mjo%2Feja%2Feja_algebra.py;h=055bbbda83d7576fda8f5e1794d759c2eee467bc;hb=f807c8c9046723a059eb4eea03dff89293510160;hp=01deb1b65140675ffbc81c9c673c21e6d4a2e77f;hpb=8e4c84f17025f6c3f2d049902a6a01f954d484ee;p=sage.d.git diff --git a/mjo/eja/eja_algebra.py b/mjo/eja/eja_algebra.py index 01deb1b..055bbbd 100644 --- a/mjo/eja/eja_algebra.py +++ b/mjo/eja/eja_algebra.py @@ -235,163 +235,6 @@ 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(): - return R.one() - if i > self.rank(): - # Guaranteed by theory - return R.zero() - - # Danger: the in-place modification is done for performance - # reasons (reconstructing a matrix with huge polynomial - # entries is slow), but I don't know how cached_method works, - # so it's highly possible that we're modifying some global - # list variable by reference, here. In other words, you - # probably shouldn't call this method twice on the same - # algebra, at the same time, in two threads - Ai_orig = A_of_x.column(i) - A_of_x.set_column(i,xr) - numerator = A_of_x.det() - A_of_x.set_column(i,Ai_orig) - - # We're relying on the theory here to ensure that each a_i is - # indeed back in R, and the added negative signs are to make - # the whole charpoly expression sum to zero. - return R(-numerator/detA) - - - @cached_method - def _charpoly_matrix_system(self): - """ - 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. - """ - 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): """ @@ -436,20 +279,21 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): r = self.rank() n = self.dimension() - # The list of coefficient polynomials a_0, a_1, a_2, ..., a_n. - a = [ self._charpoly_coeff(i) for i in range(r+1) ] + # 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) + 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 range(len(a)) ) + return (t**r + sum( a[k]*(t**k) for k in range(r) )) def inner_product(self, x, y): @@ -826,18 +670,60 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): return cls(n, field, **kwargs) @cached_method - def rank(self): + def _charpoly_coefficients(self): + r""" + The `r` polynomial coefficients of the "characteristic polynomial + of" function. """ - Return the rank of this EJA. + n = self.dimension() + var_names = [ "X" + str(z) for z in range(1,n+1) ] + R = PolynomialRing(self.base_ring(), var_names) + 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. - We first compute the polynomial "column matrices" `p_{k}` that - evaluate to `x^k` on the coordinates of `x`. Then, we begin - adding them to a matrix one at a time, and trying to solve the - system that makes `p_{0}`,`p_{1}`,..., `p_{s-1}` add up to - `p_{s}`. This will succeed only when `s` is the rank of the - algebra, as proven in a recent draft paper of mine. + 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:: @@ -918,45 +804,8 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): sage: J.rank.clear_cache() sage: J.rank() 2 - """ - n = self.dimension() - if n == 0: - return 0 - elif n == 1: - return 1 - - var_names = [ "X" + str(z) for z in range(1,n+1) ] - R = PolynomialRing(self.base_ring(), var_names) - vars = R.gens() - - 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(R, n, n, L_x_i_j) - x_powers = [ vars[k]*(L_x**k)*self.one().to_vector() - for k in range(n) ] - - # Can assume n >= 2 - M = matrix([x_powers[0]]) - old_rank = 1 - - for d in range(1,n): - M = matrix(M.rows() + [x_powers[d]]) - M.echelonize() - # TODO: we've basically solved the system here. - # We should save the echelonized matrix somehow - # so that it can be reused in the charpoly method. - new_rank = M.rank() - if new_rank == old_rank: - return new_rank - else: - old_rank = new_rank - - return n + return len(self._charpoly_coefficients()) def vector_space(self): @@ -1054,7 +903,7 @@ class HadamardEJA(FiniteDimensionalEuclideanJordanAlgebra): return x.to_vector().inner_product(y.to_vector()) -def random_eja(field=AA, nontrivial=False): +def random_eja(field=AA): """ Return a "random" finite-dimensional Euclidean Jordan Algebra. @@ -1068,21 +917,17 @@ def random_eja(field=AA, nontrivial=False): Euclidean Jordan algebra of dimension... """ - eja_classes = [HadamardEJA, - JordanSpinEJA, - RealSymmetricEJA, - ComplexHermitianEJA, - QuaternionHermitianEJA] - if not nontrivial: - eja_classes.append(TrivialEJA) - classname = choice(eja_classes) + classname = choice([TrivialEJA, + HadamardEJA, + JordanSpinEJA, + RealSymmetricEJA, + ComplexHermitianEJA, + QuaternionHermitianEJA]) return classname.random_instance(field=field) - - class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra): @staticmethod def _max_test_case_size(): @@ -1096,7 +941,7 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra): 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 @@ -1125,7 +970,7 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra): @cached_method - def rank(self): + def _charpoly_coefficients(self): r""" Override the parent method with something that tries to compute over a faster (non-extension) field. @@ -1133,7 +978,7 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra): 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).rank() + return super(MatrixEuclideanJordanAlgebra, self)._charpoly_coefficients() else: basis = ( (b/n) for (b,n) in zip(self.natural_basis(), self._basis_normalizers) ) @@ -1144,39 +989,26 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra): J = MatrixEuclideanJordanAlgebra(QQ, basis, normalize_basis=False) - return J.rank() - - @cached_method - def _charpoly_coeff(self, i): - """ - 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 zip(self.natural_basis(), - self._basis_normalizers) ) - - # Do this over the rationals and convert back at the end. - J = MatrixEuclideanJordanAlgebra(QQ, - basis, - normalize_basis=False) - (_,x,_,_) = J._charpoly_matrix_system() - p = J._charpoly_coeff(i) - # p might be missing some vars, have to substitute "optionally" - pairs = zip(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()) + 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 @@ -1194,6 +1026,9 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra): # 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. + if len(basis) == 0: + return [] + field = basis[0].base_ring() dimension = basis[0].nrows() @@ -1351,6 +1186,11 @@ class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra): sage: x.operator().matrix().is_symmetric() True + We can construct the (trivial) algebra of rank zero:: + + sage: RealSymmetricEJA(0) + Euclidean Jordan algebra of dimension 0 over Algebraic Real Field + """ @classmethod def _denormalized_basis(cls, n, field): @@ -1624,6 +1464,11 @@ class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra): sage: x.operator().matrix().is_symmetric() True + We can construct the (trivial) algebra of rank zero:: + + sage: ComplexHermitianEJA(0) + Euclidean Jordan algebra of dimension 0 over Algebraic Real Field + """ @classmethod @@ -1919,6 +1764,11 @@ class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra): sage: x.operator().matrix().is_symmetric() True + We can construct the (trivial) algebra of rank zero:: + + sage: QuaternionHermitianEJA(0) + Euclidean Jordan algebra of dimension 0 over Algebraic Real Field + """ @classmethod def _denormalized_basis(cls, n, field):