X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=mjo%2Feja%2Feja_algebra.py;h=3b5828fed55098dfd3f4c95bf2354ac4e38efc87;hb=0fc6cf97abf8e091787ebae4e9cd60534ebdbc32;hp=8ca501d6ba944a539901d384a78a5c8a0f45a7ca;hpb=3f8818f3eca0b2ea388130ff805875012cf902cb;p=sage.d.git diff --git a/mjo/eja/eja_algebra.py b/mjo/eja/eja_algebra.py index 8ca501d..3b5828f 100644 --- a/mjo/eja/eja_algebra.py +++ b/mjo/eja/eja_algebra.py @@ -3,6 +3,17 @@ Euclidean Jordan Algebras. These are formally-real Jordan Algebras; 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 @@ -14,7 +25,6 @@ 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, @@ -216,25 +226,6 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): coords = W.coordinate_vector(_mat2vec(elt)) return self.from_vector(coords) - @staticmethod - def _max_random_instance_size(): - """ - Return an integer "size" that is an upper bound on the size of - this algebra when it is used in a random test - case. Unfortunately, the term "size" is 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. - """ - raise NotImplementedError - def _repr_(self): """ Return a string representation of ``self``. @@ -588,6 +579,16 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): 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. @@ -832,16 +833,6 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): return tuple( self.random_element(thorough) for idx in range(count) ) - @classmethod - def random_instance(cls, field=AA, **kwargs): - """ - Return a random instance of this type of algebra. - - Beware, this will crash for "most instances" because the - constructor below looks wrong. - """ - n = ZZ.random_element(cls._max_random_instance_size() + 1) - return cls(n, field, **kwargs) @cached_method def _charpoly_coefficients(self): @@ -1003,32 +994,6 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): Element = FiniteDimensionalEuclideanJordanAlgebraElement - -def random_eja(field=AA): - """ - Return a "random" finite-dimensional Euclidean Jordan Algebra. - - SETUP:: - - sage: from mjo.eja.eja_algebra import random_eja - - TESTS:: - - sage: random_eja() - Euclidean Jordan algebra of dimension... - - """ - classname = choice([TrivialEJA, - HadamardEJA, - JordanSpinEJA, - RealSymmetricEJA, - ComplexHermitianEJA, - QuaternionHermitianEJA]) - return classname.random_instance(field=field) - - - - class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra): r""" Algebras whose basis consists of vectors with rational @@ -1088,12 +1053,72 @@ class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebr return tuple(map(lambda x: x.change_ring(self.base_ring()), a)) -class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra): +class ConcreteEuclideanJordanAlgebra: + r""" + A class for the Euclidean Jordan algebras that we know by name. + + These are the Jordan algebras whose basis, multiplication table, + rank, and so on are known a priori. More to the point, they are + the Euclidean Jordan algebras for which we are able to conjure up + a "random instance." + + SETUP:: + + sage: from mjo.eja.eja_algebra import ConcreteEuclideanJordanAlgebra + + TESTS: + + Our natural basis is normalized with respect to the natural 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 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: J = ConcreteEuclideanJordanAlgebra.random_instance() + sage: x = J.random_element() + sage: x.operator().matrix().is_symmetric() + True + + """ + @staticmethod 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 + """ + Return an integer "size" that is an upper bound on the size of + this algebra when it is used in a random test + case. Unfortunately, the term "size" is ambiguous -- when + dealing with `R^n` under either the Hadamard or Jordan spin + product, the "size" refers to the dimension `n`. When dealing + with a matrix algebra (real symmetric or complex/quaternion + Hermitian), it refers to the size of the matrix, which is far + less than the dimension of the underlying vector space. + + This method must be implemented in each subclass. + """ + raise NotImplementedError + + @classmethod + def random_instance(cls, field=AA, **kwargs): + """ + Return a random instance of this type of algebra. + + This method should be implemented in each subclass. + """ + from sage.misc.prandom import choice + eja_class = choice(cls.__subclasses__()) + return eja_class.random_instance(field) + + +class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra): def __init__(self, field, basis, normalize_basis=True, **kwargs): """ @@ -1108,7 +1133,8 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra): # time to ensure that it isn't a generator expression. basis = tuple(basis) - if len(basis) > 1 and normalize_basis: + algebra_dim = len(basis) + if algebra_dim > 1 and normalize_basis: # We'll need sqrt(2) to normalize the basis, and this # winds up in the multiplication table, so the whole # algebra needs to be over the field extension. @@ -1129,6 +1155,14 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra): natural_basis=basis, **kwargs) + if algebra_dim == 0: + self.one.set_cache(self.zero()) + else: + n = basis[0].nrows() + # The identity wrt (A,B) -> (AB + BA)/2 is independent of the + # details of this algebra. + self.one.set_cache(self(matrix.identity(field,n))) + @cached_method def _charpoly_coefficients(self): @@ -1268,7 +1302,8 @@ class RealMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): return M -class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra): +class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra, + ConcreteEuclideanJordanAlgebra): """ The rank-n simple EJA consisting of real symmetric n-by-n matrices, the usual symmetric Jordan product, and the trace inner @@ -1327,25 +1362,6 @@ class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra): sage: RealSymmetricEJA(3, prefix='q').gens() (q0, q1, q2, q3, q4, q5) - Our natural basis is normalized with respect to the natural inner - product unless we specify otherwise:: - - sage: set_random_seed() - sage: J = RealSymmetricEJA.random_instance() - sage: all( b.norm() == 1 for b in J.gens() ) - True - - Since our natural basis is normalized with respect to the natural - inner product, and since we know that this algebra is an EJA, any - left-multiplication operator's matrix will be symmetric because - natural->EJA basis representation is an isometry and within the EJA - the operator is self-adjoint by the Jordan axiom:: - - sage: set_random_seed() - sage: x = RealSymmetricEJA.random_instance().random_element() - sage: x.operator().matrix().is_symmetric() - True - We can construct the (trivial) algebra of rank zero:: sage: RealSymmetricEJA(0) @@ -1388,6 +1404,13 @@ class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra): 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) @@ -1431,8 +1454,7 @@ class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): Embedding is a homomorphism (isomorphism, in fact):: sage: set_random_seed() - sage: n_max = ComplexMatrixEuclideanJordanAlgebra._max_random_instance_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) @@ -1557,7 +1579,8 @@ class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): return RealMatrixEuclideanJordanAlgebra.natural_inner_product(X,Y)/2 -class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra): +class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra, + ConcreteEuclideanJordanAlgebra): """ The rank-n simple EJA consisting of complex Hermitian n-by-n matrices over the real numbers, the usual symmetric Jordan product, @@ -1608,25 +1631,6 @@ class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra): sage: ComplexHermitianEJA(2, prefix='z').gens() (z0, z1, z2, z3) - Our natural basis is normalized with respect to the natural inner - product unless we specify otherwise:: - - sage: set_random_seed() - sage: J = ComplexHermitianEJA.random_instance() - sage: all( b.norm() == 1 for b in J.gens() ) - True - - Since our natural basis is normalized with respect to the natural - inner product, and since we know that this algebra is an EJA, any - left-multiplication operator's matrix will be symmetric because - natural->EJA basis representation is an isometry and within the EJA - the operator is self-adjoint by the Jordan axiom:: - - sage: set_random_seed() - sage: x = ComplexHermitianEJA.random_instance().random_element() - sage: x.operator().matrix().is_symmetric() - True - We can construct the (trivial) algebra of rank zero:: sage: ComplexHermitianEJA(0) @@ -1696,6 +1700,17 @@ class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra): **kwargs) self.rank.set_cache(n) + @staticmethod + def _max_random_instance_size(): + return 3 # Dimension 9 + + @classmethod + def random_instance(cls, field=AA, **kwargs): + """ + Return a random instance of this type of algebra. + """ + n = ZZ.random_element(cls._max_random_instance_size() + 1) + return cls(n, field, **kwargs) class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): @staticmethod @@ -1727,8 +1742,7 @@ class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): Embedding is a homomorphism (isomorphism, in fact):: sage: set_random_seed() - sage: n_max = QuaternionMatrixEuclideanJordanAlgebra._max_random_instance_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) @@ -1860,8 +1874,9 @@ class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): return RealMatrixEuclideanJordanAlgebra.natural_inner_product(X,Y)/4 -class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra): - """ +class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra, + ConcreteEuclideanJordanAlgebra): + r""" The rank-n simple EJA consisting of self-adjoint n-by-n quaternion matrices, the usual symmetric Jordan product, and the real-part-of-trace inner product. It has dimension `2n^2 - n` over @@ -1911,25 +1926,6 @@ class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra): sage: QuaternionHermitianEJA(2, prefix='a').gens() (a0, a1, a2, a3, a4, a5) - Our natural basis is normalized with respect to the natural inner - product unless we specify otherwise:: - - sage: set_random_seed() - sage: J = QuaternionHermitianEJA.random_instance() - sage: all( b.norm() == 1 for b in J.gens() ) - True - - Since our natural basis is normalized with respect to the natural - inner product, and since we know that this algebra is an EJA, any - left-multiplication operator's matrix will be symmetric because - natural->EJA basis representation is an isometry and within the EJA - the operator is self-adjoint by the Jordan axiom:: - - sage: set_random_seed() - sage: x = QuaternionHermitianEJA.random_instance().random_element() - sage: x.operator().matrix().is_symmetric() - True - We can construct the (trivial) algebra of rank zero:: sage: QuaternionHermitianEJA(0) @@ -2000,8 +1996,24 @@ class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra): **kwargs) self.rank.set_cache(n) + @staticmethod + def _max_random_instance_size(): + r""" + The maximum rank of a random QuaternionHermitianEJA. + """ + return 2 # Dimension 6 + + @classmethod + def random_instance(cls, field=AA, **kwargs): + """ + Return a random instance of this type of algebra. + """ + n = ZZ.random_element(cls._max_random_instance_size() + 1) + return cls(n, field, **kwargs) + -class HadamardEJA(RationalBasisEuclideanJordanAlgebra): +class HadamardEJA(RationalBasisEuclideanJordanAlgebra, + ConcreteEuclideanJordanAlgebra): """ Return the Euclidean Jordan Algebra corresponding to the set `R^n` under the Hadamard product. @@ -2052,10 +2064,27 @@ class HadamardEJA(RationalBasisEuclideanJordanAlgebra): **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) + + def inner_product(self, x, y): """ Faster to reimplement than to use natural representations. @@ -2081,7 +2110,8 @@ class HadamardEJA(RationalBasisEuclideanJordanAlgebra): return x.to_vector().inner_product(y.to_vector()) -class BilinearFormEJA(RationalBasisEuclideanJordanAlgebra): +class BilinearFormEJA(RationalBasisEuclideanJordanAlgebra, + 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 = @@ -2111,6 +2141,20 @@ class BilinearFormEJA(RationalBasisEuclideanJordanAlgebra): 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:: @@ -2151,7 +2195,7 @@ class BilinearFormEJA(RationalBasisEuclideanJordanAlgebra): n = B.nrows() if not B.is_positive_definite(): - raise TypeError("matrix B is not positive-definite") + raise ValueError("bilinear form is not positive-definite") V = VectorSpace(field, n) mult_table = [[V.zero() for j in range(n)] for i in range(n)] @@ -2177,8 +2221,16 @@ class BilinearFormEJA(RationalBasisEuclideanJordanAlgebra): **kwargs) self.rank.set_cache(min(n,2)) + if n == 0: + self.one.set_cache( self.zero() ) + else: + self.one.set_cache( self.monomial(0) ) + @staticmethod def _max_random_instance_size(): + r""" + The maximum dimension of a random BilinearFormEJA. + """ return 5 @classmethod @@ -2187,9 +2239,8 @@ class BilinearFormEJA(RationalBasisEuclideanJordanAlgebra): 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) + if n.is_zero(): + B = matrix.identity(field, n) return cls(B, field, **kwargs) B11 = matrix.identity(field,1) @@ -2292,8 +2343,26 @@ class JordanSpinEJA(BilinearFormEJA): B = matrix.identity(field, n) super(JordanSpinEJA, self).__init__(B, field, **kwargs) + @staticmethod + def _max_random_instance_size(): + r""" + The maximum dimension of a random JordanSpinEJA. + """ + return 5 -class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra): + @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, + ConcreteEuclideanJordanAlgebra): """ The trivial Euclidean Jordan algebra consisting of only a zero element. @@ -2331,6 +2400,7 @@ class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra): # The rank is zero using my definition, namely the dimension of the # largest subalgebra generated by any element. self.rank.set_cache(0) + self.one.set_cache( self.zero() ) @classmethod def random_instance(cls, field=AA, **kwargs): @@ -2519,3 +2589,7 @@ class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra): y2 = pi_right(y) return (x1.inner_product(y1) + x2.inner_product(y2)) + + + +random_eja = ConcreteEuclideanJordanAlgebra.random_instance