specifically those where u^2 + v^2 = 0 implies that u = v = 0. They
are used in optimization, and have some additional nice methods beyond
what can be supported in a general Jordan Algebra.
+
+
+SETUP::
+
+ sage: from mjo.eja.eja_algebra import random_eja
+
+EXAMPLES::
+
+ sage: random_eja()
+ Euclidean Jordan algebra of dimension...
+
"""
from itertools import repeat
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.prandom import choice
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
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.
sage: J(1)
Traceback (most recent call last):
...
- ValueError: not a naturally-represented algebra element
+ ValueError: not an element of this algebra
"""
return None
def __init__(self,
field,
- mult_table,
+ multiplication_table,
+ inner_product_table,
prefix='e',
category=None,
- natural_basis=None,
- check=True):
+ 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 (JordanSpinEJA, random_eja)
+ sage: from mjo.eja.eja_algebra import (
+ ....: FiniteDimensionalEuclideanJordanAlgebra,
+ ....: JordanSpinEJA,
+ ....: random_eja)
EXAMPLES:
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::
+ The ``field`` we're given must be real with ``check_field=True``::
sage: JordanSpinEJA(2,QQbar)
Traceback (most recent call last):
...
- ValueError: field is not real
+ ValueError: scalar field is not real
+
+ The multiplication table must be square with ``check_axioms=True``::
+
+ 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:
+ if check_field:
if not field.is_subring(RR):
# Note: this does return true for the real algebraic
- # field, and any quadratic field where we've specified
- # a real embedding.
- raise ValueError('field is not real')
-
- self._natural_basis = natural_basis
+ # field, the rationals, and any quadratic field where
+ # we've specified a real embedding.
+ raise ValueError("scalar field is not real")
+
+
+ # 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 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:
category = MagmaticAlgebras(field).FiniteDimensional()
category = category.WithBasis().Unital()
- # The multiplication table had better be square
- n = len(mult_table)
-
fda = super(FiniteDimensionalEuclideanJordanAlgebra, self)
fda.__init__(field,
range(n),
# 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.
- self._multiplication_table = [
- list(map(lambda x: self.from_vector(x), ls))
- for ls in mult_table
- ]
-
- if check:
- if not self._is_commutative():
- raise ValueError("algebra is not commutative")
+ #
+ # 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 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(multiplication_table[i][j])
+ self._multiplication_table[i][j] = elt
+
+ # 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': False}
+
+ if check_axioms:
if not self._is_jordanian():
raise ValueError("Jordan identity does not hold")
if not self._inner_product_is_associative():
def _element_constructor_(self, elt):
"""
- Construct an element of this algebra from its natural
+ Construct an element of this algebra from its vector or matrix
representation.
This gets called only after the parent element _call_ method
sage: J(A)
Traceback (most recent call last):
...
- ArithmeticError: vector is not in free module
+ ValueError: not an element of this algebra
TESTS:
Ensure that we can convert any element of the two non-matrix
- simple algebras (whose natural representations are their usual
- vector representations) back and forth faithfully::
+ simple algebras (whose matrix representations are columns)
+ back and forth faithfully::
sage: set_random_seed()
sage: J = HadamardEJA.random_instance()
sage: x = J.random_element()
sage: J(x.to_vector().column()) == x
True
-
"""
- msg = "not a naturally-represented algebra element"
+ msg = "not an element of this algebra"
if elt == 0:
# The superclass implementation of random_element()
# needs to be able to coerce "0" into the algebra.
# that the integer 3 belongs to the space of 2-by-2 matrices.
raise ValueError(msg)
- natural_basis = self.natural_basis()
- basis_space = natural_basis[0].matrix_space()
- if elt not in basis_space:
+ if elt not in self.matrix_space():
raise ValueError(msg)
# Thanks for nothing! Matrix spaces aren't vector spaces in
- # Sage, so we have to figure out its natural-basis coordinates
+ # Sage, so we have to figure out its matrix-basis coordinates
# ourselves. We use the basis space's ring instead of the
# element's ring because the basis space might be an algebraic
# closure whereas the base ring of the 3-by-3 identity matrix
# could be QQ instead of QQbar.
- V = VectorSpace(basis_space.base_ring(), elt.nrows()*elt.ncols())
- W = V.span_of_basis( _mat2vec(s) for s in natural_basis )
- coords = W.coordinate_vector(_mat2vec(elt))
- return self.from_vector(coords)
+ V = VectorSpace(self.base_ring(), elt.nrows()*elt.ncols())
+ W = V.span_of_basis( _mat2vec(s) for s in self.matrix_basis() )
- @staticmethod
- def _max_test_case_size():
- """
- Return an integer "size" that is an upper bound on the size of
- this algebra when it is used in a random test
- case. Unfortunately, the term "size" is quite vague -- when
- dealing with `R^n` under either the Hadamard or Jordan spin
- product, the "size" refers to the dimension `n`. When dealing
- with a matrix algebra (real symmetric or complex/quaternion
- Hermitian), it refers to the size of the matrix, which is
- far less than the dimension of the underlying vector space.
+ try:
+ coords = W.coordinate_vector(_mat2vec(elt))
+ except ArithmeticError: # vector is not in free module
+ raise ValueError(msg)
- We default to five in this class, which is safe in `R^n`. The
- matrix algebra subclasses (or any class where the "size" is
- interpreted to be far less than the dimension) should override
- with a smaller number.
- """
- return 5
+ return self.from_vector(coords)
def _repr_(self):
"""
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"""
Whether or not this algebra's multiplication table is commutative.
This method should of course always return ``True``, unless
- this algebra was constructed with ``check=False`` and passed
- an invalid multiplication table.
+ this algebra was constructed with ``check_axioms=False`` and
+ passed an invalid multiplication table.
"""
return all( self.product_on_basis(i,j) == self.product_on_basis(i,j)
for i in range(self.dimension())
``True`` result to be truly true, you should also check
:meth:`_is_commutative`. This method should of course always
return ``True``, unless this algebra was constructed with
- ``check=False`` and passed an invalid multiplication table.
+ ``check_axioms=False`` and passed an invalid multiplication table.
"""
return all( (self.monomial(i)**2)*(self.monomial(i)*self.monomial(j))
==
associative; that is, whether or not `B(xy,z) = B(x,yz)`.
This method should of course always return ``True``, unless
- this algebra was constructed with ``check=False`` and passed
- an invalid multiplication table.
+ this algebra was constructed with ``check_axioms=False`` and
+ passed an invalid multiplication table.
"""
+
+ # Used to check whether or not something is zero in an inexact
+ # ring. This number is sufficient to allow the construction of
+ # QuaternionHermitianEJA(2, RDF) with check_axioms=True.
+ epsilon = 1e-16
+
for i in range(self.dimension()):
for j in range(self.dimension()):
for k in range(self.dimension()):
x = self.monomial(i)
y = self.monomial(j)
z = self.monomial(k)
- if (x*y).inner_product(z) != x.inner_product(y*z):
- return False
+ diff = (x*y).inner_product(z) - x.inner_product(y*z)
+
+ if self.base_ring().is_exact():
+ if diff != 0:
+ return False
+ else:
+ if diff.abs() > epsilon:
+ return False
return True
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):
+----++----+----+----+----+
"""
- 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)
- def natural_basis(self):
+ def matrix_basis(self):
"""
- Return a more-natural representation of this algebra's basis.
+ Return an (often more natural) representation of this algebras
+ basis as an ordered tuple of matrices.
- Every finite-dimensional Euclidean Jordan Algebra is a direct
- sum of five simple algebras, four of which comprise Hermitian
- matrices. This method returns the original "natural" basis
- for our underlying vector space. (Typically, the natural basis
- is used to construct the multiplication table in the first place.)
+ Every finite-dimensional Euclidean Jordan Algebra is a, up to
+ Jordan isomorphism, a direct sum of five simple
+ algebras---four of which comprise Hermitian matrices. And the
+ last type of algebra can of course be thought of as `n`-by-`1`
+ column matrices (ambiguusly called column vectors) to avoid
+ special cases. As a result, matrices (and column vectors) are
+ a natural representation format for Euclidean Jordan algebra
+ elements.
- Note that this will always return a matrix. The standard basis
- in `R^n` will be returned as `n`-by-`1` column matrices.
+ But, when we construct an algebra from a basis of matrices,
+ those matrix representations are lost in favor of coordinate
+ vectors *with respect to* that basis. We could eventually
+ convert back if we tried hard enough, but having the original
+ representations handy is valuable enough that we simply store
+ them and return them from this method.
+
+ Why implement this for non-matrix algebras? Avoiding special
+ cases for the :class:`BilinearFormEJA` pays with simplicity in
+ its own right. But mainly, we would like to be able to assume
+ that elements of a :class:`DirectSumEJA` can be displayed
+ nicely, without having to have special classes for direct sums
+ one of whose components was a matrix algebra.
SETUP::
sage: J = RealSymmetricEJA(2)
sage: J.basis()
Finite family {0: e0, 1: e1, 2: e2}
- sage: J.natural_basis()
+ sage: J.matrix_basis()
(
[1 0] [ 0 0.7071067811865475?] [0 0]
[0 0], [0.7071067811865475? 0], [0 1]
sage: J = JordanSpinEJA(2)
sage: J.basis()
Finite family {0: e0, 1: e1}
- sage: J.natural_basis()
+ sage: J.matrix_basis()
(
[1] [0]
[0], [1]
)
-
"""
- if self._natural_basis is None:
- M = self.natural_basis_space()
+ if self._matrix_basis is None:
+ M = self.matrix_space()
return tuple( M(b.to_vector()) for b in self.basis() )
else:
- return self._natural_basis
+ return self._matrix_basis
- def natural_basis_space(self):
+ def matrix_space(self):
"""
- Return the matrix space in which this algebra's natural basis
- elements live.
+ Return the matrix space in which this algebra's elements live, if
+ we think of them as matrices (including column vectors of the
+ appropriate size).
Generally this will be an `n`-by-`1` column-vector space,
except when the algebra is trivial. There it's `n`-by-`n`
- (where `n` is zero), to ensure that two elements of the
- natural basis space (empty matrices) can be multiplied.
+ (where `n` is zero), to ensure that two elements of the matrix
+ space (empty matrices) can be multiplied.
+
+ Matrix algebras override this with something more useful.
"""
if self.is_trivial():
return MatrixSpace(self.base_ring(), 0)
- elif self._natural_basis is None or len(self._natural_basis) == 0:
+ elif self._matrix_basis is None or len(self._matrix_basis) == 0:
return MatrixSpace(self.base_ring(), self.dimension(), 1)
else:
- 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()
+ return self._matrix_basis[0].matrix_space()
@cached_method
sage: actual == expected
True
+ Ensure that the cached unit element (often precomputed by
+ hand) agrees with the computed one::
+
+ sage: set_random_seed()
+ sage: J = random_eja()
+ sage: cached = J.one()
+ sage: J.one.clear_cache()
+ sage: J.one() == cached
+ True
+
"""
# We can brute-force compute the matrices of the operators
# that correspond to the basis elements of this algebra.
# appeal to the "long vectors" isometry.
oper_vecs = [ _mat2vec(g.operator().matrix()) for g in self.gens() ]
- # Now we use basis linear algebra to find the coefficients,
+ # Now we use basic linear algebra to find the coefficients,
# of the matrices-as-vectors-linear-combination, which should
# work for the original algebra basis too.
- A = matrix.column(self.base_ring(), oper_vecs)
+ A = matrix(self.base_ring(), oper_vecs)
# We used the isometry on the left-hand side already, but we
# still need to do it for the right-hand side. Recall that we
# wanted something that summed to the identity matrix.
b = _mat2vec( matrix.identity(self.base_ring(), self.dimension()) )
- # Now if there's an identity element in the algebra, this should work.
- coeffs = A.solve_right(b)
- return self.linear_combination(zip(self.gens(), coeffs))
+ # Now if there's an identity element in the algebra, this
+ # should work. We solve on the left to avoid having to
+ # transpose the matrix "A".
+ return self.from_vector(A.solve_left(b))
def peirce_decomposition(self, c):
Vector space of degree 6 and dimension 2...
sage: J1
Euclidean Jordan algebra of dimension 3...
- sage: J0.one().natural_representation()
+ sage: J0.one().to_matrix()
[0 0 0]
[0 0 0]
[0 0 1]
sage: orig_df = AA.options.display_format
sage: AA.options.display_format = 'radical'
- sage: J.from_vector(J5.basis()[0]).natural_representation()
+ sage: J.from_vector(J5.basis()[0]).to_matrix()
[ 0 0 1/2*sqrt(2)]
[ 0 0 0]
[1/2*sqrt(2) 0 0]
- sage: J.from_vector(J5.basis()[1]).natural_representation()
+ sage: J.from_vector(J5.basis()[1]).to_matrix()
[ 0 0 0]
[ 0 0 1/2*sqrt(2)]
[ 0 1/2*sqrt(2) 0]
sage: AA.options.display_format = orig_df
- sage: J1.one().natural_representation()
+ sage: J1.one().to_matrix()
[1 0 0]
[0 1 0]
[0 0 0]
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
J5 = eigspace
else:
gens = tuple( self.from_vector(b) for b in eigspace.basis() )
- subalg = FiniteDimensionalEuclideanJordanSubalgebra(self, gens)
+ subalg = FiniteDimensionalEuclideanJordanSubalgebra(self,
+ gens,
+ check_axioms=False)
if eigval == 0:
J0 = subalg
elif eigval == 1:
return tuple( self.random_element(thorough)
for idx in range(count) )
- @classmethod
- def random_instance(cls, field=AA, **kwargs):
- """
- Return a random instance of this type of algebra.
-
- Beware, this will crash for "most instances" because the
- constructor below looks wrong.
- """
- if cls is TrivialEJA:
- # The TrivialEJA class doesn't take an "n" argument because
- # there's only one.
- return cls(field)
-
- n = ZZ.random_element(cls._max_test_case_size() + 1)
- return cls(n, field, **kwargs)
@cached_method
def _charpoly_coefficients(self):
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()
Ensure that computing the rank actually works, since the ranks
of all simple algebras are known and will be cached by default::
- sage: J = HadamardEJA(4)
- sage: J.rank.clear_cache()
- sage: J.rank()
- 4
-
- ::
-
- sage: J = JordanSpinEJA(4)
- sage: J.rank.clear_cache()
- sage: J.rank()
- 2
-
- ::
-
- sage: J = RealSymmetricEJA(3)
- sage: J.rank.clear_cache()
- sage: J.rank()
- 3
-
- ::
-
- sage: J = ComplexHermitianEJA(2)
- sage: J.rank.clear_cache()
- sage: J.rank()
- 2
-
- ::
+ sage: set_random_seed() # long time
+ sage: J = random_eja() # long time
+ sage: caches = J.rank() # long time
+ sage: J.rank.clear_cache() # long time
+ sage: J.rank() == cached # long time
+ True
- sage: J = QuaternionHermitianEJA(2)
- sage: J.rank.clear_cache()
- sage: J.rank()
- 2
"""
return len(self._charpoly_coefficients())
Element = FiniteDimensionalEuclideanJordanAlgebraElement
-
-class HadamardEJA(FiniteDimensionalEuclideanJordanAlgebra):
- """
- Return the Euclidean Jordan Algebra corresponding to the set
- `R^n` under the Hadamard product.
-
- Note: this is nothing more than the Cartesian product of ``n``
- copies of the spin algebra. Once Cartesian product algebras
- are implemented, this can go.
+class RationalBasisEuclideanJordanAlgebraNg(FiniteDimensionalEuclideanJordanAlgebra):
+ r"""
+ New class for algebras whose supplied basis elements have all rational entries.
SETUP::
- sage: from mjo.eja.eja_algebra import HadamardEJA
+ sage: from mjo.eja.eja_algebra import BilinearFormEJA
EXAMPLES:
- This multiplication table can be verified by hand::
+ The supplied basis is orthonormalized by default::
- sage: J = HadamardEJA(3)
- sage: e0,e1,e2 = J.gens()
- sage: e0*e0
- e0
- sage: e0*e1
- 0
- sage: e0*e2
- 0
- sage: e1*e1
- e1
- sage: e1*e2
- 0
- sage: e2*e2
- e2
+ 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]
+ )
- TESTS:
+ """
+ def __init__(self,
+ field,
+ basis,
+ jordan_product,
+ inner_product,
+ orthonormalize=True,
+ prefix='e',
+ category=None,
+ check_field=True,
+ check_axioms=True):
- We can change the generator prefix::
+ n = len(basis)
+ vector_basis = basis
- sage: HadamardEJA(3, prefix='r').gens()
- (r0, r1, r2)
+ from sage.structure.element import is_Matrix
+ basis_is_matrices = False
- """
- 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) ]
+ 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()
- fdeja = super(HadamardEJA, self)
- fdeja.__init__(field, mult_table, **kwargs)
- self.rank.set_cache(n)
+ V = VectorSpace(field, degree)
- def inner_product(self, x, y):
- """
- Faster to reimplement than to use natural representations.
+ # 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
- SETUP::
+ 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 )
- sage: from mjo.eja.eja_algebra import HadamardEJA
+ # Normalize the "matrix" basis, too!
+ basis = vector_basis
- TESTS:
+ if basis_is_matrices:
+ from mjo.eja.eja_utils import _vec2mat
+ basis = tuple( map(_vec2mat,basis) )
- Ensure that this is the usual inner product for the algebras
- over `R^n`::
+ W = V.span_of_basis( vector_basis )
- 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
+ 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) ]
- """
- return x.to_vector().inner_product(y.to_vector())
+ # 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
+
+ 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, # 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
-def random_eja(field=AA):
- """
- Return a "random" finite-dimensional Euclidean Jordan Algebra.
+ EXAMPLES:
+
+ 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)::
+
+ sage: J = JordanSpinEJA(3)
+ sage: J._charpoly_coefficients()
+ (X1^2 - X2^2 - X3^2, -2*X1)
+ sage: a0 = J._charpoly_coefficients()[0]
+ sage: J.base_ring()
+ 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)
+ 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.
+ J = FiniteDimensionalEuclideanJordanAlgebra(QQ,
+ mult_table,
+ check_field=False,
+ check_axioms=False)
+ 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.
+
+ These are the Jordan algebras whose basis, multiplication table,
+ rank, and so on are known a priori. More to the point, they are
+ the Euclidean Jordan algebras for which we are able to conjure up
+ a "random instance."
SETUP::
- sage: from mjo.eja.eja_algebra import random_eja
+ sage: from mjo.eja.eja_algebra import ConcreteEuclideanJordanAlgebra
+
+ TESTS:
+
+ Our basis is normalized with respect to the algebra's inner
+ product, unless we specify otherwise::
- TESTS::
+ sage: set_random_seed()
+ sage: J = ConcreteEuclideanJordanAlgebra.random_instance()
+ sage: all( b.norm() == 1 for b in J.gens() )
+ True
- sage: random_eja()
- Euclidean Jordan algebra of dimension...
+ Since our basis is orthonormal with respect to the algebra's inner
+ product, and since we know that this algebra is an EJA, any
+ left-multiplication operator's matrix will be symmetric because
+ natural->EJA basis representation is an isometry and within the
+ EJA the operator is self-adjoint by the Jordan axiom::
+ sage: set_random_seed()
+ sage: J = ConcreteEuclideanJordanAlgebra.random_instance()
+ sage: x = J.random_element()
+ sage: x.operator().is_self_adjoint()
+ True
"""
- classname = choice([TrivialEJA,
- HadamardEJA,
- JordanSpinEJA,
- RealSymmetricEJA,
- ComplexHermitianEJA,
- QuaternionHermitianEJA])
- return classname.random_instance(field=field)
+ @staticmethod
+ def _max_random_instance_size():
+ """
+ Return an integer "size" that is an upper bound on the size of
+ this algebra when it is used in a random test
+ case. Unfortunately, the term "size" is ambiguous -- when
+ dealing with `R^n` under either the Hadamard or Jordan spin
+ product, the "size" refers to the dimension `n`. When dealing
+ with a matrix algebra (real symmetric or complex/quaternion
+ Hermitian), it refers to the size of the matrix, which is far
+ less than the dimension of the underlying vector space.
+
+ This method must be implemented in each subclass.
+ """
+ raise NotImplementedError
+
+ @classmethod
+ def random_instance(cls, field=AA, **kwargs):
+ """
+ Return a random instance of this type of algebra.
+ This method should be implemented in each subclass.
+ """
+ from sage.misc.prandom import choice
+ eja_class = choice(cls.__subclasses__())
+ return eja_class.random_instance(field)
class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
- @staticmethod
- def _max_test_case_size():
- # Play it safe, since this will be squared and the underlying
- # field can have dimension 4 (quaternions) too.
- return 2
def __init__(self, field, basis, normalize_basis=True, **kwargs):
"""
# time to ensure that it isn't a generator expression.
basis = tuple(basis)
- if len(basis) > 1 and normalize_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.
field = field.extension(p, 'sqrt2', embedding=RLF(2).sqrt())
basis = tuple( s.change_ring(field) for s in basis )
self._basis_normalizers = tuple(
- ~(self.natural_inner_product(s,s).sqrt()) for s in basis )
+ ~(self.matrix_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 = [[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))
- fdeja = super(MatrixEuclideanJordanAlgebra, self)
- fdeja.__init__(field, Qs, natural_basis=basis, **kwargs)
- return
+ 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
Override the parent method with something that tries to compute
over a faster (non-extension) field.
"""
- if self._basis_normalizers is None:
- # We didn't normalize, so assume that the basis we started
- # with had entries in a nice 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()
- else:
- basis = ( (b/n) for (b,n) in zip(self.natural_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.
- J = MatrixEuclideanJordanAlgebra(QQ,
- basis,
- normalize_basis=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 )
-
-
- @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
+ 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 )
@staticmethod
"""
raise NotImplementedError
-
@classmethod
- def natural_inner_product(cls,X,Y):
+ def matrix_inner_product(cls,X,Y):
Xu = cls.real_unembed(X)
Yu = cls.real_unembed(Y)
tr = (Xu*Yu).trace()
- if tr in RLF:
- # It's real already.
- return tr
-
- # Otherwise, try the thing that works for complex numbers; and
- # if that doesn't work, the thing that works for quaternions.
try:
- return tr.vector()[0] # real part, imag part is index 1
+ # Works in QQ, AA, RDF, et cetera.
+ return tr.real()
except AttributeError:
- # A quaternions doesn't have a vector() method, but does
+ # A quaternion doesn't have a real() method, but does
# have coefficient_tuple() method that returns the
# coefficients of 1, i, j, and k -- in that order.
return tr.coefficient_tuple()[0]
return M
-class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra):
+class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra,
+ ConcreteEuclideanJordanAlgebra):
"""
The rank-n simple EJA consisting of real symmetric n-by-n
matrices, the usual symmetric Jordan product, and the trace inner
The dimension of this algebra is `(n^2 + n) / 2`::
sage: set_random_seed()
- sage: n_max = RealSymmetricEJA._max_test_case_size()
+ sage: n_max = RealSymmetricEJA._max_random_instance_size()
sage: n = ZZ.random_element(1, n_max)
sage: J = RealSymmetricEJA(n)
sage: J.dimension() == (n^2 + n)/2
sage: set_random_seed()
sage: J = RealSymmetricEJA.random_instance()
sage: x,y = J.random_elements(2)
- sage: actual = (x*y).natural_representation()
- sage: X = x.natural_representation()
- sage: Y = y.natural_representation()
+ sage: actual = (x*y).to_matrix()
+ sage: X = x.to_matrix()
+ sage: Y = y.to_matrix()
sage: expected = (X*Y + Y*X)/2
sage: actual == expected
True
sage: RealSymmetricEJA(3, prefix='q').gens()
(q0, q1, q2, q3, q4, q5)
- Our natural basis is normalized with respect to the natural inner
- product unless we specify otherwise::
-
- sage: set_random_seed()
- sage: J = RealSymmetricEJA.random_instance()
- sage: all( b.norm() == 1 for b in J.gens() )
- True
-
- Since our natural basis is normalized with respect to the natural
- inner product, and since we know that this algebra is an EJA, any
- left-multiplication operator's matrix will be symmetric because
- natural->EJA basis representation is an isometry and within the EJA
- the operator is self-adjoint by the Jordan axiom::
-
- sage: set_random_seed()
- sage: x = RealSymmetricEJA.random_instance().random_element()
- sage: x.operator().matrix().is_symmetric()
- True
-
We can construct the (trivial) algebra of rank zero::
sage: RealSymmetricEJA(0)
@staticmethod
- def _max_test_case_size():
+ def _max_random_instance_size():
return 4 # Dimension 10
+ @classmethod
+ def random_instance(cls, field=AA, **kwargs):
+ """
+ Return a random instance of this type of algebra.
+ """
+ n = ZZ.random_element(cls._max_random_instance_size() + 1)
+ return cls(n, field, **kwargs)
def __init__(self, n, field=AA, **kwargs):
basis = self._denormalized_basis(n, field)
- super(RealSymmetricEJA, self).__init__(field, basis, **kwargs)
+ super(RealSymmetricEJA, self).__init__(field,
+ basis,
+ check_axioms=False,
+ **kwargs)
self.rank.set_cache(n)
Embedding is a homomorphism (isomorphism, in fact)::
sage: set_random_seed()
- sage: n_max = ComplexMatrixEuclideanJordanAlgebra._max_test_case_size()
- sage: n = ZZ.random_element(n_max)
+ sage: n = ZZ.random_element(3)
sage: F = QuadraticField(-1, 'I')
sage: X = random_matrix(F, n)
sage: Y = random_matrix(F, n)
@classmethod
- def natural_inner_product(cls,X,Y):
+ def matrix_inner_product(cls,X,Y):
"""
- Compute a natural inner product in this algebra directly from
+ Compute a matrix inner product in this algebra directly from
its real embedding.
SETUP::
sage: set_random_seed()
sage: J = ComplexHermitianEJA.random_instance()
sage: x,y = J.random_elements(2)
- sage: Xe = x.natural_representation()
- sage: Ye = y.natural_representation()
+ sage: Xe = x.to_matrix()
+ sage: Ye = y.to_matrix()
sage: X = ComplexHermitianEJA.real_unembed(Xe)
sage: Y = ComplexHermitianEJA.real_unembed(Ye)
sage: expected = (X*Y).trace().real()
- sage: actual = ComplexHermitianEJA.natural_inner_product(Xe,Ye)
+ sage: actual = ComplexHermitianEJA.matrix_inner_product(Xe,Ye)
sage: actual == expected
True
"""
- return RealMatrixEuclideanJordanAlgebra.natural_inner_product(X,Y)/2
+ return RealMatrixEuclideanJordanAlgebra.matrix_inner_product(X,Y)/2
-class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra):
+class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra,
+ ConcreteEuclideanJordanAlgebra):
"""
The rank-n simple EJA consisting of complex Hermitian n-by-n
matrices over the real numbers, the usual symmetric Jordan product,
The dimension of this algebra is `n^2`::
sage: set_random_seed()
- sage: n_max = ComplexHermitianEJA._max_test_case_size()
+ sage: n_max = ComplexHermitianEJA._max_random_instance_size()
sage: n = ZZ.random_element(1, n_max)
sage: J = ComplexHermitianEJA(n)
sage: J.dimension() == n^2
sage: set_random_seed()
sage: J = ComplexHermitianEJA.random_instance()
sage: x,y = J.random_elements(2)
- sage: actual = (x*y).natural_representation()
- sage: X = x.natural_representation()
- sage: Y = y.natural_representation()
+ sage: actual = (x*y).to_matrix()
+ sage: X = x.to_matrix()
+ sage: Y = y.to_matrix()
sage: expected = (X*Y + Y*X)/2
sage: actual == expected
True
sage: ComplexHermitianEJA(2, prefix='z').gens()
(z0, z1, z2, z3)
- Our natural basis is normalized with respect to the natural inner
- product unless we specify otherwise::
-
- sage: set_random_seed()
- sage: J = ComplexHermitianEJA.random_instance()
- sage: all( b.norm() == 1 for b in J.gens() )
- True
-
- Since our natural basis is normalized with respect to the natural
- inner product, and since we know that this algebra is an EJA, any
- left-multiplication operator's matrix will be symmetric because
- natural->EJA basis representation is an isometry and within the EJA
- the operator is self-adjoint by the Jordan axiom::
-
- sage: set_random_seed()
- sage: x = ComplexHermitianEJA.random_instance().random_element()
- sage: x.operator().matrix().is_symmetric()
- True
-
We can construct the (trivial) algebra of rank zero::
sage: ComplexHermitianEJA(0)
def __init__(self, n, field=AA, **kwargs):
basis = self._denormalized_basis(n,field)
- super(ComplexHermitianEJA,self).__init__(field, basis, **kwargs)
+ super(ComplexHermitianEJA,self).__init__(field,
+ basis,
+ check_axioms=False,
+ **kwargs)
self.rank.set_cache(n)
+ @staticmethod
+ def _max_random_instance_size():
+ return 3 # Dimension 9
+
+ @classmethod
+ def random_instance(cls, field=AA, **kwargs):
+ """
+ Return a random instance of this type of algebra.
+ """
+ n = ZZ.random_element(cls._max_random_instance_size() + 1)
+ return cls(n, field, **kwargs)
class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
@staticmethod
Embedding is a homomorphism (isomorphism, in fact)::
sage: set_random_seed()
- sage: n_max = QuaternionMatrixEuclideanJordanAlgebra._max_test_case_size()
- sage: n = ZZ.random_element(n_max)
+ sage: n = ZZ.random_element(2)
sage: Q = QuaternionAlgebra(QQ,-1,-1)
sage: X = random_matrix(Q, n)
sage: Y = random_matrix(Q, n)
@classmethod
- def natural_inner_product(cls,X,Y):
+ def matrix_inner_product(cls,X,Y):
"""
- Compute a natural inner product in this algebra directly from
+ Compute a matrix inner product in this algebra directly from
its real embedding.
SETUP::
sage: set_random_seed()
sage: J = QuaternionHermitianEJA.random_instance()
sage: x,y = J.random_elements(2)
- sage: Xe = x.natural_representation()
- sage: Ye = y.natural_representation()
+ sage: Xe = x.to_matrix()
+ sage: Ye = y.to_matrix()
sage: X = QuaternionHermitianEJA.real_unembed(Xe)
sage: Y = QuaternionHermitianEJA.real_unembed(Ye)
sage: expected = (X*Y).trace().coefficient_tuple()[0]
- sage: actual = QuaternionHermitianEJA.natural_inner_product(Xe,Ye)
+ sage: actual = QuaternionHermitianEJA.matrix_inner_product(Xe,Ye)
sage: actual == expected
True
"""
- return RealMatrixEuclideanJordanAlgebra.natural_inner_product(X,Y)/4
+ return RealMatrixEuclideanJordanAlgebra.matrix_inner_product(X,Y)/4
-class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra):
- """
+class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra,
+ ConcreteEuclideanJordanAlgebra):
+ r"""
The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
matrices, the usual symmetric Jordan product, and the
real-part-of-trace inner product. It has dimension `2n^2 - n` over
The dimension of this algebra is `2*n^2 - n`::
sage: set_random_seed()
- sage: n_max = QuaternionHermitianEJA._max_test_case_size()
+ sage: n_max = QuaternionHermitianEJA._max_random_instance_size()
sage: n = ZZ.random_element(1, n_max)
sage: J = QuaternionHermitianEJA(n)
sage: J.dimension() == 2*(n^2) - n
sage: set_random_seed()
sage: J = QuaternionHermitianEJA.random_instance()
sage: x,y = J.random_elements(2)
- sage: actual = (x*y).natural_representation()
- sage: X = x.natural_representation()
- sage: Y = y.natural_representation()
+ sage: actual = (x*y).to_matrix()
+ sage: X = x.to_matrix()
+ sage: Y = y.to_matrix()
sage: expected = (X*Y + Y*X)/2
sage: actual == expected
True
sage: QuaternionHermitianEJA(2, prefix='a').gens()
(a0, a1, a2, a3, a4, a5)
- Our natural basis is normalized with respect to the natural inner
- product unless we specify otherwise::
-
- sage: set_random_seed()
- sage: J = QuaternionHermitianEJA.random_instance()
- sage: all( b.norm() == 1 for b in J.gens() )
- True
-
- Since our natural basis is normalized with respect to the natural
- inner product, and since we know that this algebra is an EJA, any
- left-multiplication operator's matrix will be symmetric because
- natural->EJA basis representation is an isometry and within the EJA
- the operator is self-adjoint by the Jordan axiom::
-
- sage: set_random_seed()
- sage: x = QuaternionHermitianEJA.random_instance().random_element()
- sage: x.operator().matrix().is_symmetric()
- True
-
We can construct the (trivial) algebra of rank zero::
sage: QuaternionHermitianEJA(0)
def __init__(self, n, field=AA, **kwargs):
basis = self._denormalized_basis(n,field)
- super(QuaternionHermitianEJA,self).__init__(field, basis, **kwargs)
+ super(QuaternionHermitianEJA,self).__init__(field,
+ basis,
+ check_axioms=False,
+ **kwargs)
+ self.rank.set_cache(n)
+
+ @staticmethod
+ def _max_random_instance_size():
+ r"""
+ The maximum rank of a random QuaternionHermitianEJA.
+ """
+ return 2 # Dimension 6
+
+ @classmethod
+ def random_instance(cls, field=AA, **kwargs):
+ """
+ Return a random instance of this type of algebra.
+ """
+ n = ZZ.random_element(cls._max_random_instance_size() + 1)
+ return cls(n, field, **kwargs)
+
+
+class HadamardEJA(RationalBasisEuclideanJordanAlgebraNg,
+ ConcreteEuclideanJordanAlgebra):
+ """
+ Return the Euclidean Jordan Algebra corresponding to the set
+ `R^n` under the Hadamard product.
+
+ Note: this is nothing more than the Cartesian product of ``n``
+ copies of the spin algebra. Once Cartesian product algebras
+ are implemented, this can go.
+
+ SETUP::
+
+ sage: from mjo.eja.eja_algebra import HadamardEJA
+
+ EXAMPLES:
+
+ This multiplication table can be verified by hand::
+
+ sage: J = HadamardEJA(3)
+ sage: e0,e1,e2 = J.gens()
+ sage: e0*e0
+ e0
+ sage: e0*e1
+ 0
+ sage: e0*e2
+ 0
+ sage: e1*e1
+ e1
+ sage: e1*e2
+ 0
+ sage: e2*e2
+ e2
+
+ TESTS:
+
+ We can change the generator prefix::
+
+ sage: HadamardEJA(3, prefix='r').gens()
+ (r0, r1, r2)
+
+ """
+ def __init__(self, n, field=AA, **kwargs):
+ V = VectorSpace(field, n)
+ basis = V.basis()
+
+ 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,
+ basis,
+ jordan_product,
+ inner_product,
+ **kwargs)
self.rank.set_cache(n)
+ if n == 0:
+ self.one.set_cache( self.zero() )
+ else:
+ self.one.set_cache( sum(self.gens()) )
+
+ @staticmethod
+ def _max_random_instance_size():
+ r"""
+ The maximum dimension of a random HadamardEJA.
+ """
+ return 5
+
+ @classmethod
+ def random_instance(cls, field=AA, **kwargs):
+ """
+ Return a random instance of this type of algebra.
+ """
+ n = ZZ.random_element(cls._max_random_instance_size() + 1)
+ return cls(n, field, **kwargs)
+
-class BilinearFormEJA(FiniteDimensionalEuclideanJordanAlgebra):
+class BilinearFormEJA(RationalBasisEuclideanJordanAlgebraNg,
+ 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 =
- (x0*y0 + <B*x_bar,y_bar>, x0*y_bar + y0*x_bar)`` where ``B`` is a
- symmetric positive-definite "bilinear form" matrix. It has
- dimension `n` over the reals, and reduces to the ``JordanSpinEJA``
- when ``B`` is the identity matrix of order ``n-1``.
+ (<Bx,y>,y_bar>, x0*y_bar + y0*x_bar)`` where `B = 1 \times B22` is
+ a symmetric positive-definite "bilinear form" matrix. Its
+ dimension is the size of `B`, and it has rank two in dimensions
+ larger than two. It reduces to the ``JordanSpinEJA`` when `B` is
+ the identity matrix of order ``n``.
+
+ We insist that the one-by-one upper-left identity block of `B` be
+ passed in as well so that we can be passed a matrix of size zero
+ to construct a trivial algebra.
SETUP::
When no bilinear form is specified, the identity matrix is used,
and the resulting algebra is the Jordan spin algebra::
- sage: J0 = BilinearFormEJA(3)
+ sage: B = matrix.identity(AA,3)
+ sage: J0 = BilinearFormEJA(B)
sage: J1 = JordanSpinEJA(3)
sage: J0.multiplication_table() == J0.multiplication_table()
True
+ An error is raised if the matrix `B` does not correspond to a
+ positive-definite bilinear form::
+
+ sage: B = matrix.random(QQ,2,3)
+ sage: J = BilinearFormEJA(B)
+ Traceback (most recent call last):
+ ...
+ ValueError: bilinear form is not positive-definite
+ sage: B = matrix.zero(QQ,3)
+ sage: J = BilinearFormEJA(B)
+ Traceback (most recent call last):
+ ...
+ ValueError: bilinear form is not positive-definite
+
TESTS:
We can create a zero-dimensional algebra::
- sage: J = BilinearFormEJA(0)
+ sage: B = matrix.identity(AA,0)
+ sage: J = BilinearFormEJA(B)
sage: J.basis()
Finite family {}
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: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
- sage: B = M.transpose()*M
- sage: J = BilinearFormEJA(n, B=B)
+ sage: B11 = matrix.identity(QQ,1)
+ sage: B22 = M.transpose()*M
+ sage: B = block_matrix(2,2,[ [B11,0 ],
+ ....: [0, B22 ] ])
+ 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)
sage: actual == expected
True
"""
- def __init__(self, n, field=AA, B=None, **kwargs):
- if B is None:
- self._B = matrix.identity(field, max(0,n-1))
- else:
- self._B = B
+ def __init__(self, B, field=AA, **kwargs):
+ 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 = x0*y0 + (self._B*xbar).inner_product(ybar)
- zbar = y0*xbar + x0*ybar
- z = V([z0] + zbar.list())
- mult_table[i][j] = z
+
+ 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,
+ V.basis(),
+ jordan_product,
+ inner_product,
+ **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).
- fdeja = super(BilinearFormEJA, self)
- fdeja.__init__(field, mult_table, **kwargs)
self.rank.set_cache(min(n,2))
- 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::
+ if n == 0:
+ self.one.set_cache( self.zero() )
+ else:
+ self.one.set_cache( self.monomial(0) )
- sage: from mjo.eja.eja_algebra import BilinearFormEJA
+ @staticmethod
+ def _max_random_instance_size():
+ r"""
+ The maximum dimension of a random BilinearFormEJA.
+ """
+ return 5
- TESTS:
+ @classmethod
+ def random_instance(cls, field=AA, **kwargs):
+ """
+ Return a random instance of this algebra.
+ """
+ n = ZZ.random_element(cls._max_random_instance_size() + 1)
+ if n.is_zero():
+ B = matrix.identity(field, n)
+ return cls(B, field, **kwargs)
- 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::
+ B11 = matrix.identity(field,1)
+ M = matrix.random(field, n-1)
+ I = matrix.identity(field, n-1)
+ alpha = field.zero()
+ while alpha.is_zero():
+ alpha = field.random_element().abs()
+ B22 = M.transpose()*M + alpha*I
- sage: set_random_seed()
- sage: n = ZZ.random_element(2,5)
- sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
- sage: B = M.transpose()*M
- sage: J = BilinearFormEJA(n, B=B)
- sage: x = J.random_element()
- sage: y = J.random_element()
- sage: x.inner_product(y) == (x*y).trace()/2
- True
+ from sage.matrix.special import block_matrix
+ B = block_matrix(2,2, [ [B11, ZZ(0) ],
+ [ZZ(0), B22 ] ])
- """
- xvec = x.to_vector()
- xbar = xvec[1:]
- yvec = y.to_vector()
- ybar = yvec[1:]
- return x[0]*y[0] + (self._B*xbar).inner_product(ybar)
+ return cls(B, field, **kwargs)
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, n, field=AA, **kwargs):
# This is a special case of the BilinearFormEJA with the identity
# matrix as its bilinear form.
- return super(JordanSpinEJA, self).__init__(n, field, **kwargs)
+ B = matrix.identity(field, n)
+ super(JordanSpinEJA, self).__init__(B, field, **kwargs)
+
+ @staticmethod
+ def _max_random_instance_size():
+ r"""
+ The maximum dimension of a random JordanSpinEJA.
+ """
+ return 5
+
+ @classmethod
+ def random_instance(cls, field=AA, **kwargs):
+ """
+ Return a random instance of this type of algebra.
+ Needed here to override the implementation for ``BilinearFormEJA``.
+ """
+ n = ZZ.random_element(cls._max_random_instance_size() + 1)
+ return cls(n, field, **kwargs)
-class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra):
+
+class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra,
+ ConcreteEuclideanJordanAlgebra):
"""
The trivial Euclidean Jordan algebra consisting of only a zero element.
"""
def __init__(self, field=AA, **kwargs):
mult_table = []
- fdeja = super(TrivialEJA, self)
+ 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
# largest subalgebra generated by any element.
- fdeja.__init__(field, mult_table, **kwargs)
self.rank.set_cache(0)
+ self.one.set_cache( self.zero() )
+ @classmethod
+ def random_instance(cls, field=AA, **kwargs):
+ # We don't take a "size" argument so the superclass method is
+ # inappropriate for us.
+ return cls(field, **kwargs)
class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra):
r"""
SETUP::
- sage: from mjo.eja.eja_algebra import (HadamardEJA,
+ sage: from mjo.eja.eja_algebra import (random_eja,
+ ....: HadamardEJA,
....: RealSymmetricEJA,
....: DirectSumEJA)
sage: J.rank()
5
+ TESTS:
+
+ The external direct sum construction is only valid when the two factors
+ have the same base ring; an error is raised otherwise::
+
+ sage: set_random_seed()
+ sage: J1 = random_eja(AA)
+ sage: J2 = random_eja(QQ)
+ sage: J = DirectSumEJA(J1,J2)
+ Traceback (most recent call last):
+ ...
+ ValueError: algebras must share the same base field
+
"""
- def __init__(self, J1, J2, field=AA, **kwargs):
+ def __init__(self, J1, J2, **kwargs):
+ if J1.base_ring() != J2.base_ring():
+ raise ValueError("algebras must share the same base field")
+ field = J1.base_ring()
+
+ self._factors = (J1, J2)
n1 = J1.dimension()
n2 = J2.dimension()
n = n1+n2
p = (J2.monomial(i)*J2.monomial(j)).to_vector()
mult_table[n1+i][n1+j] = V([field.zero()]*n1 + p.list())
- fdeja = super(DirectSumEJA, self)
- fdeja.__init__(field, mult_table, **kwargs)
+ # 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())
+
+
+ def factors(self):
+ r"""
+ Return the pair of this algebra's factors.
+
+ SETUP::
+
+ sage: from mjo.eja.eja_algebra import (HadamardEJA,
+ ....: JordanSpinEJA,
+ ....: DirectSumEJA)
+
+ EXAMPLES::
+
+ sage: J1 = HadamardEJA(2,QQ)
+ sage: J2 = JordanSpinEJA(3,QQ)
+ sage: J = DirectSumEJA(J1,J2)
+ sage: J.factors()
+ (Euclidean Jordan algebra of dimension 2 over Rational Field,
+ Euclidean Jordan algebra of dimension 3 over Rational Field)
+
+ """
+ return self._factors
+
+ def projections(self):
+ r"""
+ Return a pair of projections onto this algebra's factors.
+
+ SETUP::
+
+ sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
+ ....: ComplexHermitianEJA,
+ ....: DirectSumEJA)
+
+ EXAMPLES::
+
+ sage: J1 = JordanSpinEJA(2)
+ sage: J2 = ComplexHermitianEJA(2)
+ sage: J = DirectSumEJA(J1,J2)
+ sage: (pi_left, pi_right) = J.projections()
+ sage: J.one().to_vector()
+ (1, 0, 1, 0, 0, 1)
+ sage: pi_left(J.one()).to_vector()
+ (1, 0)
+ sage: pi_right(J.one()).to_vector()
+ (1, 0, 0, 1)
+
+ """
+ (J1,J2) = self.factors()
+ m = J1.dimension()
+ n = J2.dimension()
+ V_basis = self.vector_space().basis()
+ # 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)
+
+ def inclusions(self):
+ r"""
+ Return the pair of inclusion maps from our factors into us.
+
+ SETUP::
+
+ sage: from mjo.eja.eja_algebra import (random_eja,
+ ....: JordanSpinEJA,
+ ....: RealSymmetricEJA,
+ ....: DirectSumEJA)
+
+ EXAMPLES::
+
+ sage: J1 = JordanSpinEJA(3)
+ sage: J2 = RealSymmetricEJA(2)
+ sage: J = DirectSumEJA(J1,J2)
+ sage: (iota_left, iota_right) = J.inclusions()
+ sage: iota_left(J1.zero()) == J.zero()
+ True
+ sage: iota_right(J2.zero()) == J.zero()
+ True
+ sage: J1.one().to_vector()
+ (1, 0, 0)
+ sage: iota_left(J1.one()).to_vector()
+ (1, 0, 0, 0, 0, 0)
+ sage: J2.one().to_vector()
+ (1, 0, 1)
+ sage: iota_right(J2.one()).to_vector()
+ (0, 0, 0, 1, 0, 1)
+ sage: J.one().to_vector()
+ (1, 0, 0, 1, 0, 1)
+
+ TESTS:
+
+ Composing a projection with the corresponding inclusion should
+ produce the identity map, and mismatching them should produce
+ the zero map::
+
+ sage: set_random_seed()
+ sage: J1 = random_eja()
+ sage: J2 = random_eja()
+ sage: J = DirectSumEJA(J1,J2)
+ sage: (iota_left, iota_right) = J.inclusions()
+ sage: (pi_left, pi_right) = J.projections()
+ sage: pi_left*iota_left == J1.one().operator()
+ True
+ sage: pi_right*iota_right == J2.one().operator()
+ True
+ sage: (pi_left*iota_right).is_zero()
+ True
+ sage: (pi_right*iota_left).is_zero()
+ True
+
+ """
+ (J1,J2) = self.factors()
+ m = J1.dimension()
+ n = J2.dimension()
+ V_basis = self.vector_space().basis()
+ # 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)
+
+ def inner_product(self, x, y):
+ r"""
+ The standard Cartesian inner-product.
+
+ We project ``x`` and ``y`` onto our factors, and add up the
+ inner-products from the subalgebras.
+
+ SETUP::
+
+
+ sage: from mjo.eja.eja_algebra import (HadamardEJA,
+ ....: QuaternionHermitianEJA,
+ ....: DirectSumEJA)
+
+ EXAMPLE::
+
+ sage: J1 = HadamardEJA(3,QQ)
+ sage: J2 = QuaternionHermitianEJA(2,QQ,normalize_basis=False)
+ sage: J = DirectSumEJA(J1,J2)
+ sage: x1 = J1.one()
+ sage: x2 = x1
+ sage: y1 = J2.one()
+ sage: y2 = y1
+ sage: x1.inner_product(x2)
+ 3
+ sage: y1.inner_product(y2)
+ 2
+ sage: J.one().inner_product(J.one())
+ 5
+
+ """
+ (pi_left, pi_right) = self.projections()
+ x1 = pi_left(x)
+ x2 = pi_right(x)
+ y1 = pi_left(y)
+ y2 = pi_right(y)
+
+ return (x1.inner_product(y1) + x2.inner_product(y2))
+
+
+
+random_eja = ConcreteEuclideanJordanAlgebra.random_instance
projects / sage.d.git / blobdiff