From de451def1161cc9dfefcfc125523029881cb160a Mon Sep 17 00:00:00 2001 From: Michael Orlitzky Date: Sun, 6 Dec 2020 11:03:58 -0500 Subject: [PATCH] eja: begin dropping support for vector representations. --- mjo/eja/eja_algebra.py | 111 ++++++++++-------------------- mjo/eja/eja_element_subalgebra.py | 49 +++++++------ 2 files changed, 61 insertions(+), 99 deletions(-) diff --git a/mjo/eja/eja_algebra.py b/mjo/eja/eja_algebra.py index 1250fbd..ea62da8 100644 --- a/mjo/eja/eja_algebra.py +++ b/mjo/eja/eja_algebra.py @@ -113,60 +113,36 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): # we see in things like x = 1*e1 + 2*e2. vector_basis = basis - from sage.structure.element import is_Matrix - basis_is_matrices = False - degree = 0 if n > 0: - if is_Matrix(basis[0]): - if basis[0].is_square(): - # TODO: this ugly is_square() hack works around the problem - # of passing to_matrix()ed vectors in as the basis from a - # subalgebra. They aren't REALLY matrices, at least not of - # the type that we assume here... Ugh. - basis_is_matrices = True - from mjo.eja.eja_utils import _vec2mat - vector_basis = tuple( map(_mat2vec,basis) ) - degree = basis[0].nrows()**2 - else: - # convert from column matrices to vectors, yuck - basis = tuple( map(_mat2vec,basis) ) - vector_basis = basis - degree = basis[0].degree() - else: - degree = basis[0].degree() + # Works on both column and square matrices... + degree = len(basis[0].list()) - # Build an ambient space that fits... + # Build an ambient space that fits our matrix basis when + # written out as "long vectors." V = VectorSpace(field, degree) - # We overwrite the name "vector_basis" in a second, but never modify it - # in place, to this effectively makes a copy of it. - deortho_vector_basis = vector_basis + # The matrix that will hole the orthonormal -> unorthonormal + # coordinate transformation. self._deortho_matrix = None if orthonormalize: - from mjo.eja.eja_utils import gram_schmidt - if basis_is_matrices: - vector_ip = lambda x,y: inner_product(_vec2mat(x), _vec2mat(y)) - vector_basis = gram_schmidt(vector_basis, vector_ip) - else: - vector_basis = gram_schmidt(vector_basis, inner_product) + # Save a copy of the un-orthonormalized basis for later. + # Convert it to ambient V (vector) coordinates while we're + # at it, because we'd have to do it later anyway. + deortho_vector_basis = tuple( V(b.list()) for b in basis ) - # Normalize the "matrix" basis, too! - basis = vector_basis - - if basis_is_matrices: - basis = tuple( map(_vec2mat,basis) ) + from mjo.eja.eja_utils import gram_schmidt + basis = gram_schmidt(basis, inner_product) - # Save the matrix "basis" for later... this is the last time we'll - # reference it in this constructor. - if basis_is_matrices: - self._matrix_basis = basis - else: - MS = MatrixSpace(self.base_ring(), degree, 1) - self._matrix_basis = tuple( MS(b) for b in basis ) + # Save the (possibly orthonormalized) matrix basis for + # later... + self._matrix_basis = basis - # Now create the vector space for the algebra... + # Now create the vector space for the algebra, which will have + # its own set of non-ambient coordinates (in terms of the + # supplied basis). + vector_basis = tuple( V(b.list()) for b in basis ) W = V.span_of_basis( vector_basis, check=check_axioms) if orthonormalize: @@ -183,36 +159,24 @@ 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._multiplication_table = [ [0 for j in range(i+1)] for i in range(n) ] + self._multiplication_table = [ [0 for j in range(i+1)] + for i in range(n) ] - print("vector_basis:") - print(vector_basis) # Note: the Jordan and inner-products are defined in terms # of the ambient basis. It's important that their arguments # are in ambient coordinates as well. for i in range(n): for j in range(i+1): # ortho basis w.r.t. ambient coords - q_i = vector_basis[i] - q_j = vector_basis[j] - - if basis_is_matrices: - q_i = _vec2mat(q_i) - q_j = _vec2mat(q_j) + q_i = basis[i] + q_j = basis[j] elt = jordan_product(q_i, q_j) ip = inner_product(q_i, q_j) - if basis_is_matrices: - # do another mat2vec because the multiplication - # table is in terms of vectors - elt = _mat2vec(elt) - - # TODO: the jordan product turns things back into - # matrices here even if they're supposed to be - # vectors. ugh. Can we get rid of vectors all together - # please? - elt = W.coordinate_vector(elt) + # The jordan product returns a matrixy answer, so we + # have to convert it to the algebra coordinates. + 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 @@ -2222,11 +2186,16 @@ class HadamardEJA(ConcreteEJA): """ def __init__(self, n, **kwargs): - def jordan_product(x,y): - P = x.parent() - return P(tuple( xi*yi for (xi,yi) in zip(x,y) )) - def inner_product(x,y): - return x.inner_product(y) + if n == 0: + jordan_product = lambda x,y: x + inner_product = lambda x,y: x + else: + def jordan_product(x,y): + P = x.parent() + return P( xi*yi for (xi,yi) in zip(x,y) ) + + def inner_product(x,y): + return (x.T*y)[0,0] # New defaults for keyword arguments. Don't orthonormalize # because our basis is already orthonormal with respect to our @@ -2237,12 +2206,8 @@ class HadamardEJA(ConcreteEJA): if "orthonormalize" not in kwargs: kwargs["orthonormalize"] = False if "check_axioms" not in kwargs: kwargs["check_axioms"] = False - - standard_basis = FreeModule(ZZ, n).basis() - super(HadamardEJA, self).__init__(standard_basis, - jordan_product, - inner_product, - **kwargs) + column_basis = tuple( b.column() for b in FreeModule(ZZ, n).basis() ) + super().__init__(column_basis, jordan_product, inner_product, **kwargs) self.rank.set_cache(n) if n == 0: diff --git a/mjo/eja/eja_element_subalgebra.py b/mjo/eja/eja_element_subalgebra.py index dceb3b4..846c13b 100644 --- a/mjo/eja/eja_element_subalgebra.py +++ b/mjo/eja/eja_element_subalgebra.py @@ -9,36 +9,33 @@ class FiniteDimensionalEJAElementSubalgebra(FiniteDimensionalEJASubalgebra): def __init__(self, elt, orthonormalize=True, **kwargs): superalgebra = elt.parent() + # TODO: going up to the superalgebra dimension here is + # overkill. We should append p vectors as rows to a matrix + # and continually rref() it until the rank stops going + # up. When n=10 but the dimension of the algebra is 1, that + # can save a shitload of time (especially over AA). powers = tuple( elt**k for k in range(superalgebra.dimension()) ) power_vectors = ( p.to_vector() for p in powers ) P = matrix(superalgebra.base_ring(), power_vectors) - if orthonormalize: - basis = powers # let god sort 'em out - else: - # Echelonize the matrix ourselves, because otherwise the - # call to P.pivot_rows() below can choose a non-optimal - # row-reduction algorithm. In particular, scaling can - # 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 - if superalgebra.base_ring() is not QQ: - algo = "scaled_partial_pivoting" - P.echelonize(algorithm=algo) - - # In this case, we just need to figure out which elements - # of the "powers" list are redundant... First compute the - # vector subspace spanned by the powers of the given - # element. - - # Figure out which powers form a linearly-independent set. - ind_rows = P.pivot_rows() - - # Pick those out of the list of all powers. - basis = tuple(map(powers.__getitem__, ind_rows)) + # Echelonize the matrix ourselves, because otherwise the + # call to P.pivot_rows() below can choose a non-optimal + # row-reduction algorithm. In particular, scaling can + # 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 + if superalgebra.base_ring() is not QQ: + algo = "scaled_partial_pivoting" + P.echelonize(algorithm=algo) + + # Figure out which powers form a linearly-independent set. + ind_rows = P.pivot_rows() + + # Pick those out of the list of all powers. + basis = tuple(map(powers.__getitem__, ind_rows)) super().__init__(superalgebra, -- 2.43.2