X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=mjo%2Feja%2Feja_element.py;h=e30dbb13f39d06b6d49c6bfd34c829ab30d6112d;hb=e4d568e25c62d79a2dbe34b81ee4dc21edf09316;hp=a5880c4c4c3fc083726fbf2adcef6ee071e7cd63;hpb=fe2af66b109e9487a59f21d5b67bb5c4aafdc98d;p=sage.d.git diff --git a/mjo/eja/eja_element.py b/mjo/eja/eja_element.py index a5880c4..e30dbb1 100644 --- a/mjo/eja/eja_element.py +++ b/mjo/eja/eja_element.py @@ -2,15 +2,10 @@ from sage.matrix.constructor import matrix from sage.modules.free_module import VectorSpace from sage.modules.with_basis.indexed_element import IndexedFreeModuleElement -# TODO: make this unnecessary somehow. -from sage.misc.lazy_import import lazy_import -lazy_import('mjo.eja.eja_algebra', 'FiniteDimensionalEuclideanJordanAlgebra') -lazy_import('mjo.eja.eja_element_subalgebra', - 'FiniteDimensionalEuclideanJordanElementSubalgebra') -from mjo.eja.eja_operator import FiniteDimensionalEuclideanJordanAlgebraOperator +from mjo.eja.eja_operator import FiniteDimensionalEJAOperator from mjo.eja.eja_utils import _mat2vec -class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): +class FiniteDimensionalEJAElement(IndexedFreeModuleElement): """ An element of a Euclidean Jordan algebra. """ @@ -232,9 +227,9 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): Ditto for the quaternions:: - sage: J = QuaternionHermitianEJA(3) + sage: J = QuaternionHermitianEJA(2) sage: J.one().inner_product(J.one()) - 3 + 2 TESTS: @@ -345,6 +340,8 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): sage: from mjo.eja.eja_algebra import (JordanSpinEJA, ....: TrivialEJA, + ....: RealSymmetricEJA, + ....: ComplexHermitianEJA, ....: random_eja) EXAMPLES:: @@ -392,6 +389,37 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): sage: x,y = J.random_elements(2) sage: (x*y).det() == x.det()*y.det() True + + The determinant in matrix algebras is just the usual determinant:: + + sage: set_random_seed() + sage: X = matrix.random(QQ,3) + sage: X = X + X.T + sage: J1 = RealSymmetricEJA(3) + sage: J2 = RealSymmetricEJA(3,field=QQ,orthonormalize=False) + sage: expected = X.det() + sage: actual1 = J1(X).det() + sage: actual2 = J2(X).det() + sage: actual1 == expected + True + sage: actual2 == expected + True + + :: + + sage: set_random_seed() + sage: J1 = ComplexHermitianEJA(2) + sage: J2 = ComplexHermitianEJA(2,field=QQ,orthonormalize=False) + sage: X = matrix.random(GaussianIntegers(), 2) + sage: X = X + X.H + sage: expected = AA(X.det()) + sage: actual1 = J1(J1.real_embed(X)).det() + sage: actual2 = J2(J2.real_embed(X)).det() + sage: expected == actual1 + True + sage: expected == actual2 + True + """ P = self.parent() r = P.rank() @@ -416,8 +444,13 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): ALGORITHM: - We appeal to the quadratic representation as in Koecher's - Theorem 12 in Chapter III, Section 5. + In general we appeal to the quadratic representation as in + Koecher's Theorem 12 in Chapter III, Section 5. But if the + parent algebra's "characteristic polynomial of" coefficients + happen to be cached, then we use Proposition II.2.4 in Faraut + and Korányi which gives a formula for the inverse based on the + characteristic polynomial and the Cayley-Hamilton theorem for + Euclidean Jordan algebras:: SETUP:: @@ -487,26 +520,31 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): ....: x.operator().inverse()(J.one()) == x.inverse() ) True - Proposition II.2.4 in Faraut and Korányi gives a formula for - the inverse based on the characteristic polynomial and the - Cayley-Hamilton theorem for Euclidean Jordan algebras:: + Check that the fast (cached) and slow algorithms give the same + answer:: - sage: set_random_seed() - sage: J = ComplexHermitianEJA(3) - sage: x = J.random_element() - sage: while not x.is_invertible(): - ....: x = J.random_element() - sage: r = J.rank() - sage: a = x.characteristic_polynomial().coefficients(sparse=False) - sage: expected = (-1)^(r+1)/x.det() - sage: expected *= sum( a[i+1]*x^i for i in range(r) ) - sage: x.inverse() == expected + sage: set_random_seed() # long time + sage: J = random_eja(field=QQ, orthonormalize=False) # long time + sage: x = J.random_element() # long time + sage: while not x.is_invertible(): # long time + ....: x = J.random_element() # long time + sage: slow = x.inverse() # long time + sage: _ = J._charpoly_coefficients() # long time + sage: fast = x.inverse() # long time + sage: slow == fast # long time True - """ if not self.is_invertible(): raise ValueError("element is not invertible") + if self.parent()._charpoly_coefficients.is_in_cache(): + # We can invert using our charpoly if it will be fast to + # compute. If the coefficients are cached, our rank had + # better be too! + r = self.parent().rank() + a = self.characteristic_polynomial().coefficients(sparse=False) + return (-1)**(r+1)*sum(a[i+1]*self**i for i in range(r))/self.det() + return (~self.quadratic_representation())(self) @@ -520,9 +558,14 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): zero, but we need the characteristic polynomial for the determinant. The minimal polynomial is a lot easier to get, so we use Corollary 2 in Chapter V of Koecher to check - whether or not the paren't algebra's zero element is a root + whether or not the parent algebra's zero element is a root of this element's minimal polynomial. + That is... unless the coefficients of our algebra's + "characteristic polynomial of" function are already cached! + In that case, we just use the determinant (which will be fast + as a result). + Beware that we can't use the superclass method, because it relies on the algebra being associative. @@ -546,6 +589,18 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): sage: (not J.is_trivial()) and J.zero().is_invertible() False + Test that the fast (cached) and slow algorithms give the same + answer:: + + sage: set_random_seed() # long time + sage: J = random_eja(field=QQ, orthonormalize=False) # long time + sage: x = J.random_element() # long time + sage: slow = x.is_invertible() # long time + sage: _ = J._charpoly_coefficients() # long time + sage: fast = x.is_invertible() # long time + sage: slow == fast # long time + True + """ if self.is_zero(): if self.parent().is_trivial(): @@ -553,6 +608,11 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): else: return False + if self.parent()._charpoly_coefficients.is_in_cache(): + # The determinant will be quicker than computing the minimal + # polynomial from scratch, most likely. + return (not self.det().is_zero()) + # In fact, we only need to know if the constant term is non-zero, # so we can pass in the field's zero element instead. zero = self.base_ring().zero() @@ -770,10 +830,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): ALGORITHM: - For now, we skip the messy minimal polynomial computation - and instead return the dimension of the vector space spanned - by the powers of this element. The latter is a bit more - straightforward to compute. + ......... SETUP:: @@ -821,12 +878,59 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): True """ - if self.is_zero() and not self.parent().is_trivial(): + n = self.parent().dimension() + + if n == 0: + # The minimal polynomial is an empty product, i.e. the + # constant polynomial "1" having degree zero. + return 0 + elif self.is_zero(): # The minimal polynomial of zero in a nontrivial algebra - # is "t"; in a trivial algebra it's "1" by convention - # (it's an empty product). + # is "t", and is of degree one. return 1 - return self.subalgebra_generated_by().dimension() + elif n == 1: + # If this is a nonzero element of a nontrivial algebra, it + # has degree at least one. It follows that, in an algebra + # of dimension one, the degree must be actually one. + return 1 + + # BEWARE: The subalgebra_generated_by() method uses the result + # of this method to construct a basis for the subalgebra. That + # means, in particular, that we cannot implement this method + # as ``self.subalgebra_generated_by().dimension()``. + + # Algorithm: keep appending (vector representations of) powers + # self as rows to a matrix and echelonizing it. When its rank + # stops increasing, we've reached a redundancy. + + # Given the special cases above, we can assume that "self" is + # nonzero, the algebra is nontrivial, and that its dimension + # is at least two. + M = matrix([(self.parent().one()).to_vector()]) + old_rank = 1 + + # Specifying the row-reduction algorithm can e.g. help over + # AA because it avoids the RecursionError that gets thrown + # when we have to look too hard for a root. + # + # Beware: QQ supports an entirely different set of "algorithm" + # keywords than do AA and RR. + algo = None + from sage.rings.all import QQ + if self.parent().base_ring() is not QQ: + algo = "scaled_partial_pivoting" + + for d in range(1,n): + M = matrix(M.rows() + [(self**d).to_vector()]) + M.echelonize(algo) + new_rank = M.rank() + if new_rank == old_rank: + return new_rank + else: + old_rank = new_rank + + return n + def left_matrix(self): @@ -932,7 +1036,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): sage: n_max = RealSymmetricEJA._max_random_instance_size() sage: n = ZZ.random_element(1, n_max) sage: J1 = RealSymmetricEJA(n) - sage: J2 = RealSymmetricEJA(n,normalize_basis=False) + sage: J2 = RealSymmetricEJA(n,orthonormalize=False) sage: X = random_matrix(AA,n) sage: X = X*X.transpose() sage: x1 = J1(X) @@ -953,20 +1057,21 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): # in the "normal" case without us having to think about it. return self.operator().minimal_polynomial() - A = self.subalgebra_generated_by(orthonormalize_basis=False) + A = self.subalgebra_generated_by(orthonormalize=False) return A(self).operator().minimal_polynomial() - def natural_representation(self): + def to_matrix(self): """ - Return a more-natural representation of this element. + Return an (often more natural) representation of this element as a + matrix. - 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" representation of this element as a Hermitian - matrix, if it has one. If not, you get the usual representation. + Every finite-dimensional Euclidean Jordan Algebra is a direct + sum of five simple algebras, four of which comprise Hermitian + matrices. This method returns a "natural" matrix + representation of this element as either a Hermitian matrix or + column vector. SETUP:: @@ -978,7 +1083,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): sage: J = ComplexHermitianEJA(3) sage: J.one() e0 + e3 + e8 - sage: J.one().natural_representation() + sage: J.one().to_matrix() [1 0 0 0 0 0] [0 1 0 0 0 0] [0 0 1 0 0 0] @@ -988,31 +1093,27 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): :: - sage: J = QuaternionHermitianEJA(3) + sage: J = QuaternionHermitianEJA(2) sage: J.one() - e0 + e5 + e14 - sage: J.one().natural_representation() - [1 0 0 0 0 0 0 0 0 0 0 0] - [0 1 0 0 0 0 0 0 0 0 0 0] - [0 0 1 0 0 0 0 0 0 0 0 0] - [0 0 0 1 0 0 0 0 0 0 0 0] - [0 0 0 0 1 0 0 0 0 0 0 0] - [0 0 0 0 0 1 0 0 0 0 0 0] - [0 0 0 0 0 0 1 0 0 0 0 0] - [0 0 0 0 0 0 0 1 0 0 0 0] - [0 0 0 0 0 0 0 0 1 0 0 0] - [0 0 0 0 0 0 0 0 0 1 0 0] - [0 0 0 0 0 0 0 0 0 0 1 0] - [0 0 0 0 0 0 0 0 0 0 0 1] + e0 + e5 + sage: J.one().to_matrix() + [1 0 0 0 0 0 0 0] + [0 1 0 0 0 0 0 0] + [0 0 1 0 0 0 0 0] + [0 0 0 1 0 0 0 0] + [0 0 0 0 1 0 0 0] + [0 0 0 0 0 1 0 0] + [0 0 0 0 0 0 1 0] + [0 0 0 0 0 0 0 1] """ - B = self.parent().natural_basis() - W = self.parent().natural_basis_space() + B = self.parent().matrix_basis() + W = self.parent().matrix_space() # This is just a manual "from_vector()", but of course # matrix spaces aren't vector spaces in sage, so they # don't have a from_vector() method. - return W.linear_combination(zip(B,self.to_vector())) + return W.linear_combination( zip(B, self.to_vector()) ) def norm(self): @@ -1065,10 +1166,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): P = self.parent() left_mult_by_self = lambda y: self*y L = P.module_morphism(function=left_mult_by_self, codomain=P) - return FiniteDimensionalEuclideanJordanAlgebraOperator( - P, - P, - L.matrix() ) + return FiniteDimensionalEJAOperator(P, P, L.matrix() ) def quadratic_representation(self, other=None): @@ -1260,13 +1358,13 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): [(0, f2), (1, f0)] """ - A = self.subalgebra_generated_by(orthonormalize_basis=True) + A = self.subalgebra_generated_by(orthonormalize=True) result = [] for (evalue, proj) in A(self).operator().spectral_decomposition(): result.append( (evalue, proj(A.one()).superalgebra_element()) ) return result - def subalgebra_generated_by(self, orthonormalize_basis=False): + def subalgebra_generated_by(self, **kwargs): """ Return the associative subalgebra of the parent EJA generated by this element. @@ -1313,7 +1411,14 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement): True """ - return FiniteDimensionalEuclideanJordanElementSubalgebra(self, orthonormalize_basis) + from mjo.eja.eja_subalgebra import FiniteDimensionalEJASubalgebra + powers = tuple( self**k for k in range(self.degree()) ) + A = FiniteDimensionalEJASubalgebra(self.parent(), + powers, + associative=True, + **kwargs) + A.one.set_cache(A(self.parent().one())) + return A def subalgebra_idempotent(self):