coords = W.coordinate_vector(_mat2vec(elt))
return self.from_vector(coords)
+ @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
def _repr_(self):
"""
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::
-
- sage: set_random_seed()
- sage: J = random_eja()
- sage: J._a_regular_element().is_regular()
- True
-
- """
- 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
-
-
- @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)
-
-
-
- @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)
-
-
- @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())
-
-
@cached_method
def characteristic_polynomial(self):
"""
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 sum( a[k]*(t**k) for k in range(len(a)) )
+ return (t**r + sum( a[k]*(t**k) for k in range(r) ))
def inner_product(self, x, y):
"""
Return the matrix space in which this algebra's natural basis
elements live.
+
+ 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.
"""
- if self._natural_basis is None or len(self._natural_basis) == 0:
+ if self.is_trivial():
+ return MatrixSpace(self.base_ring(), 0)
+ elif 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 tuple( self.random_element() 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 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()
+ var_names = [ "X" + str(z) for z in range(1,n+1) ]
+ R = PolynomialRing(self.base_ring(), var_names)
+ 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) )
+
+ 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.
- 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.
+ 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::
sage: J.rank.clear_cache()
sage: J.rank()
2
-
"""
- 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):
+class HadamardEJA(FiniteDimensionalEuclideanJordanAlgebra):
"""
Return the Euclidean Jordan Algebra corresponding to the set
`R^n` under the Hadamard product.
return x.to_vector().inner_product(y.to_vector())
-def random_eja(field=AA, nontrivial=False):
+def random_eja(field=AA):
"""
Return a "random" finite-dimensional Euclidean Jordan Algebra.
Euclidean Jordan algebra of dimension...
"""
- eja_classes = KnownRankEJA.__subclasses__()
- if nontrivial:
- eja_classes.remove(TrivialEJA)
- classname = choice(eja_classes)
+ classname = choice([TrivialEJA,
+ HadamardEJA,
+ JordanSpinEJA,
+ RealSymmetricEJA,
+ ComplexHermitianEJA,
+ QuaternionHermitianEJA])
return classname.random_instance(field=field)
-
-
class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
@staticmethod
def _max_test_case_size():
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.
+ # Used in this class's fast _charpoly_coefficients() override.
self._basis_normalizers = None
# We're going to loop through this a few times, so now's a good
@cached_method
- def rank(self):
+ def _charpoly_coefficients(self):
r"""
Override the parent method with something that tries to compute
over a faster (non-extension) field.
if self._basis_normalizers is None:
# We didn't normalize, so assume that the basis we started
# with had entries in a nice field.
- return super(MatrixEuclideanJordanAlgebra, self).rank()
+ return super(MatrixEuclideanJordanAlgebra, self)._charpoly_coefficients()
else:
basis = ( (b/n) for (b,n) in zip(self.natural_basis(),
self._basis_normalizers) )
J = MatrixEuclideanJordanAlgebra(QQ,
basis,
normalize_basis=False)
- return J.rank()
-
- @cached_method
- def _charpoly_coeff(self, i):
- """
- 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.
- 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())
+ 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
# 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()
return M
-class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra, KnownRankEJA):
+class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra):
"""
The rank-n simple EJA consisting of real symmetric n-by-n
matrices, the usual symmetric Jordan product, and the trace inner
sage: x.operator().matrix().is_symmetric()
True
+ We can construct the (trivial) algebra of rank zero::
+
+ sage: RealSymmetricEJA(0)
+ Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
+
"""
@classmethod
def _denormalized_basis(cls, n, field):
return RealMatrixEuclideanJordanAlgebra.natural_inner_product(X,Y)/2
-class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra, KnownRankEJA):
+class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra):
"""
The rank-n simple EJA consisting of complex Hermitian n-by-n
matrices over the real numbers, the usual symmetric Jordan product,
sage: x.operator().matrix().is_symmetric()
True
+ We can construct the (trivial) algebra of rank zero::
+
+ sage: ComplexHermitianEJA(0)
+ Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
+
"""
@classmethod
return RealMatrixEuclideanJordanAlgebra.natural_inner_product(X,Y)/4
-class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra,
- KnownRankEJA):
+class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra):
"""
The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
matrices, the usual symmetric Jordan product, and the
sage: x.operator().matrix().is_symmetric()
True
+ We can construct the (trivial) algebra of rank zero::
+
+ sage: QuaternionHermitianEJA(0)
+ Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
+
"""
@classmethod
def _denormalized_basis(cls, n, field):
self.rank.set_cache(n)
-class BilinearFormEJA(FiniteDimensionalEuclideanJordanAlgebra, KnownRankEJA):
+class BilinearFormEJA(FiniteDimensionalEuclideanJordanAlgebra):
r"""
The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
with the half-trace inner product and jordan product ``x*y =
return super(JordanSpinEJA, self).__init__(n, field, **kwargs)
-class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra, KnownRankEJA):
+class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra):
"""
The trivial Euclidean Jordan algebra consisting of only a zero element.
projects / sage.d.git / blobdiff