X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=mjo%2Feja%2Feja_algebra.py;h=d4ee9b022d8524752d3e15edf37284c2381ce509;hb=ec7dbfb6ce0054f55280412e43870b4019abc40c;hp=362b03a6695a00061034cb9168ba76021a9e2daf;hpb=a7cf81951d0cb51ea40d9362a75385204596df42;p=sage.d.git diff --git a/mjo/eja/eja_algebra.py b/mjo/eja/eja_algebra.py index 362b03a..d4ee9b0 100644 --- a/mjo/eja/eja_algebra.py +++ b/mjo/eja/eja_algebra.py @@ -5,9 +5,9 @@ are used in optimization, and have some additional nice methods beyond what can be supported in a general Jordan Algebra. """ -from sage.algebras.finite_dimensional_algebras.finite_dimensional_algebra import FiniteDimensionalAlgebra from sage.algebras.quatalg.quaternion_algebra import QuaternionAlgebra from sage.categories.finite_dimensional_algebras_with_basis import FiniteDimensionalAlgebrasWithBasis +from sage.combinat.free_module import CombinatorialFreeModule from sage.matrix.constructor import matrix from sage.misc.cachefunc import cached_method from sage.misc.prandom import choice @@ -17,53 +17,16 @@ from sage.rings.number_field.number_field import QuadraticField from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing from sage.rings.rational_field import QQ from sage.structure.element import is_Matrix -from sage.structure.category_object import normalize_names from mjo.eja.eja_element import FiniteDimensionalEuclideanJordanAlgebraElement -from mjo.eja.eja_utils import _vec2mat, _mat2vec - -class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra): - @staticmethod - def __classcall_private__(cls, - field, - mult_table, - rank, - names='e', - assume_associative=False, - category=None, - natural_basis=None): - n = len(mult_table) - mult_table = [b.base_extend(field) for b in mult_table] - for b in mult_table: - b.set_immutable() - if not (is_Matrix(b) and b.dimensions() == (n, n)): - raise ValueError("input is not a multiplication table") - mult_table = tuple(mult_table) - - cat = FiniteDimensionalAlgebrasWithBasis(field) - cat.or_subcategory(category) - if assume_associative: - cat = cat.Associative() - - names = normalize_names(n, names) - - fda = super(FiniteDimensionalEuclideanJordanAlgebra, cls) - return fda.__classcall__(cls, - field, - mult_table, - rank, - assume_associative=assume_associative, - names=names, - category=cat, - natural_basis=natural_basis) - +from mjo.eja.eja_utils import _mat2vec +class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): def __init__(self, field, mult_table, rank, - names='e', - assume_associative=False, + prefix='e', category=None, natural_basis=None): """ @@ -85,12 +48,88 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra): """ self._rank = rank self._natural_basis = natural_basis - self._multiplication_table = mult_table + + if category is None: + category = FiniteDimensionalAlgebrasWithBasis(field).Unital() fda = super(FiniteDimensionalEuclideanJordanAlgebra, self) fda.__init__(field, - mult_table, - names=names, + range(len(mult_table)), + prefix=prefix, category=category) + self.print_options(bracket='') + + # The multiplication table we're given is necessarily in terms + # of vectors, because we don't have an algebra yet for + # anything to be an element of. However, it's faster in the + # 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 = matrix( + [ map(lambda x: self.from_vector(x), ls) + for ls in mult_table ] ) + self._multiplication_table.set_immutable() + + + def _element_constructor_(self, elt): + """ + Construct an element of this algebra from its natural + representation. + + This gets called only after the parent element _call_ method + fails to find a coercion for the argument. + + SETUP:: + + sage: from mjo.eja.eja_algebra import (JordanSpinEJA, + ....: RealCartesianProductEJA, + ....: RealSymmetricEJA) + + EXAMPLES: + + The identity in `S^n` is converted to the identity in the EJA:: + + sage: J = RealSymmetricEJA(3) + sage: I = matrix.identity(QQ,3) + sage: J(I) == J.one() + True + + This skew-symmetric matrix can't be represented in the EJA:: + + sage: J = RealSymmetricEJA(3) + sage: A = matrix(QQ,3, lambda i,j: i-j) + sage: J(A) + Traceback (most recent call last): + ... + ArithmeticError: vector is not in free module + + 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:: + + sage: set_random_seed() + sage: J = RealCartesianProductEJA(5) + sage: x = J.random_element() + sage: J(x.to_vector().column()) == x + True + sage: J = JordanSpinEJA(5) + sage: x = J.random_element() + sage: J(x.to_vector().column()) == x + True + + """ + natural_basis = self.natural_basis() + if elt not in natural_basis[0].matrix_space(): + raise ValueError("not a naturally-represented algebra element") + + # Thanks for nothing! Matrix spaces aren't vector + # spaces in Sage, so we have to figure out its + # natural-basis coordinates ourselves. + V = VectorSpace(elt.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) def _repr_(self): @@ -111,9 +150,12 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra): Euclidean Jordan algebra of degree 3 over Real Double Field """ + # TODO: change this to say "dimension" and fix all the tests. fmt = "Euclidean Jordan algebra of degree {} over {}" - return fmt.format(self.degree(), self.base_ring()) + 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): """ @@ -153,7 +195,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra): """ z = self._a_regular_element() V = self.vector_space() - V1 = V.span_of_basis( (z**k).vector() for k in range(self.rank()) ) + V1 = V.span_of_basis( (z**k).to_vector() for k in range(self.rank()) ) b = (V1.basis() + V1.complement().basis()) return V.span_of_basis(b) @@ -205,11 +247,14 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra): n = self.dimension() # Construct a new algebra over a multivariate polynomial ring... - names = ['X' + str(i) for i in range(1,n+1)] + names = tuple('X' + str(i) for i in range(1,n+1)) R = PolynomialRing(self.base_ring(), names) - J = FiniteDimensionalEuclideanJordanAlgebra(R, - self._multiplication_table, - r) + # Hack around the fact that our multiplication table is in terms of + # algebra elements but the constructor wants it in terms of vectors. + vmt = [ tuple([ self._multiplication_table[i,j].to_vector() + for j in range(self._multiplication_table.nrows()) ]) + for i in range(self._multiplication_table.ncols()) ] + J = FiniteDimensionalEuclideanJordanAlgebra(R, tuple(vmt), r) idmat = matrix.identity(J.base_ring(), n) @@ -231,17 +276,17 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra): # 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 = J(W(R.gens())) + x = J.from_vector(W(R.gens())) # Handle the zeroth power separately, because computing # the unit element in J is mathematically suspect. - x0 = W.coordinates(self.one().vector()) - l1 = [ matrix.column(x0) ] - l1 += [ matrix.column(W.coordinates((x**k).vector())) + x0 = W.coordinate_vector(self.one().to_vector()) + l1 = [ x0.column() ] + l1 += [ W.coordinate_vector((x**k).to_vector()).column() for k in range(1,r) ] l2 = [idmat.column(k-1).column() for k in range(r+1, n+1)] A_of_x = matrix.block(R, 1, n, (l1 + l2)) - xr = W.coordinates((x**r).vector()) + xr = W.coordinate_vector((x**r).to_vector()) return (A_of_x, x, xr, A_of_x.det()) @@ -271,7 +316,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra): sage: J = JordanSpinEJA(3) sage: p = J.characteristic_polynomial(); p X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2 - sage: xvec = J.one().vector() + sage: xvec = J.one().to_vector() sage: p(*xvec) t^2 - 2*t + 1 @@ -337,6 +382,41 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra): return x.trace_inner_product(y) + def multiplication_table(self): + """ + Return a readable matrix representation of this algebra's + multiplication table. The (i,j)th entry in the matrix contains + the product of the ith basis element with the jth. + + This is not extraordinarily useful, but it overrides a superclass + method that would otherwise just crash and complain about the + algebra being infinite. + + SETUP:: + + sage: from mjo.eja.eja_algebra import (JordanSpinEJA, + ....: RealCartesianProductEJA) + + EXAMPLES:: + + sage: J = RealCartesianProductEJA(3) + sage: J.multiplication_table() + [e0 0 0] + [ 0 e1 0] + [ 0 0 e2] + + :: + + sage: J = JordanSpinEJA(3) + sage: J.multiplication_table() + [e0 e1 e2] + [e1 e0 0] + [e2 0 e0] + + """ + return self._multiplication_table + + def natural_basis(self): """ Return a more-natural representation of this algebra's basis. @@ -359,7 +439,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra): sage: J = RealSymmetricEJA(2) sage: J.basis() - Family (e0, e1, e2) + Finite family {0: e0, 1: e1, 2: e2} sage: J.natural_basis() ( [1 0] [0 1] [0 0] @@ -370,7 +450,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra): sage: J = JordanSpinEJA(2) sage: J.basis() - Family (e0, e1) + Finite family {0: e0, 1: e1} sage: J.natural_basis() ( [1] [0] @@ -379,7 +459,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra): """ if self._natural_basis is None: - return tuple( b.vector().column() for b in self.basis() ) + return tuple( b.to_vector().column() for b in self.basis() ) else: return self._natural_basis @@ -400,7 +480,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra): sage: J.one() e0 + e1 + e2 + e3 + e4 - TESTS:: + TESTS: The identity element acts like the identity:: @@ -442,7 +522,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra): # Now if there's an identity element in the algebra, this should work. coeffs = A.solve_right(b) - return self.linear_combination(zip(coeffs,self.gens())) + return self.linear_combination(zip(self.gens(), coeffs)) def rank(self): @@ -520,7 +600,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra): Vector space of dimension 3 over Rational Field """ - return self.zero().vector().parent().ambient_vector_space() + return self.zero().to_vector().parent().ambient_vector_space() Element = FiniteDimensionalEuclideanJordanAlgebraElement @@ -559,18 +639,13 @@ class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra): e2 """ - @staticmethod - def __classcall_private__(cls, n, field=QQ): - # The FiniteDimensionalAlgebra constructor takes a list of - # matrices, the ith representing right multiplication by the ith - # basis element in the vector space. So if e_1 = (1,0,0), then - # right (Hadamard) multiplication of x by e_1 picks out the first - # component of x; and likewise for the ith basis element e_i. - Qs = [ matrix(field, n, n, lambda k,j: 1*(k == j == i)) - for i in xrange(n) ] - - fdeja = super(RealCartesianProductEJA, cls) - return fdeja.__classcall_private__(cls, field, Qs, rank=n) + def __init__(self, n, field=QQ): + V = VectorSpace(field, n) + mult_table = [ [ V.gen(i)*(i == j) for j in range(n) ] + for i in range(n) ] + + fdeja = super(RealCartesianProductEJA, self) + return fdeja.__init__(field, mult_table, rank=n) def inner_product(self, x, y): return _usual_ip(x,y) @@ -741,10 +816,7 @@ def _multiplication_table_from_matrix_basis(basis): 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. A reordered copy - of the basis is also returned to work around the fact that - the ``span()`` in this function will change the order of the basis - from what we think it is, to... something else. + 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 @@ -755,30 +827,15 @@ def _multiplication_table_from_matrix_basis(basis): dimension = basis[0].nrows() V = VectorSpace(field, dimension**2) - W = V.span( _mat2vec(s) for s in basis ) - - # Taking the span above reorders our basis (thanks, jerk!) so we - # need to put our "matrix basis" in the same order as the - # (reordered) vector basis. - S = tuple( _vec2mat(b) for b in W.basis() ) - - Qs = [] - for s in S: - # Brute force the multiplication-by-s matrix by looping - # through all elements of the basis and doing the computation - # to find out what the corresponding row should be. BEWARE: - # these multiplication tables won't be symmetric! It therefore - # becomes REALLY IMPORTANT that the underlying algebra - # constructor uses ROW vectors and not COLUMN vectors. That's - # why we're computing rows here and not columns. - Q_rows = [] - for t in S: - this_row = _mat2vec((s*t + t*s)/2) - Q_rows.append(W.coordinates(this_row)) - Q = matrix(field, W.dimension(), Q_rows) - Qs.append(Q) - - return (Qs, S) + W = V.span_of_basis( _mat2vec(s) for s in basis ) + n = len(basis) + mult_table = [[W.zero() for j in range(n)] for i in range(n)] + for i in range(n): + for j in range(n): + mat_entry = (basis[i]*basis[j] + basis[j]*basis[i])/2 + mult_table[i][j] = W.coordinate_vector(_mat2vec(mat_entry)) + + return mult_table def _embed_complex_matrix(M): @@ -1009,7 +1066,7 @@ def _unembed_quaternion_matrix(M): # The usual inner product on R^n. def _usual_ip(x,y): - return x.vector().inner_product(y.vector()) + return x.to_vector().inner_product(y.to_vector()) # The inner product used for the real symmetric simple EJA. # We keep it as a separate function because e.g. the complex @@ -1043,12 +1100,12 @@ class RealSymmetricEJA(FiniteDimensionalEuclideanJordanAlgebra): TESTS: - The degree of this algebra is `(n^2 + n) / 2`:: + The dimension of this algebra is `(n^2 + n) / 2`:: sage: set_random_seed() sage: n = ZZ.random_element(1,5) sage: J = RealSymmetricEJA(n) - sage: J.degree() == (n^2 + n)/2 + sage: J.dimension() == (n^2 + n)/2 True The Jordan multiplication is what we think it is:: @@ -1068,17 +1125,15 @@ class RealSymmetricEJA(FiniteDimensionalEuclideanJordanAlgebra): True """ - @staticmethod - def __classcall_private__(cls, n, field=QQ): + def __init__(self, n, field=QQ): S = _real_symmetric_basis(n, field=field) - (Qs, T) = _multiplication_table_from_matrix_basis(S) + Qs = _multiplication_table_from_matrix_basis(S) - fdeja = super(RealSymmetricEJA, cls) - return fdeja.__classcall_private__(cls, - field, - Qs, - rank=n, - natural_basis=T) + fdeja = super(RealSymmetricEJA, self) + return fdeja.__init__(field, + Qs, + rank=n, + natural_basis=S) def inner_product(self, x, y): return _matrix_ip(x,y) @@ -1097,12 +1152,12 @@ class ComplexHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra): TESTS: - The degree of this algebra is `n^2`:: + The dimension of this algebra is `n^2`:: sage: set_random_seed() sage: n = ZZ.random_element(1,5) sage: J = ComplexHermitianEJA(n) - sage: J.degree() == n^2 + sage: J.dimension() == n^2 True The Jordan multiplication is what we think it is:: @@ -1122,17 +1177,16 @@ class ComplexHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra): True """ - @staticmethod - def __classcall_private__(cls, n, field=QQ): + def __init__(self, n, field=QQ): S = _complex_hermitian_basis(n) - (Qs, T) = _multiplication_table_from_matrix_basis(S) + Qs = _multiplication_table_from_matrix_basis(S) + + fdeja = super(ComplexHermitianEJA, self) + return fdeja.__init__(field, + Qs, + rank=n, + natural_basis=S) - fdeja = super(ComplexHermitianEJA, cls) - return fdeja.__classcall_private__(cls, - field, - Qs, - rank=n, - natural_basis=T) def inner_product(self, x, y): # Since a+bi on the diagonal is represented as @@ -1158,12 +1212,12 @@ class QuaternionHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra): TESTS: - The degree of this algebra is `n^2`:: + The dimension of this algebra is `n^2`:: sage: set_random_seed() sage: n = ZZ.random_element(1,5) sage: J = QuaternionHermitianEJA(n) - sage: J.degree() == 2*(n^2) - n + sage: J.dimension() == 2*(n^2) - n True The Jordan multiplication is what we think it is:: @@ -1183,17 +1237,15 @@ class QuaternionHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra): True """ - @staticmethod - def __classcall_private__(cls, n, field=QQ): + def __init__(self, n, field=QQ): S = _quaternion_hermitian_basis(n) - (Qs, T) = _multiplication_table_from_matrix_basis(S) + Qs = _multiplication_table_from_matrix_basis(S) - fdeja = super(QuaternionHermitianEJA, cls) - return fdeja.__classcall_private__(cls, - field, - Qs, - rank=n, - natural_basis=T) + fdeja = super(QuaternionHermitianEJA, self) + return fdeja.__init__(field, + Qs, + rank=n, + natural_basis=S) def inner_product(self, x, y): # Since a+bi+cj+dk on the diagonal is represented as @@ -1241,26 +1293,28 @@ class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra): 0 """ - @staticmethod - def __classcall_private__(cls, n, field=QQ): - Qs = [] - id_matrix = matrix.identity(field, n) - for i in xrange(n): - ei = id_matrix.column(i) - Qi = matrix.zero(field, n) - Qi.set_row(0, ei) - Qi.set_column(0, ei) - Qi += matrix.diagonal(n, [ei[0]]*n) - # The addition of the diagonal matrix adds an extra ei[0] in the - # upper-left corner of the matrix. - Qi[0,0] = Qi[0,0] * ~field(2) - Qs.append(Qi) + def __init__(self, n, field=QQ): + 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:] + # z = x*y + z0 = x.inner_product(y) + zbar = y0*xbar + x0*ybar + z = V([z0] + zbar.list()) + mult_table[i][j] = z # The rank of the spin algebra is two, unless we're in a # one-dimensional ambient space (because the rank is bounded by # the ambient dimension). - fdeja = super(JordanSpinEJA, cls) - return fdeja.__classcall_private__(cls, field, Qs, rank=min(n,2)) + fdeja = super(JordanSpinEJA, self) + return fdeja.__init__(field, mult_table, rank=min(n,2)) def inner_product(self, x, y): return _usual_ip(x,y)