From: Michael Orlitzky Date: Sat, 28 Nov 2020 04:14:31 +0000 (-0500) Subject: eja: more work on realizing the new constructor. X-Git-Url: http://gitweb.michael.orlitzky.com/?p=sage.d.git;a=commitdiff_plain;h=8207ab8350d947b557bf6682da9e3c3ffe638523 eja: more work on realizing the new constructor. --- diff --git a/mjo/eja/eja_algebra.py b/mjo/eja/eja_algebra.py index 9eba669..14666c2 100644 --- a/mjo/eja/eja_algebra.py +++ b/mjo/eja/eja_algebra.py @@ -34,7 +34,9 @@ from mjo.eja.eja_operator import FiniteDimensionalEuclideanJordanAlgebraOperator from mjo.eja.eja_utils import _mat2vec class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): - + r""" + The lowest-level class for representing a Euclidean Jordan algebra. + """ def _coerce_map_from_base_ring(self): """ Disable the map from the base ring into the algebra. @@ -61,13 +63,23 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): def __init__(self, field, - mult_table, + multiplication_table, + inner_product_table, prefix='e', category=None, matrix_basis=None, check_field=True, check_axioms=True): """ + INPUT: + + * field -- the scalar field for this algebra (must be real) + + * multiplication_table -- the multiplication table for this + algebra's implicit basis. Only the lower-triangular portion + of the table is used, since the multiplication is assumed + to be commutative. + SETUP:: sage: from mjo.eja.eja_algebra import ( @@ -85,6 +97,24 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): sage: x*y == y*x True + An error is raised if the Jordan product is not commutative:: + + sage: JP = ((1,2),(0,0)) + sage: IP = ((1,0),(0,1)) + sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,JP,IP) + Traceback (most recent call last): + ... + ValueError: Jordan product is not commutative + + An error is raised if the inner-product is not commutative:: + + sage: JP = ((1,0),(0,1)) + sage: IP = ((1,2),(0,0)) + sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,JP,IP) + Traceback (most recent call last): + ... + ValueError: inner-product is not commutative + TESTS: The ``field`` we're given must be real with ``check_field=True``:: @@ -96,11 +126,26 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): The multiplication table must be square with ``check_axioms=True``:: - sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,((),())) + sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,((),()),((1,),)) Traceback (most recent call last): ... ValueError: multiplication table is not square + The multiplication and inner-product tables must be the same + size (and in particular, the inner-product table must also be + square) with ``check_axioms=True``:: + + sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,((1,),),(())) + Traceback (most recent call last): + ... + ValueError: multiplication and inner-product tables are + different sizes + sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,((1,),),((1,2),)) + Traceback (most recent call last): + ... + ValueError: multiplication and inner-product tables are + different sizes + """ if check_field: if not field.is_subring(RR): @@ -109,12 +154,35 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): # we've specified a real embedding. raise ValueError("scalar field is not real") - # The multiplication table had better be square - n = len(mult_table) + + # The multiplication and inner-product tables should be square + # if the user wants us to verify them. And we verify them as + # soon as possible, because we want to exploit their symmetry. + n = len(multiplication_table) if check_axioms: - if not all( len(l) == n for l in mult_table ): + if not all( len(l) == n for l in multiplication_table ): raise ValueError("multiplication table is not square") + # If the multiplication table is square, we can check if + # the inner-product table is square by comparing it to the + # multiplication table's dimensions. + msg = "multiplication and inner-product tables are different sizes" + if not len(inner_product_table) == n: + raise ValueError(msg) + + if not all( len(l) == n for l in inner_product_table ): + raise ValueError(msg) + + if not all( multiplication_table[j][i] + == multiplication_table[i][j] + for i in range(n) + for j in range(i+1) ): + raise ValueError("Jordan product is not commutative") + if not all( inner_product_table[j][i] + == inner_product_table[i][j] + for i in range(n) + for j in range(i+1) ): + raise ValueError("inner-product is not commutative") self._matrix_basis = matrix_basis if category is None: @@ -134,19 +202,28 @@ 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. + # + # Note: we take advantage of symmetry here, and only store + # the lower-triangular portion of the table. self._multiplication_table = [ [ self.vector_space().zero() - for i in range(n) ] - for j in range(n) ] - # take advantage of symmetry + for j in range(i+1) ] + for i in range(n) ] + for i in range(n): for j in range(i+1): - elt = self.from_vector(mult_table[i][j]) + elt = self.from_vector(multiplication_table[i][j]) self._multiplication_table[i][j] = elt - self._multiplication_table[j][i] = elt + + # Save our inner product as a matrix, since the efficiency of + # matrix multiplication will usually outweigh the fact that we + # have to store a redundant upper- or lower-triangular part. + # Pre-cache the fact that these are Hermitian (real symmetric, + # in fact) in case some e.g. matrix multiplication routine can + # take advantage of it. + self._inner_product_matrix = matrix(field, inner_product_table) + self._inner_product_matrix._cache = {'hermitian': False} 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(): @@ -253,7 +330,12 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): return fmt.format(self.dimension(), self.base_ring()) def product_on_basis(self, i, j): - return self._multiplication_table[i][j] + # We only stored the lower-triangular portion of the + # multiplication table. + if j <= i: + return self._multiplication_table[i][j] + else: + return self._multiplication_table[j][i] def _is_commutative(self): r""" @@ -509,10 +591,20 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): +----++----+----+----+----+ """ - M = list(self._multiplication_table) # copy - for i in range(len(M)): - # M had better be "square" + n = self.dimension() + M = [ [ self.zero() for j in range(n) ] + for i in range(n) ] + for i in range(n): + for j in range(i+1): + M[i][j] = self._multiplication_table[i][j] + M[j][i] = M[i][j] + + for i in range(n): + # Prepend the left "header" column entry Can't do this in + # the loop because it messes up the symmetry. M[i] = [self.monomial(i)] + M[i] + + # Prepend the header row. M = [["*"] + list(self.gens())] + M return table(M, header_row=True, header_column=True, frame=True) @@ -1026,6 +1118,27 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): Element = FiniteDimensionalEuclideanJordanAlgebraElement class RationalBasisEuclideanJordanAlgebraNg(FiniteDimensionalEuclideanJordanAlgebra): + r""" + New class for algebras whose supplied basis elements have all rational entries. + + SETUP:: + + sage: from mjo.eja.eja_algebra import BilinearFormEJA + + EXAMPLES: + + The supplied basis is orthonormalized by default:: + + sage: B = matrix(QQ, [[1, 0, 0], [0, 25, -32], [0, -32, 41]]) + sage: J = BilinearFormEJA(B) + sage: J.matrix_basis() + ( + [1] [ 0] [ 0] + [0] [1/5] [32/5] + [0], [ 0], [ 5] + ) + + """ def __init__(self, field, basis, @@ -1054,26 +1167,26 @@ class RationalBasisEuclideanJordanAlgebraNg(FiniteDimensionalEuclideanJordanAlge V = VectorSpace(field, degree) - # Compute this from "Q" (obtained from Gram-Schmidt) below as - # R = Q.solve_right(A), where the rows of "Q" are the - # orthonormalized vector_basis and and the rows of "A" are the - # original vector_basis. - self._deorthonormalization_matrix = None + # If we were asked to orthonormalize, and if the orthonormal + # basis is different from the given one, then we also want to + # compute multiplication and inner-product tables for the + # deorthonormalized basis. These can be used later to + # construct a deorthonormalized copy of this algebra over QQ + # in which several operations are much faster. + self._deortho_multiplication_table = None + self._deortho_inner_product_table = None if orthonormalize: from mjo.eja.eja_utils import gram_schmidt - A = matrix(field, vector_basis) vector_basis = gram_schmidt(vector_basis, inner_product) W = V.span_of_basis( vector_basis ) - Q = matrix(field, vector_basis) - # A = QR <==> A.T == R.T*Q.T - # So, Q.solve_right() is equivalent to the Q.T.solve_left() - # that we want. - self._deorthonormalization_matrix = Q.solve_right(A) + + # Normalize the "matrix" basis, too! + basis = vector_basis if basis_is_matrices: from mjo.eja.eja_utils import _vec2mat - basis = tuple( map(_vec2mat,vector_basis) ) + basis = tuple( map(_vec2mat,basis) ) W = V.span_of_basis( vector_basis ) @@ -1107,8 +1220,6 @@ class RationalBasisEuclideanJordanAlgebraNg(FiniteDimensionalEuclideanJordanAlge 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() @@ -1117,6 +1228,7 @@ class RationalBasisEuclideanJordanAlgebraNg(FiniteDimensionalEuclideanJordanAlge super().__init__(field, mult_table, + ip_table, prefix, category, basis, # matrix basis @@ -1287,7 +1399,7 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra): 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)] + ip_table = [[field.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): @@ -1303,14 +1415,9 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra): except: pass - try: - # HACK PART DEUX - self._inner_product_matrix = matrix(field,ip_table) - except: - pass - super(MatrixEuclideanJordanAlgebra, self).__init__(field, mult_table, + ip_table, matrix_basis=basis, **kwargs) @@ -2272,7 +2379,9 @@ class BilinearFormEJA(RationalBasisEuclideanJordanAlgebraNg, 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:: + choice of basis, as does anything utilizing the bilinear form + matrix. We opt not to orthonormalize the basis, because if we + did, we would have to normalize the `s_{i}` in a similar manner:: sage: set_random_seed() sage: n = ZZ.random_element(5) @@ -2281,10 +2390,10 @@ class BilinearFormEJA(RationalBasisEuclideanJordanAlgebraNg, sage: B22 = M.transpose()*M sage: B = block_matrix(2,2,[ [B11,0 ], ....: [0, B22 ] ]) - sage: J = BilinearFormEJA(B) + sage: J = BilinearFormEJA(B, orthonormalize=False) 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())) + sage: sis = [ J( V([0] + (M.inverse()*ei).list()).column() ) ....: for ei in eis ] sage: actual = [ sis[i]*sis[j] ....: for i in range(n-1) @@ -2468,9 +2577,10 @@ class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra, """ def __init__(self, field=AA, **kwargs): mult_table = [] - self._inner_product_matrix = matrix(field,0) + ip_table = [] super(TrivialEJA, self).__init__(field, mult_table, + ip_table, check_axioms=False, **kwargs) # The rank is zero using my definition, namely the dimension of the @@ -2545,8 +2655,12 @@ class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra): 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 = None super(DirectSumEJA, self).__init__(field, mult_table, + ip_table, check_axioms=False, **kwargs) self.rank.set_cache(J1.rank() + J2.rank()) diff --git a/mjo/eja/eja_subalgebra.py b/mjo/eja/eja_subalgebra.py index e7e559f..6eca475 100644 --- a/mjo/eja/eja_subalgebra.py +++ b/mjo/eja/eja_subalgebra.py @@ -191,7 +191,6 @@ class FiniteDimensionalEuclideanJordanSubalgebra(FiniteDimensionalEuclideanJorda product_vector = V.from_vector(product.to_vector()) mult_table[i][j] = W.coordinate_vector(product_vector) - self._inner_product_matrix = matrix(field, ip_table) matrix_basis = tuple( b.to_matrix() for b in basis ) @@ -200,6 +199,7 @@ class FiniteDimensionalEuclideanJordanSubalgebra(FiniteDimensionalEuclideanJorda fdeja = super(FiniteDimensionalEuclideanJordanSubalgebra, self) fdeja.__init__(field, mult_table, + ip_table, prefix=prefix, category=category, matrix_basis=matrix_basis,