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: scalar field is not real
+ sage: JordanSpinEJA(2, QQbar, check_field=False)
+ Euclidean Jordan algebra of dimension 2 over Algebraic Field
+
+ The multiplication table must be square with ``check_axioms=True``::
- sage: JordanSpinEJA(2,QQbar)
+ sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,((),()),((1,),))
Traceback (most recent call last):
...
- ValueError: field is not real
+ 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)
+
+ # Check commutativity of the Jordan product (symmetry of
+ # the multiplication table) and the commutativity of the
+ # inner-product (symmetry of the inner-product table)
+ # first if we're going to check them at all.. This has to
+ # be done before we define product_on_basis(), because
+ # that method assumes that self._multiplication_table is
+ # symmetric. And it has to be done before we build
+ # self._inner_product_matrix, because the process used to
+ # construct it assumes symmetry as well.
+ 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()
fda = super(FiniteDimensionalEuclideanJordanAlgebra, self)
fda.__init__(field,
- range(len(mult_table)),
+ range(n),
prefix=prefix,
category=category)
self.print_options(bracket='')
# 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
- ]
+ #
+ # 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
+
+ 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.
+ ip_matrix_constructor = lambda i,j: inner_product_table[i][j] if j <= i else inner_product_table[j][i]
+ self._inner_product_matrix = matrix(field, n, ip_matrix_constructor)
+ self._inner_product_matrix._cache = {'hermitian': True}
+ self._inner_product_matrix.set_immutable()
+
+ if check_axioms:
+ if not self._is_jordanian():
+ raise ValueError("Jordan identity does not hold")
+ if not self._inner_product_is_associative():
+ raise ValueError("inner product is not 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)
+ #
+ # We pass check=False because the matrix basis is "guaranteed"
+ # to be linearly independent... right? Ha ha.
+ V = VectorSpace(self.base_ring(), elt.nrows()*elt.ncols())
+ W = V.span_of_basis( (_mat2vec(s) for s in self.matrix_basis()),
+ check=False)
+
+ try:
+ coords = W.coordinate_vector(_mat2vec(elt))
+ except ArithmeticError: # vector is not in free module
+ raise ValueError(msg)
+ 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]
-
- def _a_regular_element(self):
- """
- Guess a regular element. Needed to compute the basis for our
- characteristic polynomial coefficients.
-
- SETUP::
-
- sage: from mjo.eja.eja_algebra import random_eja
-
- TESTS:
-
- Ensure that this hacky method succeeds for every algebra that we
- know how to construct::
+ # 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]
- sage: set_random_seed()
- sage: J = random_eja()
- sage: J._a_regular_element().is_regular()
- True
+ 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_axioms=False`` and
+ passed an invalid multiplication table.
"""
- gs = self.gens()
- z = self.sum( (i+1)*gs[i] for i in range(len(gs)) )
- if not z.is_regular():
- raise ValueError("don't know a regular element")
- return z
-
+ return all( self.product_on_basis(i,j) == self.product_on_basis(i,j)
+ for i in range(self.dimension())
+ for j in range(self.dimension()) )
- @cached_method
- def _charpoly_basis_space(self):
- """
- Return the vector space spanned by the basis used in our
- characteristic polynomial coefficients. This is used not only to
- compute those coefficients, but also any time we need to
- evaluate the coefficients (like when we compute the trace or
- determinant).
- """
- z = self._a_regular_element()
- # Don't use the parent vector space directly here in case this
- # happens to be a subalgebra. In that case, we would be e.g.
- # two-dimensional but span_of_basis() would expect three
- # coordinates.
- V = VectorSpace(self.base_ring(), self.vector_space().dimension())
- basis = [ (z**k).to_vector() for k in range(self.rank()) ]
- V1 = V.span_of_basis( basis )
- b = (V1.basis() + V1.complement().basis())
- return V.span_of_basis(b)
+ def _is_jordanian(self):
+ r"""
+ Whether or not this algebra's multiplication table respects the
+ Jordan identity `(x^{2})(xy) = x(x^{2}y)`.
+
+ We only check one arrangement of `x` and `y`, so for a
+ ``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_axioms=False`` and passed an invalid multiplication table.
+ """
+ return all( (self.monomial(i)**2)*(self.monomial(i)*self.monomial(j))
+ ==
+ (self.monomial(i))*((self.monomial(i)**2)*self.monomial(j))
+ for i in range(self.dimension())
+ for j in range(self.dimension()) )
+
+ def _inner_product_is_associative(self):
+ r"""
+ Return whether or not this algebra's inner product `B` is
+ 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_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
- @cached_method
- def _charpoly_coeff(self, i):
- """
- Return the coefficient polynomial "a_{i}" of this algebra's
- general characteristic polynomial.
-
- Having this be a separate cached method lets us compute and
- store the trace/determinant (a_{r-1} and a_{0} respectively)
- separate from the entire characteristic polynomial.
- """
- (A_of_x, x, xr, detA) = self._charpoly_matrix_system()
- R = A_of_x.base_ring()
-
- if i == self.rank():
- return R.one()
- if i > self.rank():
- # Guaranteed by theory
- return R.zero()
-
- # Danger: the in-place modification is done for performance
- # reasons (reconstructing a matrix with huge polynomial
- # entries is slow), but I don't know how cached_method works,
- # so it's highly possible that we're modifying some global
- # list variable by reference, here. In other words, you
- # probably shouldn't call this method twice on the same
- # algebra, at the same time, in two threads
- Ai_orig = A_of_x.column(i)
- A_of_x.set_column(i,xr)
- numerator = A_of_x.det()
- A_of_x.set_column(i,Ai_orig)
-
- # We're relying on the theory here to ensure that each a_i is
- # indeed back in R, and the added negative signs are to make
- # the whole charpoly expression sum to zero.
- return R(-numerator/detA)
+ 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)
+ 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
- @cached_method
- def _charpoly_matrix_system(self):
- """
- Compute the matrix whose entries A_ij are polynomials in
- X1,...,XN, the vector ``x`` of variables X1,...,XN, the vector
- corresponding to `x^r` and the determinent of the matrix A =
- [A_ij]. In other words, all of the fixed (cachable) data needed
- to compute the coefficients of the characteristic polynomial.
- """
- r = self.rank()
- n = self.dimension()
-
- # Turn my vector space into a module so that "vectors" can
- # have multivatiate polynomial entries.
- names = tuple('X' + str(i) for i in range(1,n+1))
- R = PolynomialRing(self.base_ring(), names)
-
- # Using change_ring() on the parent's vector space doesn't work
- # here because, in a subalgebra, that vector space has a basis
- # and change_ring() tries to bring the basis along with it. And
- # that doesn't work unless the new ring is a PID, which it usually
- # won't be.
- V = FreeModule(R,n)
-
- # Now let x = (X1,X2,...,Xn) be the vector whose entries are
- # indeterminates...
- x = V(names)
-
- # And figure out the "left multiplication by x" matrix in
- # that setting.
- lmbx_cols = []
- monomial_matrices = [ self.monomial(i).operator().matrix()
- for i in range(n) ] # don't recompute these!
- for k in range(n):
- ek = self.monomial(k).to_vector()
- lmbx_cols.append(
- sum( x[i]*(monomial_matrices[i]*ek)
- for i in range(n) ) )
- Lx = matrix.column(R, lmbx_cols)
-
- # Now we can compute powers of x "symbolically"
- x_powers = [self.one().to_vector(), x]
- for d in range(2, r+1):
- x_powers.append( Lx*(x_powers[-1]) )
-
- idmat = matrix.identity(R, n)
-
- W = self._charpoly_basis_space()
- W = W.change_ring(R.fraction_field())
-
- # Starting with the standard coordinates x = (X1,X2,...,Xn)
- # and then converting the entries to W-coordinates allows us
- # to pass in the standard coordinates to the charpoly and get
- # back the right answer. Specifically, with x = (X1,X2,...,Xn),
- # we have
- #
- # W.coordinates(x^2) eval'd at (standard z-coords)
- # =
- # W-coords of (z^2)
- # =
- # W-coords of (standard coords of x^2 eval'd at std-coords of z)
- #
- # We want the middle equivalent thing in our matrix, but use
- # the first equivalent thing instead so that we can pass in
- # standard coordinates.
- x_powers = [ W.coordinate_vector(xp) for xp in x_powers ]
- l2 = [idmat.column(k-1) for k in range(r+1, n+1)]
- A_of_x = matrix.column(R, n, (x_powers[:r] + l2))
- return (A_of_x, x, x_powers[r], A_of_x.det())
-
+ return True
@cached_method
- def characteristic_polynomial(self):
+ def characteristic_polynomial_of(self):
"""
- Return a characteristic polynomial that works for all elements
- of this algebra.
+ Return the algebra's "characteristic polynomial of" function,
+ which is itself a multivariate polynomial that, when evaluated
+ at the coordinates of some algebra element, returns that
+ element's characteristic polynomial.
The resulting polynomial has `n+1` variables, where `n` is the
dimension of this algebra. The first `n` variables correspond to
Alizadeh, Example 11.11::
sage: J = JordanSpinEJA(3)
- sage: p = J.characteristic_polynomial(); p
+ sage: p = J.characteristic_polynomial_of(); p
X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
sage: xvec = J.one().to_vector()
sage: p(*xvec)
any argument::
sage: J = TrivialEJA()
- sage: J.characteristic_polynomial()
+ sage: J.characteristic_polynomial_of()
1
"""
r = self.rank()
n = self.dimension()
- # The list of coefficient polynomials a_0, a_1, a_2, ..., a_n.
- a = [ self._charpoly_coeff(i) for i in range(r+1) ]
+ # The list of coefficient polynomials a_0, a_1, a_2, ..., a_(r-1).
+ a = self._charpoly_coefficients()
# We go to a bit of trouble here to reorder the
# indeterminates, so that it's easier to evaluate the
# characteristic polynomial at x's coordinates and get back
# something in terms of t, which is what we want.
- R = a[0].parent()
S = PolynomialRing(self.base_ring(),'t')
t = S.gen(0)
- S = PolynomialRing(S, R.variable_names())
- t = S(t)
+ if r > 0:
+ R = a[0].parent()
+ S = PolynomialRing(S, R.variable_names())
+ t = S(t)
+
+ 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)
- return sum( a[k]*(t**k) for k in range(len(a)) )
+ EXAMPLES::
+ sage: J = HadamardEJA(2)
+ sage: J.coordinate_polynomial_ring()
+ Multivariate Polynomial Ring in X1, X2...
+ sage: J = RealSymmetricEJA(3,QQ,orthonormalize=False)
+ 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, 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.
- 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.)
+ 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.
- Note that this will always return a matrix. The standard basis
- in `R^n` will be returned as `n`-by-`1` column matrices.
+ 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):
- """
- Return the matrix space in which this algebra's natural basis
- elements live.
+ def matrix_space(self):
"""
- if self._natural_basis is None or len(self._natural_basis) == 0:
- return MatrixSpace(self.base_ring(), self.dimension(), 1)
- else:
- return self._natural_basis[0].matrix_space()
+ 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 matrix
+ space (empty matrices) can be multiplied.
- @staticmethod
- def natural_inner_product(X,Y):
+ Matrix algebras override this with something more useful.
"""
- 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()
+ if self.is_trivial():
+ return MatrixSpace(self.base_ring(), 0)
+ elif self._matrix_basis is None or len(self._matrix_basis) == 0:
+ return MatrixSpace(self.base_ring(), self.dimension(), 1)
+ else:
+ 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().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]).to_matrix()
+ [ 0 0 1/2*sqrt(2)]
+ [ 0 0 0]
+ [1/2*sqrt(2) 0 0]
+ 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().to_matrix()
+ [1 0 0]
+ [0 1 0]
+ [0 0 0]
TESTS:
sage: J1.superalgebra() == J and J1.dimension() == J.dimension()
True
- The identity elements in the two subalgebras are the
- projections onto their respective subspaces of the
- superalgebra's identity element::
+ The decomposition is into eigenspaces, and its components are
+ therefore necessarily orthogonal. Moreover, the identity
+ elements in the two subalgebras are the projections onto their
+ respective subspaces of the superalgebra's identity element::
sage: set_random_seed()
sage: J = random_eja()
....: x = J.random_element()
sage: c = x.subalgebra_idempotent()
sage: J0,J5,J1 = J.peirce_decomposition(c)
+ sage: ipsum = 0
+ sage: for (w,y,z) in zip(J0.basis(), J5.basis(), J1.basis()):
+ ....: w = w.superalgebra_element()
+ ....: y = J.from_vector(y)
+ ....: z = z.superalgebra_element()
+ ....: ipsum += w.inner_product(y).abs()
+ ....: ipsum += w.inner_product(z).abs()
+ ....: ipsum += y.inner_product(z).abs()
+ sage: ipsum
+ 0
sage: J1(c) == J1.one()
True
sage: J0(J.one() - c) == J0.one()
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 (J0, J5, J1)
- def random_elements(self, count):
+ def random_element(self, thorough=False):
+ r"""
+ Return a random element of this algebra.
+
+ Our algebra superclass method only returns a linear
+ combination of at most two basis elements. We instead
+ want the vector space "random element" method that
+ returns a more diverse selection.
+
+ INPUT:
+
+ - ``thorough`` -- (boolean; default False) whether or not we
+ should generate irrational coefficients for the random
+ element when our base ring is irrational; this slows the
+ algebra operations to a crawl, but any truly random method
+ should include them
+
+ """
+ # For a general base ring... maybe we can trust this to do the
+ # right thing? Unlikely, but.
+ V = self.vector_space()
+ v = V.random_element()
+
+ if self.base_ring() is AA:
+ # The "random element" method of the algebraic reals is
+ # stupid at the moment, and only returns integers between
+ # -2 and 2, inclusive:
+ #
+ # https://trac.sagemath.org/ticket/30875
+ #
+ # Instead, we implement our own "random vector" method,
+ # and then coerce that into the algebra. We use the vector
+ # space degree here instead of the dimension because a
+ # subalgebra could (for example) be spanned by only two
+ # vectors, each with five coordinates. We need to
+ # generate all five coordinates.
+ if thorough:
+ v *= QQbar.random_element().real()
+ else:
+ v *= QQ.random_element()
+
+ return self.from_vector(V.coordinate_vector(v))
+
+ def random_elements(self, count, thorough=False):
"""
Return ``count`` random elements as a tuple.
+ INPUT:
+
+ - ``thorough`` -- (boolean; default False) whether or not we
+ should generate irrational coefficients for the random
+ elements when our base ring is irrational; this slows the
+ algebra operations to a crawl, but any truly random method
+ should include them
+
SETUP::
sage: from mjo.eja.eja_algebra import JordanSpinEJA
True
"""
- return tuple( self.random_element() for idx in range(count) )
+ return tuple( self.random_element(thorough)
+ for idx in range(count) )
+
@cached_method
- def rank(self):
+ def _charpoly_coefficients(self):
+ r"""
+ The `r` polynomial coefficients of the "characteristic polynomial
+ of" function.
"""
- Return the rank of this EJA.
+ n = self.dimension()
+ R = self.coordinate_polynomial_ring()
+ vars = R.gens()
+ F = R.fraction_field()
- ALGORITHM:
+ def L_x_i_j(i,j):
+ # From a result in my book, these are the entries of the
+ # basis representation of L_x.
+ return sum( vars[k]*self.monomial(k).operator().matrix()[i,j]
+ for k in range(n) )
- We first compute the polynomial "column matrices" `p_{k}` that
- evaluate to `x^k` on the coordinates of `x`. Then, we begin
- adding them to a matrix one at a time, and trying to solve the
- system that makes `p_{0}`,`p_{1}`,..., `p_{s-1}` add up to
- `p_{s}`. This will succeed only when `s` is the rank of the
- algebra, as proven in a recent draft paper of mine.
+ L_x = matrix(F, n, n, L_x_i_j)
+
+ r = None
+ if self.rank.is_in_cache():
+ r = self.rank()
+ # There's no need to pad the system with redundant
+ # columns if we *know* they'll be redundant.
+ n = r
+
+ # Compute an extra power in case the rank is equal to
+ # the dimension (otherwise, we would stop at x^(r-1)).
+ x_powers = [ (L_x**k)*self.one().to_vector()
+ for k in range(n+1) ]
+ A = matrix.column(F, x_powers[:n])
+ AE = A.extended_echelon_form()
+ E = AE[:,n:]
+ A_rref = AE[:,:n]
+ if r is None:
+ r = A_rref.rank()
+ b = x_powers[r]
+
+ # The theory says that only the first "r" coefficients are
+ # nonzero, and they actually live in the original polynomial
+ # ring and not the fraction field. We negate them because
+ # in the actual characteristic polynomial, they get moved
+ # to the other side where x^r lives.
+ return -A_rref.solve_right(E*b).change_ring(R)[:r]
+
+ @cached_method
+ def rank(self):
+ r"""
+ Return the rank of this EJA.
+
+ This is a cached method because we know the rank a priori for
+ all of the algebras we can construct. Thus we can avoid the
+ expensive ``_charpoly_coefficients()`` call unless we truly
+ need to compute the whole characteristic polynomial.
SETUP::
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: J = QuaternionHermitianEJA(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
"""
- n = self.dimension()
- if n == 0:
- return 0
- elif n == 1:
- return 1
-
- var_names = [ "X" + str(z) for z in range(1,n+1) ]
- R = PolynomialRing(self.base_ring(), var_names)
- vars = R.gens()
-
- def L_x_i_j(i,j):
- # From a result in my book, these are the entries of the
- # basis representation of L_x.
- return sum( vars[k]*self.monomial(k).operator().matrix()[i,j]
- for k in range(n) )
-
- L_x = matrix(R, n, n, L_x_i_j)
- x_powers = [ vars[k]*(L_x**k)*self.one().to_vector()
- for k in range(n) ]
-
- # Can assume n >= 2
- M = matrix([x_powers[0]])
- old_rank = 1
-
- for d in range(1,n):
- M = matrix(M.rows() + [x_powers[d]])
- M.echelonize()
- # TODO: we've basically solved the system here.
- # We should save the echelonized matrix somehow
- # so that it can be reused in the charpoly method.
- new_rank = M.rank()
- if new_rank == old_rank:
- return new_rank
- else:
- old_rank = new_rank
-
- return n
+ return len(self._charpoly_coefficients())
def vector_space(self):
Element = FiniteDimensionalEuclideanJordanAlgebraElement
-
-class KnownRankEJA(object):
- """
- A class for algebras that we actually know we can construct. The
- main issue is that, for most of our methods to make sense, we need
- to know the rank of our algebra. Thus we can't simply generate a
- "random" algebra, or even check that a given basis and product
- satisfy the axioms; because even if everything looks OK, we wouldn't
- know the rank we need to actuallty build the thing.
-
- Not really a subclass of FDEJA because doing that causes method
- resolution errors, e.g.
-
- TypeError: Error when calling the metaclass bases
- Cannot create a consistent method resolution
- order (MRO) for bases FiniteDimensionalEuclideanJordanAlgebra,
- KnownRankEJA
-
- """
- @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.
-
- 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
-
- @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)
-
-
-class HadamardEJA(FiniteDimensionalEuclideanJordanAlgebra, KnownRankEJA):
- """
- 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 RationalBasisEuclideanJordanAlgebra(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
-
- TESTS:
-
- We can change the generator prefix::
-
- sage: HadamardEJA(3, prefix='r').gens()
- (r0, r1, r2)
+ 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, 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) ]
+ def __init__(self,
+ field,
+ basis,
+ jordan_product,
+ inner_product,
+ orthonormalize=True,
+ prefix='e',
+ category=None,
+ check_field=True,
+ check_axioms=True):
- fdeja = super(HadamardEJA, self)
- fdeja.__init__(field, mult_table, **kwargs)
- self.rank.set_cache(n)
+ if check_field:
+ # Abuse the check_field parameter to check that the entries of
+ # out basis (in ambient coordinates) are in the field QQ.
+ if not all( all(b_i in QQ for b_i in b.list()) for b in basis ):
+ raise TypeError("basis not rational")
- def inner_product(self, x, y):
- """
- Faster to reimplement than to use natural representations.
+ n = len(basis)
+ vector_basis = basis
+
+ from sage.structure.element import is_Matrix
+ basis_is_matrices = False
+
+ degree = 0
+ if n > 0:
+ if is_Matrix(basis[0]):
+ basis_is_matrices = True
+ from mjo.eja.eja_utils import _vec2mat
+ vector_basis = tuple( map(_mat2vec,basis) )
+ degree = basis[0].nrows()**2
+ else:
+ degree = basis[0].degree()
+
+ V = VectorSpace(field, degree)
+
+ # 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._rational_algebra = None
+
+ if orthonormalize:
+ if self.base_ring() is not QQ:
+ # There's no point in constructing the extra algebra if this
+ # one is already rational. If the original basis is rational
+ # but normalization would make it irrational, then this whole
+ # constructor will just fail anyway as it tries to stick an
+ # irrational number into a rational algebra.
+ #
+ # Note: the same Jordan and inner-products work here,
+ # because they are necessarily defined with respect to
+ # ambient coordinates and not any particular basis.
+ self._rational_algebra = RationalBasisEuclideanJordanAlgebra(
+ QQ,
+ basis,
+ jordan_product,
+ inner_product,
+ orthonormalize=False,
+ prefix=prefix,
+ category=category,
+ check_field=False,
+ check_axioms=False)
+
+ # Compute the deorthonormalized tables before we orthonormalize
+ # the given basis. The "check" parameter here guarantees that
+ # the basis is linearly-independent.
+ W = V.span_of_basis( vector_basis, check=check_axioms)
+
+ # 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)
+
+ # We overwrite the name "vector_basis" in a second, but never modify it
+ # in place, to this effectively makes a copy of it.
+ deortho_vector_basis = vector_basis
+ self._deortho_matrix = None
+
+ if orthonormalize:
+ from mjo.eja.eja_utils import gram_schmidt
+ if basis_is_matrices:
+ vector_ip = lambda x,y: inner_product(_vec2mat(x), _vec2mat(y))
+ vector_basis = gram_schmidt(vector_basis, vector_ip)
+ else:
+ vector_basis = gram_schmidt(vector_basis, inner_product)
+
+ # Normalize the "matrix" basis, too!
+ basis = vector_basis
+
+ if basis_is_matrices:
+ basis = tuple( map(_vec2mat,basis) )
+
+ W = V.span_of_basis( vector_basis, check=check_axioms)
+
+ # Now "W" is the vector space of our algebra coordinates. The
+ # variables "X1", "X2",... refer to the entries of vectors in
+ # W. Thus to convert back and forth between the orthonormal
+ # coordinates and the given ones, we need to stick the original
+ # basis in W.
+ U = V.span_of_basis( deortho_vector_basis, check=check_axioms)
+ self._deortho_matrix = matrix( U.coordinate_vector(q)
+ for q in vector_basis )
+
+ # If the superclass constructor is going to verify the
+ # symmetry of this table, it has better at least be
+ # square...
+ if check_axioms:
+ 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) ]
+ else:
+ mult_table = [ [0 for j in range(i+1)] for i in range(n) ]
+ ip_table = [ [0 for j in range(i+1)] for i in range(n) ]
- SETUP::
+ # 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
+ ip_table[i][j] = ip
+ if check_axioms:
+ # The tables are square if we're verifying that they
+ # are commutative.
+ mult_table[j][i] = elt
+ 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 )
- sage: from mjo.eja.eja_algebra import HadamardEJA
+ super().__init__(field,
+ mult_table,
+ ip_table,
+ prefix,
+ category,
+ basis, # matrix basis
+ check_field,
+ check_axioms)
- TESTS:
+ @cached_method
+ def _charpoly_coefficients(self):
+ r"""
+ SETUP::
- Ensure that this is the usual inner product for the algebras
- over `R^n`::
+ sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
+ ....: JordanSpinEJA)
- 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
+ EXAMPLES:
- """
- return x.to_vector().inner_product(y.to_vector())
+ 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 or self._rational_algebra is None:
+ # There's no need to construct *another* algebra over the
+ # rationals if this one is already over the
+ # rationals. Likewise, if we never orthonormalized our
+ # basis, we might as well just use the given one.
+ superclass = super(RationalBasisEuclideanJordanAlgebra, self)
+ return superclass._charpoly_coefficients()
+
+ # Do the computation over the rationals. The answer will be
+ # the same, because all we've done is a change of basis.
+ # Then, change back from QQ to our real base ring
+ a = ( a_i.change_ring(self.base_ring())
+ for a_i in self._rational_algebra._charpoly_coefficients() )
+
+ # Now convert the coordinate variables back to the
+ # deorthonormalized ones.
+ R = self.coordinate_polynomial_ring()
+ from sage.modules.free_module_element import vector
+ X = vector(R, R.gens())
+ BX = self._deortho_matrix*X
+
+ subs_dict = { X[i]: BX[i] for i in range(len(X)) }
+ return tuple( a_i.subs(subs_dict) for a_i in a )
+
+class ConcreteEuclideanJordanAlgebra(RationalBasisEuclideanJordanAlgebra):
+ r"""
+ A class for the Euclidean Jordan algebras that we know by name.
-def random_eja(field=AA, nontrivial=False):
- """
- Return a "random" finite-dimensional Euclidean Jordan Algebra.
+ 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
-
- TESTS::
-
- sage: random_eja()
- Euclidean Jordan algebra of dimension...
-
- """
- eja_classes = [HadamardEJA,
- JordanSpinEJA,
- RealSymmetricEJA,
- ComplexHermitianEJA,
- QuaternionHermitianEJA]
- if not nontrivial:
- eja_classes.append(TrivialEJA)
- classname = choice(eja_classes)
- return classname.random_instance(field=field)
+ 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::
+ sage: set_random_seed()
+ sage: J = ConcreteEuclideanJordanAlgebra.random_instance()
+ sage: all( b.norm() == 1 for b in J.gens() )
+ True
+ 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
+ """
-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):
+ def _max_random_instance_size():
"""
- Compared to the superclass constructor, we take a basis instead of
- a multiplication table because the latter can be computed in terms
- of the former when the product is known (like it is here).
- """
- # Used in this class's fast _charpoly_coeff() override.
- self._basis_normalizers = None
-
- # We're going to loop through this a few times, so now's a good
- # time to ensure that it isn't a generator expression.
- basis = tuple(basis)
-
- if len(basis) > 1 and normalize_basis:
- # We'll need sqrt(2) to normalize the basis, and this
- # winds up in the multiplication table, so the whole
- # algebra needs to be over the field extension.
- R = PolynomialRing(field, 'z')
- z = R.gen()
- p = z**2 - 2
- if p.is_irreducible():
- field = field.extension(p, 'sqrt2', embedding=RLF(2).sqrt())
- basis = tuple( s.change_ring(field) for s in basis )
- self._basis_normalizers = tuple(
- ~(self.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)
-
- fdeja = super(MatrixEuclideanJordanAlgebra, self)
- fdeja.__init__(field, Qs, natural_basis=basis, **kwargs)
- return
-
+ 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.
- @cached_method
- def rank(self):
- r"""
- Override the parent method with something that tries to compute
- over a faster (non-extension) field.
+ This method must be implemented in each subclass.
"""
- if self._basis_normalizers is None:
- # We didn't normalize, so assume that the basis we started
- # with had entries in a nice field.
- return super(MatrixEuclideanJordanAlgebra, self).rank()
- 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)
- return J.rank()
+ raise NotImplementedError
- @cached_method
- def _charpoly_coeff(self, i):
- """
- Override the parent method with something that tries to compute
- over a faster (non-extension) field.
+ @classmethod
+ def random_instance(cls, *args, **kwargs):
"""
- if self._basis_normalizers is None:
- # We didn't normalize, so assume that the basis we started
- # with had entries in a nice field.
- return super(MatrixEuclideanJordanAlgebra, self)._charpoly_coeff(i)
- 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.
- J = MatrixEuclideanJordanAlgebra(QQ,
- basis,
- normalize_basis=False)
- (_,x,_,_) = J._charpoly_matrix_system()
- p = J._charpoly_coeff(i)
- # p might be missing some vars, have to substitute "optionally"
- pairs = zip(x.base_ring().gens(), self._basis_normalizers)
- substitutions = { v: v*c for (v,c) in pairs }
- result = p.subs(substitutions)
-
- # The result of "subs" can be either a coefficient-ring
- # element or a polynomial. Gotta handle both cases.
- if result in QQ:
- return self.base_ring()(result)
- else:
- return result.change_ring(self.base_ring())
-
+ Return a random instance of this type of algebra.
- @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.
- 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))
+ This method should be implemented in each subclass.
+ """
+ from sage.misc.prandom import choice
+ eja_class = choice(cls.__subclasses__())
- return mult_table
+ # These all bubble up to the RationalBasisEuclideanJordanAlgebra
+ # superclass constructor, so any (kw)args valid there are also
+ # valid here.
+ return eja_class.random_instance(*args, **kwargs)
+class MatrixEuclideanJordanAlgebra:
@staticmethod
def real_embed(M):
"""
"""
raise NotImplementedError
+ @staticmethod
+ def jordan_product(X,Y):
+ return (X*Y + Y*X)/2
@classmethod
- def natural_inner_product(cls,X,Y):
+ def trace_inner_product(cls,X,Y):
Xu = cls.real_unembed(X)
Yu = cls.real_unembed(Y)
tr = (Xu*Yu).trace()
- 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, KnownRankEJA):
+class RealSymmetricEJA(ConcreteEuclideanJordanAlgebra,
+ RealMatrixEuclideanJordanAlgebra):
"""
The rank-n simple EJA consisting of real symmetric n-by-n
matrices, the usual symmetric Jordan product, and the trace inner
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
+ We can construct the (trivial) algebra of rank zero::
- 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
+ sage: RealSymmetricEJA(0)
+ Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
"""
@classmethod
else:
Sij = Eij + Eij.transpose()
S.append(Sij)
- return S
+ return tuple(S)
@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,
+ self.jordan_product,
+ self.trace_inner_product,
+ **kwargs)
self.rank.set_cache(n)
+ self.one.set_cache(self(matrix.identity(field,n)))
class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
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 trace_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.trace_inner_product(Xe,Ye)
sage: actual == expected
True
"""
- return RealMatrixEuclideanJordanAlgebra.natural_inner_product(X,Y)/2
+ return RealMatrixEuclideanJordanAlgebra.trace_inner_product(X,Y)/2
-class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra, KnownRankEJA):
+class ComplexHermitianEJA(ConcreteEuclideanJordanAlgebra,
+ ComplexMatrixEuclideanJordanAlgebra):
"""
The rank-n simple EJA consisting of complex Hermitian n-by-n
matrices over the real numbers, the usual symmetric Jordan product,
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
+ We can construct the (trivial) algebra of rank zero::
- 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
+ sage: ComplexHermitianEJA(0)
+ Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
"""
# Since we embedded these, we can drop back to the "field" that we
# started with instead of the complex extension "F".
- return ( s.change_ring(field) for s in S )
+ return tuple( s.change_ring(field) for s in S )
def __init__(self, n, field=AA, **kwargs):
basis = self._denormalized_basis(n,field)
- super(ComplexHermitianEJA,self).__init__(field, basis, **kwargs)
+ super(ComplexHermitianEJA, self).__init__(field,
+ basis,
+ self.jordan_product,
+ self.trace_inner_product,
+ **kwargs)
self.rank.set_cache(n)
+ # TODO: pre-cache the identity!
+
+ @staticmethod
+ def _max_random_instance_size():
+ 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 trace_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.trace_inner_product(Xe,Ye)
sage: actual == expected
True
"""
- return RealMatrixEuclideanJordanAlgebra.natural_inner_product(X,Y)/4
+ return RealMatrixEuclideanJordanAlgebra.trace_inner_product(X,Y)/4
-class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra,
- KnownRankEJA):
- """
+class QuaternionHermitianEJA(ConcreteEuclideanJordanAlgebra,
+ QuaternionMatrixEuclideanJordanAlgebra):
+ r"""
The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
matrices, the usual symmetric Jordan product, and the
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::
+ We can construct the (trivial) algebra of rank zero::
- sage: set_random_seed()
- sage: x = QuaternionHermitianEJA.random_instance().random_element()
- sage: x.operator().matrix().is_symmetric()
- True
+ sage: QuaternionHermitianEJA(0)
+ Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
"""
@classmethod
# Since we embedded these, we can drop back to the "field" that we
# started with instead of the quaternion algebra "Q".
- return ( s.change_ring(field) for s in S )
+ return tuple( s.change_ring(field) for s in S )
def __init__(self, n, field=AA, **kwargs):
basis = self._denormalized_basis(n,field)
- super(QuaternionHermitianEJA,self).__init__(field, basis, **kwargs)
+ super(QuaternionHermitianEJA, self).__init__(field,
+ basis,
+ self.jordan_product,
+ self.trace_inner_product,
+ **kwargs)
+ self.rank.set_cache(n)
+ # TODO: cache one()!
+
+ @staticmethod
+ def _max_random_instance_size():
+ 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(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)
+
+ # Don't orthonormalize because our basis is already
+ # orthonormal with respect to our inner-product. But also
+ # don't pass check_field=False here, because the user can pass
+ # in a field!
+ super(HadamardEJA, self).__init__(field,
+ basis,
+ jordan_product,
+ inner_product,
+ orthonormalize=False,
+ check_axioms=False,
+ **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, KnownRankEJA):
+class BilinearFormEJA(ConcreteEuclideanJordanAlgebra):
r"""
The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
with the half-trace inner product and jordan product ``x*y =
- (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)
+ # This is a special case of the BilinearFormEJA with the
+ # identity matrix as its bilinear form.
+ B = matrix.identity(field, n)
+
+ # Don't orthonormalize because our basis is already
+ # orthonormal with respect to our inner-product. But
+ # also don't pass check_field=False here, because the
+ # user can pass in a field!
+ super(JordanSpinEJA, self).__init__(B,
+ field,
+ orthonormalize=False,
+ check_axioms=False,
+ **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, KnownRankEJA):
+class TrivialEJA(ConcreteEuclideanJordanAlgebra):
"""
The trivial Euclidean Jordan algebra consisting of only a zero element.
"""
def __init__(self, field=AA, **kwargs):
- mult_table = []
- fdeja = super(TrivialEJA, self)
+ jordan_product = lambda x,y: x
+ inner_product = lambda x,y: field(0)
+ basis = ()
+ super(TrivialEJA, self).__init__(field,
+ basis,
+ jordan_product,
+ inner_product,
+ **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"""
+ The external (orthogonal) direct sum of two other Euclidean Jordan
+ algebras. Essentially the Cartesian product of its two factors.
+ Every Euclidean Jordan algebra decomposes into an orthogonal
+ direct sum of simple Euclidean Jordan algebras, so no generality
+ is lost by providing only this construction.
+
+ SETUP::
+
+ sage: from mjo.eja.eja_algebra import (random_eja,
+ ....: HadamardEJA,
+ ....: RealSymmetricEJA,
+ ....: DirectSumEJA)
+
+ EXAMPLES::
+
+ sage: J1 = HadamardEJA(2)
+ sage: J2 = RealSymmetricEJA(3)
+ sage: J = DirectSumEJA(J1,J2)
+ sage: J.dimension()
+ 8
+ 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,orthonormalize=False)
+ sage: J = DirectSumEJA(J1,J2)
+ Traceback (most recent call last):
+ ...
+ ValueError: algebras must share the same base field
+
+ """
+ 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
+ V = VectorSpace(field, n)
+ mult_table = [ [ V.zero() for j in range(i+1) ]
+ for i in range(n) ]
+ for i in range(n1):
+ for j in range(i+1):
+ p = (J1.monomial(i)*J1.monomial(j)).to_vector()
+ mult_table[i][j] = V(p.list() + [field.zero()]*n2)
+
+ for i in range(n2):
+ for j in range(i+1):
+ 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 = [ [ field.zero() for j in range(i+1) ]
+ for i in range(n) ]
+ 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,orthonormalize=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