X-Git-Url: https://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=mjo%2Feja%2Feja_algebra.py;h=9eba6699819956048747f648fd18bc80847fbbb6;hb=c9db5fa5d2f30833e4da85d32418048d34483e43;hp=adcf75503d2a55add4c32dc56380396abd1f8d9d;hpb=f50d52bf5024f319c8757ece8ad95edd0b49270d;p=sage.d.git diff --git a/mjo/eja/eja_algebra.py b/mjo/eja/eja_algebra.py index adcf755..9eba669 100644 --- a/mjo/eja/eja_algebra.py +++ b/mjo/eja/eja_algebra.py @@ -1025,6 +1025,103 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): Element = FiniteDimensionalEuclideanJordanAlgebraElement +class RationalBasisEuclideanJordanAlgebraNg(FiniteDimensionalEuclideanJordanAlgebra): + def __init__(self, + field, + basis, + jordan_product, + inner_product, + orthonormalize=True, + prefix='e', + category=None, + check_field=True, + check_axioms=True): + + n = len(basis) + vector_basis = basis + + from sage.structure.element import is_Matrix + basis_is_matrices = False + + degree = 0 + if n > 0: + if is_Matrix(basis[0]): + basis_is_matrices = True + vector_basis = tuple( map(_mat2vec,basis) ) + degree = basis[0].nrows()**2 + else: + degree = basis[0].degree() + + 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 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) + + if basis_is_matrices: + from mjo.eja.eja_utils import _vec2mat + basis = tuple( map(_vec2mat,vector_basis) ) + + W = V.span_of_basis( vector_basis ) + + mult_table = [ [0 for i in range(n)] for j in range(n) ] + ip_table = [ [0 for i in range(n)] for j in range(n) ] + + # 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) + + 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) + + elt = W.coordinate_vector(elt) + mult_table[i][j] = elt + mult_table[j][i] = elt + 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() + else: + basis = tuple( x.column() for x in basis ) + + super().__init__(field, + mult_table, + prefix, + category, + basis, # matrix basis + check_field, + check_axioms) class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra): r""" @@ -2040,7 +2137,7 @@ class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra, return cls(n, field, **kwargs) -class HadamardEJA(RationalBasisEuclideanJordanAlgebra, +class HadamardEJA(RationalBasisEuclideanJordanAlgebraNg, ConcreteEuclideanJordanAlgebra): """ Return the Euclidean Jordan Algebra corresponding to the set @@ -2083,22 +2180,17 @@ class HadamardEJA(RationalBasisEuclideanJordanAlgebra, """ def __init__(self, n, field=AA, **kwargs): V = VectorSpace(field, n) - mult_table = [ [ V.gen(i)*(i == j) for j in range(n) ] - for i in range(n) ] + basis = V.basis() - # Inner products are real numbers and not algebra - # elements, so once we turn the algebra element - # into a vector in inner_product(), we never go - # back. As a result -- contrary to what we do with - # self._multiplication_table -- we store the inner - # product table as a plain old matrix and not as - # an algebra operator. - ip_table = matrix.identity(field,n) - self._inner_product_matrix = ip_table + def jordan_product(x,y): + return V([ xi*yi for (xi,yi) in zip(x,y) ]) + def inner_product(x,y): + return x.inner_product(y) super(HadamardEJA, self).__init__(field, - mult_table, - check_axioms=False, + basis, + jordan_product, + inner_product, **kwargs) self.rank.set_cache(n) @@ -2123,7 +2215,7 @@ class HadamardEJA(RationalBasisEuclideanJordanAlgebra, return cls(n, field, **kwargs) -class BilinearFormEJA(RationalBasisEuclideanJordanAlgebra, +class BilinearFormEJA(RationalBasisEuclideanJordanAlgebraNg, ConcreteEuclideanJordanAlgebra): r""" The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)`` @@ -2204,39 +2296,28 @@ class BilinearFormEJA(RationalBasisEuclideanJordanAlgebra, True """ def __init__(self, B, field=AA, **kwargs): - n = B.nrows() - if not B.is_positive_definite(): raise ValueError("bilinear form is not positive-definite") + n = B.nrows() V = VectorSpace(field, n) - mult_table = [[V.zero() for j in range(n)] for i in range(n)] - for i in range(n): - for j in range(n): - x = V.gen(i) - y = V.gen(j) - x0 = x[0] - xbar = x[1:] - y0 = y[0] - ybar = y[1:] - z0 = (B*x).inner_product(y) - zbar = y0*xbar + x0*ybar - z = V([z0] + zbar.list()) - mult_table[i][j] = z - - # Inner products are real numbers and not algebra - # elements, so once we turn the algebra element - # into a vector in inner_product(), we never go - # back. As a result -- contrary to what we do with - # self._multiplication_table -- we store the inner - # product table as a plain old matrix and not as - # an algebra operator. - ip_table = B - self._inner_product_matrix = ip_table + + def inner_product(x,y): + return (B*x).inner_product(y) + + def jordan_product(x,y): + x0 = x[0] + xbar = x[1:] + y0 = y[0] + ybar = y[1:] + z0 = inner_product(x,y) + zbar = y0*xbar + x0*ybar + return V([z0] + zbar.list()) super(BilinearFormEJA, self).__init__(field, - mult_table, - check_axioms=False, + V.basis(), + jordan_product, + inner_product, **kwargs) # The rank of this algebra is two, unless we're in a