X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=mjo%2Feja%2Feja_algebra.py;h=8ca501d6ba944a539901d384a78a5c8a0f45a7ca;hb=3f8818f3eca0b2ea388130ff805875012cf902cb;hp=b9389da25e69722f04633c8324998d8d4834608a;hpb=5f522a9ac121a1229c0a9c62244b18ca67e65b12;p=sage.d.git diff --git a/mjo/eja/eja_algebra.py b/mjo/eja/eja_algebra.py index b9389da..8ca501d 100644 --- a/mjo/eja/eja_algebra.py +++ b/mjo/eja/eja_algebra.py @@ -217,7 +217,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): return self.from_vector(coords) @staticmethod - def _max_test_case_size(): + def _max_random_instance_size(): """ Return an integer "size" that is an upper bound on the size of this algebra when it is used in a random test @@ -233,7 +233,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): interpreted to be far less than the dimension) should override with a smaller number. """ - return 5 + raise NotImplementedError def _repr_(self): """ @@ -599,19 +599,20 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): # 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): @@ -839,12 +840,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): 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) + n = ZZ.random_element(cls._max_random_instance_size() + 1) return cls(n, field, **kwargs) @cached_method @@ -1094,7 +1090,7 @@ class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebr class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra): @staticmethod - def _max_test_case_size(): + def _max_random_instance_size(): # Play it safe, since this will be squared and the underlying # field can have dimension 4 (quaternions) too. return 2 @@ -1306,7 +1302,7 @@ class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra): 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 @@ -1389,7 +1385,7 @@ class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra): @staticmethod - def _max_test_case_size(): + def _max_random_instance_size(): return 4 # Dimension 10 @@ -1435,7 +1431,7 @@ class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): Embedding is a homomorphism (isomorphism, in fact):: sage: set_random_seed() - sage: n_max = ComplexMatrixEuclideanJordanAlgebra._max_test_case_size() + sage: n_max = ComplexMatrixEuclideanJordanAlgebra._max_random_instance_size() sage: n = ZZ.random_element(n_max) sage: F = QuadraticField(-1, 'I') sage: X = random_matrix(F, n) @@ -1587,7 +1583,7 @@ class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra): 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 @@ -1731,7 +1727,7 @@ class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): Embedding is a homomorphism (isomorphism, in fact):: sage: set_random_seed() - sage: n_max = QuaternionMatrixEuclideanJordanAlgebra._max_test_case_size() + sage: n_max = QuaternionMatrixEuclideanJordanAlgebra._max_random_instance_size() sage: n = ZZ.random_element(n_max) sage: Q = QuaternionAlgebra(QQ,-1,-1) sage: X = random_matrix(Q, n) @@ -1890,7 +1886,7 @@ class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra): 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 @@ -2056,6 +2052,10 @@ class HadamardEJA(RationalBasisEuclideanJordanAlgebra): **kwargs) self.rank.set_cache(n) + @staticmethod + def _max_random_instance_size(): + return 5 + def inner_product(self, x, y): """ Faster to reimplement than to use natural representations. @@ -2085,10 +2085,15 @@ class BilinearFormEJA(RationalBasisEuclideanJordanAlgebra): 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 + , 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``. + (,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:: @@ -2100,7 +2105,8 @@ class BilinearFormEJA(RationalBasisEuclideanJordanAlgebra): 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 @@ -2109,7 +2115,8 @@ class BilinearFormEJA(RationalBasisEuclideanJordanAlgebra): 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 {} @@ -2121,8 +2128,11 @@ class BilinearFormEJA(RationalBasisEuclideanJordanAlgebra): 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) 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())) @@ -2136,11 +2146,12 @@ class BilinearFormEJA(RationalBasisEuclideanJordanAlgebra): 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): + self._B = B + n = B.nrows() + + if not B.is_positive_definite(): + raise TypeError("matrix B is not positive-definite") V = VectorSpace(field, n) mult_table = [[V.zero() for j in range(n)] for i in range(n)] @@ -2152,7 +2163,7 @@ class BilinearFormEJA(RationalBasisEuclideanJordanAlgebra): xbar = x[1:] y0 = y[0] ybar = y[1:] - z0 = x0*y0 + (self._B*xbar).inner_product(ybar) + z0 = (B*x).inner_product(y) zbar = y0*xbar + x0*ybar z = V([z0] + zbar.list()) mult_table[i][j] = z @@ -2166,6 +2177,35 @@ class BilinearFormEJA(RationalBasisEuclideanJordanAlgebra): **kwargs) self.rank.set_cache(min(n,2)) + @staticmethod + def _max_random_instance_size(): + return 5 + + @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 == 0: + # Special case needed since we use (n-1) below. + B = matrix.identity(field, 0) + return cls(B, field, **kwargs) + + 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 + + from sage.matrix.special import block_matrix + B = block_matrix(2,2, [ [B11, ZZ(0) ], + [ZZ(0), B22 ] ]) + + return cls(B, field, **kwargs) + def inner_product(self, x, y): r""" Half of the trace inner product. @@ -2185,21 +2225,15 @@ class BilinearFormEJA(RationalBasisEuclideanJordanAlgebra): paper:: 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: J = BilinearFormEJA.random_instance() + sage: n = J.dimension() sage: x = J.random_element() sage: y = J.random_element() - sage: x.inner_product(y) == (x*y).trace()/2 + sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2) True """ - 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 (self._B*x.to_vector()).inner_product(y.to_vector()) class JordanSpinEJA(BilinearFormEJA): @@ -2255,7 +2289,8 @@ class JordanSpinEJA(BilinearFormEJA): def __init__(self, n, field=AA, **kwargs): # This is a special case of the BilinearFormEJA with the identity # matrix as its bilinear form. - super(JordanSpinEJA, self).__init__(n, field, **kwargs) + B = matrix.identity(field, n) + super(JordanSpinEJA, self).__init__(B, field, **kwargs) class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra): @@ -2297,6 +2332,11 @@ class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra): # largest subalgebra generated by any element. self.rank.set_cache(0) + @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""" @@ -2324,6 +2364,7 @@ class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra): """ def __init__(self, J1, J2, field=AA, **kwargs): + self._factors = (J1, J2) n1 = J1.dimension() n2 = J2.dimension() n = n1+n2 @@ -2345,3 +2386,136 @@ class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra): 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() + n = J1.dimension() + pi_left = lambda x: J1.from_vector(x.to_vector()[:n]) + pi_right = lambda x: J2.from_vector(x.to_vector()[n:]) + 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 (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) + + """ + (J1,J2) = self.factors() + n = J1.dimension() + V_basis = self.vector_space().basis() + I1 = matrix.column(self.base_ring(), V_basis[:n]) + I2 = matrix.column(self.base_ring(), V_basis[n:]) + iota_left = lambda x: self.from_vector(I1*x.to_vector()) + iota_right = lambda x: self.from_vector(I2*+x.to_vector()) + 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) + sage: J2 = QuaternionHermitianEJA(2,QQ,normalize_basis=False) + sage: J = DirectSumEJA(J1,J2) + sage: x1 = J1.one() + sage: x2 = x1 + sage: y1 = J2.one() + sage: y2 = y1 + sage: x1.inner_product(x2) + 3 + sage: y1.inner_product(y2) + 2 + sage: J.one().inner_product(J.one()) + 5 + + """ + (pi_left, pi_right) = self.projections() + x1 = pi_left(x) + x2 = pi_right(x) + y1 = pi_left(y) + y2 = pi_right(y) + + return (x1.inner_product(y1) + x2.inner_product(y2))