from sage.matrix.constructor import matrix
from sage.matrix.matrix_space import MatrixSpace
from sage.misc.cachefunc import cached_method
-from sage.misc.lazy_import import lazy_import
from sage.misc.table import table
from sage.modules.free_module import FreeModule, VectorSpace
from sage.rings.all import (ZZ, QQ, AA, QQbar, RR, RLF, CLF,
PolynomialRing,
QuadraticField)
from mjo.eja.eja_element import FiniteDimensionalEuclideanJordanAlgebraElement
-lazy_import('mjo.eja.eja_subalgebra',
- 'FiniteDimensionalEuclideanJordanSubalgebra')
from mjo.eja.eja_operator import FiniteDimensionalEuclideanJordanAlgebraOperator
from mjo.eja.eja_utils import _mat2vec
return (t**r + sum( a[k]*(t**k) for k in range(r) ))
+ def coordinate_polynomial_ring(self):
+ r"""
+ The multivariate polynomial ring in which this algebra's
+ :meth:`characteristic_polynomial_of` lives.
+
+ SETUP::
+
+ sage: from mjo.eja.eja_algebra import (HadamardEJA,
+ ....: RealSymmetricEJA)
+
+ EXAMPLES::
+
+ sage: J = HadamardEJA(2)
+ sage: J.coordinate_polynomial_ring()
+ Multivariate Polynomial Ring in X1, X2...
+ sage: J = RealSymmetricEJA(3,QQ)
+ sage: J.coordinate_polynomial_ring()
+ Multivariate Polynomial Ring in X1, X2, X3, X4, X5, X6...
+
+ """
+ var_names = tuple( "X%d" % z for z in range(1, self.dimension()+1) )
+ return PolynomialRing(self.base_ring(), var_names)
def inner_product(self, x, y):
"""
SETUP::
- sage: from mjo.eja.eja_algebra import random_eja
+ sage: from mjo.eja.eja_algebra import (random_eja,
+ ....: HadamardEJA,
+ ....: BilinearFormEJA)
EXAMPLES:
sage: (x*y).inner_product(z) == y.inner_product(x*z)
True
+ TESTS:
+
+ Ensure that this is the usual inner product for the algebras
+ over `R^n`::
+
+ sage: set_random_seed()
+ sage: J = HadamardEJA.random_instance()
+ sage: x,y = J.random_elements(2)
+ sage: actual = x.inner_product(y)
+ sage: expected = x.to_vector().inner_product(y.to_vector())
+ sage: actual == expected
+ True
+
+ Ensure that this is one-half of the trace inner-product in a
+ BilinearFormEJA that isn't just the reals (when ``n`` isn't
+ one). This is in Faraut and Koranyi, and also my "On the
+ symmetry..." paper::
+
+ sage: set_random_seed()
+ sage: J = BilinearFormEJA.random_instance()
+ sage: n = J.dimension()
+ sage: x = J.random_element()
+ sage: y = J.random_element()
+ sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
+ True
"""
- X = x.natural_representation()
- Y = y.natural_representation()
- return self.natural_inner_product(X,Y)
+ B = self._inner_product_matrix
+ return (B*x.to_vector()).inner_product(y.to_vector())
def is_trivial(self):
return self._natural_basis[0].matrix_space()
- @staticmethod
- def natural_inner_product(X,Y):
- """
- Compute the inner product of two naturally-represented elements.
-
- For example in the real symmetric matrix EJA, this will compute
- the trace inner-product of two n-by-n symmetric matrices. The
- default should work for the real cartesian product EJA, the
- Jordan spin EJA, and the real symmetric matrices. The others
- will have to be overridden.
- """
- return (X.conjugate_transpose()*Y).trace()
-
-
@cached_method
def one(self):
"""
if not c.is_idempotent():
raise ValueError("element is not idempotent: %s" % c)
+ from mjo.eja.eja_subalgebra import FiniteDimensionalEuclideanJordanSubalgebra
+
# Default these to what they should be if they turn out to be
# trivial, because eigenspaces_left() won't return eigenvalues
# corresponding to trivial spaces (e.g. it returns only the
of" function.
"""
n = self.dimension()
- var_names = [ "X" + str(z) for z in range(1,n+1) ]
- R = PolynomialRing(self.base_ring(), var_names)
+ R = self.coordinate_polynomial_ring()
vars = R.gens()
F = R.fraction_field()
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
~(self.natural_inner_product(s,s).sqrt()) for s in basis )
basis = tuple(s*c for (s,c) in zip(basis,self._basis_normalizers))
- Qs = self.multiplication_table_from_matrix_basis(basis)
+ # 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 = [[W.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.natural_inner_product(basis[i],
+ basis[j])
+ except:
+ pass
+
+ try:
+ # HACK PART DEUX
+ self._inner_product_matrix = matrix(field,ip_table)
+ except:
+ pass
super(MatrixEuclideanJordanAlgebra, self).__init__(field,
- Qs,
+ mult_table,
natural_basis=basis,
**kwargs)
return tuple( a_i.subs(subs_dict) for a_i in a )
- @staticmethod
- def multiplication_table_from_matrix_basis(basis):
- """
- At least three of the five simple Euclidean Jordan algebras have the
- symmetric multiplication (A,B) |-> (AB + BA)/2, where the
- multiplication on the right is matrix multiplication. Given a basis
- for the underlying matrix space, this function returns a
- multiplication table (obtained by looping through the basis
- elements) for an algebra of those matrices.
- """
- # In S^2, for example, we nominally have four coordinates even
- # though the space is of dimension three only. The vector space V
- # is supposed to hold the entire long vector, and the subspace W
- # of V will be spanned by the vectors that arise from symmetric
- # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
- if len(basis) == 0:
- return []
-
- field = basis[0].base_ring()
- dimension = basis[0].nrows()
-
- V = VectorSpace(field, dimension**2)
- W = V.span_of_basis( _mat2vec(s) for s in basis )
- n = len(basis)
- mult_table = [[W.zero() for j in range(n)] for i in range(n)]
- for i in range(n):
- for j in range(n):
- mat_entry = (basis[i]*basis[j] + basis[j]*basis[i])/2
- mult_table[i][j] = W.coordinate_vector(_mat2vec(mat_entry))
-
- return mult_table
-
-
@staticmethod
def real_embed(M):
"""
"""
raise NotImplementedError
-
@classmethod
def natural_inner_product(cls,X,Y):
Xu = cls.real_unembed(X)
mult_table = [ [ V.gen(i)*(i == j) for j in range(n) ]
for i in range(n) ]
+ # 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
+
super(HadamardEJA, self).__init__(field,
mult_table,
check_axioms=False,
return cls(n, field, **kwargs)
- def inner_product(self, x, y):
- """
- Faster to reimplement than to use natural representations.
-
- SETUP::
-
- sage: from mjo.eja.eja_algebra import HadamardEJA
-
- TESTS:
-
- Ensure that this is the usual inner product for the algebras
- over `R^n`::
-
- sage: set_random_seed()
- sage: J = HadamardEJA.random_instance()
- sage: x,y = J.random_elements(2)
- sage: X = x.natural_representation()
- sage: Y = y.natural_representation()
- sage: x.inner_product(y) == J.natural_inner_product(X,Y)
- True
-
- """
- return x.to_vector().inner_product(y.to_vector())
-
-
class BilinearFormEJA(RationalBasisEuclideanJordanAlgebra,
ConcreteEuclideanJordanAlgebra):
r"""
True
"""
def __init__(self, B, field=AA, **kwargs):
- self._B = B
n = B.nrows()
if not B.is_positive_definite():
z = V([z0] + zbar.list())
mult_table[i][j] = z
- # The rank of this algebra is two, unless we're in a
- # one-dimensional ambient space (because the rank is bounded
- # by the ambient dimension).
+ # 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
+
super(BilinearFormEJA, self).__init__(field,
mult_table,
check_axioms=False,
**kwargs)
+
+ # The rank of this algebra is two, unless we're in a
+ # one-dimensional ambient space (because the rank is bounded
+ # by the ambient dimension).
self.rank.set_cache(min(n,2))
if n == 0:
return cls(B, field, **kwargs)
- def inner_product(self, x, y):
- r"""
- Half of the trace inner product.
-
- This is defined so that the special case of the Jordan spin
- algebra gets the usual inner product.
-
- SETUP::
-
- sage: from mjo.eja.eja_algebra import BilinearFormEJA
-
- TESTS:
-
- Ensure that this is one-half of the trace inner-product when
- the algebra isn't just the reals (when ``n`` isn't one). This
- is in Faraut and Koranyi, and also my "On the symmetry..."
- paper::
-
- sage: set_random_seed()
- sage: J = BilinearFormEJA.random_instance()
- sage: n = J.dimension()
- sage: x = J.random_element()
- sage: y = J.random_element()
- sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
- True
-
- """
- return (self._B*x.to_vector()).inner_product(y.to_vector())
-
class JordanSpinEJA(BilinearFormEJA):
"""
sage: set_random_seed()
sage: J = JordanSpinEJA.random_instance()
sage: x,y = J.random_elements(2)
- sage: X = x.natural_representation()
- sage: Y = y.natural_representation()
- sage: x.inner_product(y) == J.natural_inner_product(X,Y)
+ sage: actual = x.inner_product(y)
+ sage: expected = x.to_vector().inner_product(y.to_vector())
+ sage: actual == expected
True
"""
"""
def __init__(self, field=AA, **kwargs):
mult_table = []
+ self._inner_product_matrix = matrix(field,0)
super(TrivialEJA, self).__init__(field,
mult_table,
check_axioms=False,
"""
(J1,J2) = self.factors()
- n = J1.dimension()
+ m = J1.dimension()
+ n = J2.dimension()
V_basis = self.vector_space().basis()
- P1 = matrix(self.base_ring(), V_basis[:n])
- P2 = matrix(self.base_ring(), V_basis[n:])
+ # Need to specify the dimensions explicitly so that we don't
+ # wind up with a zero-by-zero matrix when we want e.g. a
+ # zero-by-two matrix (important for composing things).
+ P1 = matrix(self.base_ring(), m, m+n, V_basis[:m])
+ P2 = matrix(self.base_ring(), n, m+n, V_basis[m:])
pi_left = FiniteDimensionalEuclideanJordanAlgebraOperator(self,J1,P1)
pi_right = FiniteDimensionalEuclideanJordanAlgebraOperator(self,J2,P2)
return (pi_left, pi_right)
"""
(J1,J2) = self.factors()
- n = J1.dimension()
+ m = J1.dimension()
+ n = J2.dimension()
V_basis = self.vector_space().basis()
- I1 = matrix.column(self.base_ring(), V_basis[:n])
- I2 = matrix.column(self.base_ring(), V_basis[n:])
+ # Need to specify the dimensions explicitly so that we don't
+ # wind up with a zero-by-zero matrix when we want e.g. a
+ # two-by-zero matrix (important for composing things).
+ I1 = matrix.column(self.base_ring(), m, m+n, V_basis[:m])
+ I2 = matrix.column(self.base_ring(), n, m+n, V_basis[m:])
iota_left = FiniteDimensionalEuclideanJordanAlgebraOperator(J1,self,I1)
iota_right = FiniteDimensionalEuclideanJordanAlgebraOperator(J2,self,I2)
return (iota_left, iota_right)
projects / sage.d.git / blobdiff