sage: J = HadamardEJA(2)
sage: J.coordinate_polynomial_ring()
Multivariate Polynomial Ring in X1, X2...
- sage: J = RealSymmetricEJA(3,QQ)
+ sage: J = RealSymmetricEJA(3,QQ,orthonormalize=False)
sage: J.coordinate_polynomial_ring()
Multivariate Polynomial Ring in X1, X2, X3, X4, X5, X6...
if n > 0:
if is_Matrix(basis[0]):
basis_is_matrices = True
+ from mjo.eja.eja_utils import _vec2mat
vector_basis = tuple( map(_mat2vec,basis) )
degree = basis[0].nrows()**2
else:
if self._deortho_inner_product_table is not None:
self._deortho_inner_product_table = tuple(map(tuple, self._deortho_inner_product_table))
+ # 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
+ self._deortho_matrix = None
+
if orthonormalize:
from mjo.eja.eja_utils import gram_schmidt
- vector_basis = gram_schmidt(vector_basis, inner_product)
+ 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)
+
W = V.span_of_basis( vector_basis )
# Normalize the "matrix" basis, too!
basis = vector_basis
if basis_is_matrices:
- from mjo.eja.eja_utils import _vec2mat
basis = tuple( map(_vec2mat,basis) )
W = V.span_of_basis( vector_basis )
+ # Now "W" is the vector space of our algebra coordinates. The
+ # variables "X1", "X2",... refer to the entries of vectors in
+ # W. Thus to convert back and forth between the orthonormal
+ # coordinates and the given ones, we need to stick the original
+ # basis in W.
+ U = V.span_of_basis( deortho_vector_basis )
+ self._deortho_matrix = matrix( U.coordinate_vector(q)
+ for q in vector_basis )
+
# TODO: use symmetry
mult_table = [ [0 for j in range(n)] for i in range(n) ]
ip_table = [ [0 for j in range(n)] for i in range(n) ]
- # Note: the Jordan and inner- products are defined in terms
+ # 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):
EXAMPLES:
- The returned coefficients should be the same as if we'd never
- orthonormalized the basis to begin with::
-
- sage: B = matrix(QQ, [[1, 0, 0],
- ....: [0, 25, -32],
- ....: [0, -32, 41] ])
- sage: J1 = BilinearFormEJA(B)
- sage: J2 = BilinearFormEJA(B,QQ,orthonormalize=False)
- sage: J1._charpoly_coefficients()
- (X1^2 - 25*X2^2 + 64*X2*X3 - 41*X3^2, -2*X1)
- sage: J2._charpoly_coefficients()
- (X1^2 - 25*X2^2 + 64*X2*X3 - 41*X3^2, -2*X1)
-
The base ring of the resulting polynomial coefficients is what
it should be, and not the rationals (unless the algebra was
already over the rationals)::
Algebraic Real Field
sage: a0.base_ring()
Algebraic Real Field
+
"""
if self.base_ring() is QQ:
# There's no need to construct *another* algebra over the
J = FiniteDimensionalEuclideanJordanAlgebra(QQ,
self._deortho_multiplication_table,
self._deortho_inner_product_table)
- a = J._charpoly_coefficients()
- return tuple(map(lambda x: x.change_ring(self.base_ring()), a))
-class ConcreteEuclideanJordanAlgebra:
+ # Change back from QQ to our real base ring
+ a = ( a_i.change_ring(self.base_ring())
+ for a_i in J._charpoly_coefficients() )
+
+ # Now convert the coordinate variables back to the
+ # deorthonormalized ones.
+ R = self.coordinate_polynomial_ring()
+ from sage.modules.free_module_element import vector
+ X = vector(R, R.gens())
+ BX = self._deortho_matrix*X
+
+ subs_dict = { X[i]: BX[i] for i in range(len(X)) }
+ return tuple( a_i.subs(subs_dict) for a_i in a )
+
+class ConcreteEuclideanJordanAlgebra(RationalBasisEuclideanJordanAlgebra):
r"""
A class for the Euclidean Jordan algebras that we know by name.
raise NotImplementedError
@classmethod
- def random_instance(cls, field=AA, **kwargs):
+ def random_instance(cls, *args, **kwargs):
"""
Return a random instance of this type of algebra.
"""
from sage.misc.prandom import choice
eja_class = choice(cls.__subclasses__())
- return eja_class.random_instance(field)
-
-
-class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
-
- def __init__(self, field, basis, normalize_basis=True, **kwargs):
- """
- Compared to the superclass constructor, we take a basis instead of
- 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_coefficients() override.
- self._basis_normalizers = None
-
- # We're going to loop through this a few times, so now's a good
- # time to ensure that it isn't a generator expression.
- basis = tuple(basis)
-
- algebra_dim = len(basis)
- degree = 0 # size of the matrices
- if algebra_dim > 0:
- degree = basis[0].nrows()
-
- if algebra_dim > 1 and normalize_basis:
- # We'll need sqrt(2) to normalize the basis, and this
- # winds up in the multiplication table, so the whole
- # algebra needs to be over the field extension.
- R = PolynomialRing(field, 'z')
- z = R.gen()
- p = z**2 - 2
- if p.is_irreducible():
- field = field.extension(p, 'sqrt2', embedding=RLF(2).sqrt())
- basis = tuple( s.change_ring(field) for s in basis )
- self._basis_normalizers = tuple(
- ~(self.matrix_inner_product(s,s).sqrt()) for s in basis )
- basis = tuple(s*c for (s,c) in zip(basis,self._basis_normalizers))
-
- # Now compute the multiplication and inner product tables.
- # We have to do this *after* normalizing the basis, because
- # scaling affects the answers.
- V = VectorSpace(field, degree**2)
- 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 = [[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):
- mat_entry = (basis[i]*basis[j] + basis[j]*basis[i])/2
- mult_table[i][j] = W.coordinate_vector(_mat2vec(mat_entry))
-
- try:
- # HACK: ignore the error here if we don't need the
- # inner product (as is the case when we construct
- # a dummy QQ-algebra for fast charpoly coefficients.
- ip_table[i][j] = self.matrix_inner_product(basis[i],
- basis[j])
- except:
- pass
-
- super(MatrixEuclideanJordanAlgebra, self).__init__(field,
- mult_table,
- ip_table,
- matrix_basis=basis,
- **kwargs)
-
- if algebra_dim == 0:
- self.one.set_cache(self.zero())
- else:
- n = basis[0].nrows()
- # The identity wrt (A,B) -> (AB + BA)/2 is independent of the
- # details of this algebra.
- self.one.set_cache(self(matrix.identity(field,n)))
-
- @cached_method
- def _charpoly_coefficients(self):
- r"""
- Override the parent method with something that tries to compute
- over a faster (non-extension) field.
- """
- if self._basis_normalizers is None or self.base_ring() is QQ:
- # We didn't normalize, or the basis we started with had
- # entries in a nice field already. Just compute the thing.
- return super(MatrixEuclideanJordanAlgebra, self)._charpoly_coefficients()
-
- basis = ( (b/n) for (b,n) in zip(self.matrix_basis(),
- self._basis_normalizers) )
-
- # Do this over the rationals and convert back at the end.
- # Only works because we know the entries of the basis are
- # integers. The argument ``check_axioms=False`` is required
- # because the trace inner-product method for this
- # class is a stub and can't actually be checked.
- J = MatrixEuclideanJordanAlgebra(QQ,
- basis,
- normalize_basis=False,
- check_field=False,
- check_axioms=False)
- 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 )
+ # These all bubble up to the RationalBasisEuclideanJordanAlgebra
+ # superclass constructor, so any (kw)args valid there are also
+ # valid here.
+ return eja_class.random_instance(*args, **kwargs)
+class MatrixEuclideanJordanAlgebra:
@staticmethod
def real_embed(M):
"""
"""
raise NotImplementedError
+ @staticmethod
+ def jordan_product(X,Y):
+ return (X*Y + Y*X)/2
+
@classmethod
- def matrix_inner_product(cls,X,Y):
+ def trace_inner_product(cls,X,Y):
Xu = cls.real_unembed(X)
Yu = cls.real_unembed(Y)
tr = (Xu*Yu).trace()
return M
-class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra,
- ConcreteEuclideanJordanAlgebra):
+class RealSymmetricEJA(ConcreteEuclideanJordanAlgebra,
+ RealMatrixEuclideanJordanAlgebra):
"""
The rank-n simple EJA consisting of real symmetric n-by-n
matrices, the usual symmetric Jordan product, and the trace inner
basis = self._denormalized_basis(n, field)
super(RealSymmetricEJA, self).__init__(field,
basis,
- check_axioms=False,
+ self.jordan_product,
+ self.trace_inner_product,
**kwargs)
self.rank.set_cache(n)
+ self.one.set_cache(self(matrix.identity(field,n)))
class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
@classmethod
- def matrix_inner_product(cls,X,Y):
+ def trace_inner_product(cls,X,Y):
"""
Compute a matrix inner product in this algebra directly from
its real embedding.
sage: X = ComplexHermitianEJA.real_unembed(Xe)
sage: Y = ComplexHermitianEJA.real_unembed(Ye)
sage: expected = (X*Y).trace().real()
- sage: actual = ComplexHermitianEJA.matrix_inner_product(Xe,Ye)
+ sage: actual = ComplexHermitianEJA.trace_inner_product(Xe,Ye)
sage: actual == expected
True
"""
- return RealMatrixEuclideanJordanAlgebra.matrix_inner_product(X,Y)/2
+ return RealMatrixEuclideanJordanAlgebra.trace_inner_product(X,Y)/2
-class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra,
- ConcreteEuclideanJordanAlgebra):
+class ComplexHermitianEJA(ConcreteEuclideanJordanAlgebra,
+ ComplexMatrixEuclideanJordanAlgebra):
"""
The rank-n simple EJA consisting of complex Hermitian n-by-n
matrices over the real numbers, the usual symmetric Jordan product,
# Since we embedded these, we can drop back to the "field" that we
# started with instead of the complex extension "F".
- return ( s.change_ring(field) for s in S )
+ return tuple( s.change_ring(field) for s in S )
def __init__(self, n, field=AA, **kwargs):
basis = self._denormalized_basis(n,field)
- super(ComplexHermitianEJA,self).__init__(field,
- basis,
- check_axioms=False,
- **kwargs)
+ super(ComplexHermitianEJA, self).__init__(field,
+ basis,
+ self.jordan_product,
+ self.trace_inner_product,
+ **kwargs)
self.rank.set_cache(n)
+ # TODO: pre-cache the identity!
@staticmethod
def _max_random_instance_size():
@classmethod
- def matrix_inner_product(cls,X,Y):
+ def trace_inner_product(cls,X,Y):
"""
Compute a matrix inner product in this algebra directly from
its real embedding.
sage: X = QuaternionHermitianEJA.real_unembed(Xe)
sage: Y = QuaternionHermitianEJA.real_unembed(Ye)
sage: expected = (X*Y).trace().coefficient_tuple()[0]
- sage: actual = QuaternionHermitianEJA.matrix_inner_product(Xe,Ye)
+ sage: actual = QuaternionHermitianEJA.trace_inner_product(Xe,Ye)
sage: actual == expected
True
"""
- return RealMatrixEuclideanJordanAlgebra.matrix_inner_product(X,Y)/4
+ return RealMatrixEuclideanJordanAlgebra.trace_inner_product(X,Y)/4
-class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra,
- ConcreteEuclideanJordanAlgebra):
+class QuaternionHermitianEJA(ConcreteEuclideanJordanAlgebra,
+ QuaternionMatrixEuclideanJordanAlgebra):
r"""
The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
matrices, the usual symmetric Jordan product, and the
# Since we embedded these, we can drop back to the "field" that we
# started with instead of the quaternion algebra "Q".
- return ( s.change_ring(field) for s in S )
+ return tuple( s.change_ring(field) for s in S )
def __init__(self, n, field=AA, **kwargs):
basis = self._denormalized_basis(n,field)
- super(QuaternionHermitianEJA,self).__init__(field,
- basis,
- check_axioms=False,
- **kwargs)
+ super(QuaternionHermitianEJA, self).__init__(field,
+ basis,
+ self.jordan_product,
+ self.trace_inner_product,
+ **kwargs)
self.rank.set_cache(n)
+ # TODO: cache one()!
@staticmethod
def _max_random_instance_size():
return cls(n, field, **kwargs)
-class HadamardEJA(RationalBasisEuclideanJordanAlgebra,
- ConcreteEuclideanJordanAlgebra):
+class HadamardEJA(ConcreteEuclideanJordanAlgebra):
"""
Return the Euclidean Jordan Algebra corresponding to the set
`R^n` under the Hadamard product.
return cls(n, field, **kwargs)
-class BilinearFormEJA(RationalBasisEuclideanJordanAlgebra,
- ConcreteEuclideanJordanAlgebra):
+class BilinearFormEJA(ConcreteEuclideanJordanAlgebra):
r"""
The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
with the half-trace inner product and jordan product ``x*y =
return cls(n, field, **kwargs)
-class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra,
- ConcreteEuclideanJordanAlgebra):
+class TrivialEJA(ConcreteEuclideanJordanAlgebra):
"""
The trivial Euclidean Jordan algebra consisting of only a zero element.
"""
def __init__(self, field=AA, **kwargs):
- mult_table = []
- ip_table = []
+ jordan_product = lambda x,y: x
+ inner_product = lambda x,y: field(0)
+ basis = ()
super(TrivialEJA, self).__init__(field,
- mult_table,
- ip_table,
- check_axioms=False,
+ basis,
+ jordan_product,
+ inner_product,
**kwargs)
# The rank is zero using my definition, namely the dimension of the
# largest subalgebra generated by any element.
EXAMPLE::
sage: J1 = HadamardEJA(3,QQ)
- sage: J2 = QuaternionHermitianEJA(2,QQ,normalize_basis=False)
+ sage: J2 = QuaternionHermitianEJA(2,QQ,orthonormalize=False)
sage: J = DirectSumEJA(J1,J2)
sage: x1 = J1.one()
sage: x2 = x1
projects / sage.d.git / commitdiff