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.
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 (
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``::
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):
# 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:
# 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
+
+ self._multiplication_table = tuple(map(tuple, self._multiplication_table))
+
+ # 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': True}
+ self._inner_product_matrix.set_immutable()
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():
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"""
+----++----+----+----+----+
"""
- 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)
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,
n = len(basis)
vector_basis = basis
- from sage.matrix.matrix import is_Matrix
+ from sage.structure.element import is_Matrix
basis_is_matrices = False
degree = 0
V = VectorSpace(field, degree)
- self._deorthonormalization_matrix = matrix.identity(field,n)
+ # 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:
- A = matrix(field, vector_basis)
- # uh oh, this is only the "usual" inner product
- Q,R = A.gram_schmidt(orthonormal=True)
- self._deorthonormalization_matrix = R.inverse().transpose()
- vector_basis = Q.rows()
+ # Compute the deorthonormalized tables before we orthonormalize
+ # the given basis.
W = V.span_of_basis( vector_basis )
+
+ # TODO: use symmetry
+ self._deortho_multiplication_table = [ [0 for j in range(n)]
+ for i in range(n) ]
+ self._deortho_inner_product_table = [ [0 for j in range(n)]
+ for i 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):
+ # given 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)
+
+ # TODO: use symmetry
+ elt = W.coordinate_vector(elt)
+ self._deortho_multiplication_table[i][j] = elt
+ self._deortho_multiplication_table[j][i] = elt
+ self._deortho_inner_product_table[i][j] = ip
+ self._deortho_inner_product_table[j][i] = ip
+
+ if self._deortho_multiplication_table is not None:
+ self._deortho_multiplication_table = tuple(map(tuple, self._deortho_multiplication_table))
+ if self._deortho_inner_product_table is not None:
+ self._deortho_inner_product_table = tuple(map(tuple, self._deortho_inner_product_table))
+
+ if orthonormalize:
+ from mjo.eja.eja_utils import gram_schmidt
+ 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,vector_basis) )
+ basis = tuple( map(_vec2mat,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) ]
+ # 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
+ # 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):
- # do another mat2vec because the multiplication
- # table is in terms of vectors
- elt = _mat2vec(jordan_product(basis[i],basis[j]))
+ # 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)
+
+ # TODO: use symmetry
elt = W.coordinate_vector(elt)
mult_table[i][j] = elt
mult_table[j][i] = elt
- ip = inner_product(basis[i],basis[j])
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,
+ ip_table,
prefix,
category,
- basis,
+ basis, # matrix basis
check_field,
check_axioms)
-class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
- r"""
- Algebras whose basis consists of vectors with rational
- entries. Equivalently, algebras whose multiplication tables
- contain only rational coefficients.
-
- When an EJA has a basis that can be made rational, we can speed up
- the computation of its characteristic polynomial by doing it over
- ``QQ``. All of the named EJA constructors that we provide fall
- into this category.
- """
@cached_method
def _charpoly_coefficients(self):
r"""
- Override the parent method with something that tries to compute
- over a faster (non-extension) field.
-
SETUP::
- sage: from mjo.eja.eja_algebra import JordanSpinEJA
+ sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
+ ....: JordanSpinEJA)
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
# rationals if this one is already over the rationals.
- superclass = super(RationalBasisEuclideanJordanAlgebra, self)
+ superclass = super(RationalBasisEuclideanJordanAlgebraNg, self)
return superclass._charpoly_coefficients()
- mult_table = tuple(
- tuple(map(lambda x: x.to_vector(), ls))
- for ls in self._multiplication_table
- )
-
# Do the computation over the rationals. The answer will be
- # the same, because our basis coordinates are (essentially)
- # rational.
+ # the same, because all we've done is a change of basis.
J = FiniteDimensionalEuclideanJordanAlgebra(QQ,
- mult_table,
- check_field=False,
- check_axioms=False)
+ 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:
r"""
A class for the Euclidean Jordan algebras that we know by name.
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):
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)
return cls(n, field, **kwargs)
-class HadamardEJA(RationalBasisEuclideanJordanAlgebra,
+class HadamardEJA(RationalBasisEuclideanJordanAlgebraNg,
ConcreteEuclideanJordanAlgebra):
"""
Return the Euclidean Jordan Algebra corresponding to the set
"""
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)
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)``
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)
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)
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
"""
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
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())
projects / sage.d.git / blobdiff