X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=mjo%2Fcone%2Fcone.py;h=cbbfe9b692f3bdf36f68f202deedd467c684cef0;hb=115ecffcebcd0b7f86358aa38ddb52f1115e966c;hp=2296e3fc010db71091e530b211c73ada05962f81;hpb=3b659b1d0440daf3ff7bd8cf3cf53f90523a1609;p=sage.d.git diff --git a/mjo/cone/cone.py b/mjo/cone/cone.py index 2296e3f..cbbfe9b 100644 --- a/mjo/cone/cone.py +++ b/mjo/cone/cone.py @@ -1,236 +1,350 @@ -# Sage doesn't load ~/.sage/init.sage during testing (sage -t), so we -# have to explicitly mangle our sitedir here so that "mjo.cone" -# resolves. -from os.path import abspath -from site import addsitedir -addsitedir(abspath('../../')) - from sage.all import * +from sage.geometry.cone import is_Cone - -def random_cone(min_dim=None, max_dim=None, min_rays=None, max_rays=None): +def is_positive_on(L,K): r""" - Generate a random rational convex polyhedral cone. + Determine whether or not ``L`` is positive on ``K``. + + We say that ``L`` is positive on a closed convex cone ``K`` if + `L\left\lparen x \right\rparen` belongs to ``K`` for all `x` in + ``K``. This property need only be checked for generators of ``K``. - Lower and upper bounds may be provided for both the dimension of the - ambient space and the number of generating rays of the cone. Any - parameters left unspecified will be chosen randomly. + To reliably check whether or not ``L`` is positive, its base ring + must be either exact (for example, the rationals) or ``SR``. An + exact ring is more reliable, but in some cases a matrix whose + entries contain symbolic constants like ``e`` and ``pi`` will work. INPUT: - - ``min_dim`` (default: random) -- The minimum dimension of the ambient - lattice. + - ``L`` -- A matrix over either an exact ring or ``SR``. - - ``max_dim`` (default: random) -- The maximum dimension of the ambient - lattice. + - ``K`` -- A polyhedral closed convex cone. - - ``min_rays`` (default: random) -- The minimum number of generating rays - of the cone. + OUTPUT: - - ``max_rays`` (default: random) -- The maximum number of generating rays - of the cone. + If the base ring of ``L`` is exact, then ``True`` will be returned if + and only if ``L`` is positive on ``K``. - OUTPUT: + If the base ring of ``L`` is ``SR``, then the situation is more + complicated: - A new, randomly generated cone. + - ``True`` will be returned if it can be proven that ``L`` + is positive on ``K``. + - ``False`` will be returned if it can be proven that ``L`` + is not positive on ``K``. + - ``False`` will also be returned if we can't decide; specifically + if we arrive at a symbolic inequality that cannot be resolved. + + EXAMPLES: + + Nonnegative matrices are positive operators on the nonnegative + orthant:: + + sage: K = Cone([(1,0,0),(0,1,0),(0,0,1)]) + sage: L = random_matrix(QQ,3).apply_map(abs) + sage: is_positive_on(L,K) + True TESTS: - It's hard to test the output of a random process, but we can at - least make sure that we get a cone back:: + The identity operator is always positive:: - sage: from sage.geometry.cone import is_Cone - sage: K = random_cone() - sage: is_Cone(K) # long time + sage: set_random_seed() + sage: K = random_cone(max_ambient_dim=8) + sage: L = identity_matrix(K.lattice_dim()) + sage: is_positive_on(L,K) True - """ + The "zero" operator is always positive:: - def random_min_max(l,u): - r""" - We need to handle four cases to prevent us from doing - something stupid like having an upper bound that's lower than - our lower bound. And we would need to repeat all of that logic - for the dimension/rays, so we consolidate it here. - """ - if l is None and u is None: - # They're both random, just return a random nonnegative - # integer. - return ZZ.random_element().abs() - - if l is not None and u is not None: - # Both were specified. Again, just make up a number and - # return it. If the user wants to give us u < l then he - # can have an exception. - return ZZ.random_element(l,u) - - if l is not None and u is None: - # In this case, we're generating the upper bound randomly - # GIVEN A LOWER BOUND. So we add a random nonnegative - # integer to the given lower bound. - u = l + ZZ.random_element().abs() - return ZZ.random_element(l,u) - - # Here we must be in the only remaining case, where we are - # given an upper bound but no lower bound. We might as well - # use zero. - return ZZ.random_element(0,u) - - d = random_min_max(min_dim, max_dim) - r = random_min_max(min_rays, max_rays) - - L = ToricLattice(d) - rays = [L.random_element() for i in range(0,r)] - - # We pass the lattice in case there are no rays. - return Cone(rays, lattice=L) - - -def lyapunov_rank(K): - r""" - Compute the Lyapunov (or bilinearity) rank of this cone. + sage: K = random_cone(max_ambient_dim=8) + sage: R = K.lattice().vector_space().base_ring() + sage: L = zero_matrix(R, K.lattice_dim()) + sage: is_positive_on(L,K) + True - The Lyapunov rank of a cone can be thought of in (mainly) two ways: + Everything in ``K.positive_operators_gens()`` should be + positive on ``K``:: - 1. The dimension of the Lie algebra of the automorphism group of the - cone. + sage: K = random_cone(max_ambient_dim=5) + sage: all([ is_positive_on(L,K) # long time + ....: for L in K.positive_operators_gens() ]) # long time + True + sage: all([ is_positive_on(L.change_ring(SR),K) # long time + ....: for L in K.positive_operators_gens() ]) # long time + True - 2. The dimension of the linear space of all Lyapunov-like - transformations on the cone. + """ + if not is_Cone(K): + raise TypeError('K must be a Cone') + if not L.base_ring().is_exact() and not L.base_ring() is SR: + raise ValueError('base ring of operator L is neither SR nor exact') + + if L.base_ring().is_exact(): + # This should be way faster than computing the dual and + # checking a bunch of inequalities, but it doesn't work if + # ``L*x`` is symbolic. For example, ``e in Cone([(1,)])`` + # is true, but returns ``False``. + return all([ L*x in K for x in K ]) + else: + # Fall back to inequality-checking when the entries of ``L`` + # might be symbolic. + return all([ s*(L*x) >= 0 for x in K for s in K.dual() ]) + + +def is_cross_positive_on(L,K): + r""" + Determine whether or not ``L`` is cross-positive on ``K``. + + We say that ``L`` is cross-positive on a closed convex cone``K`` if + `\left\langle L\left\lparenx\right\rparen,s\right\rangle \ge 0` for + all pairs `\left\langle x,s \right\rangle` in the complementarity + set of ``K``. This property need only be checked for generators of + ``K`` and its dual. + + To reliably check whether or not ``L`` is cross-positive, its base + ring must be either exact (for example, the rationals) or ``SR``. An + exact ring is more reliable, but in some cases a matrix whose + entries contain symbolic constants like ``e`` and ``pi`` will work. INPUT: - A closed, convex polyhedral cone. + - ``L`` -- A matrix over either an exact ring or ``SR``. + + - ``K`` -- A polyhedral closed convex cone. OUTPUT: - An integer representing the Lyapunov rank of the cone. If the - dimension of the ambient vector space is `n`, then the Lyapunov rank - will be between `1` and `n` inclusive; however a rank of `n-1` is - not possible (see the first reference). + If the base ring of ``L`` is exact, then ``True`` will be returned if + and only if ``L`` is cross-positive on ``K``. - .. note:: + If the base ring of ``L`` is ``SR``, then the situation is more + complicated: - In the references, the cones are always assumed to be proper. We - do not impose this restriction. + - ``True`` will be returned if it can be proven that ``L`` + is cross-positive on ``K``. + - ``False`` will be returned if it can be proven that ``L`` + is not cross-positive on ``K``. + - ``False`` will also be returned if we can't decide; specifically + if we arrive at a symbolic inequality that cannot be resolved. - .. seealso:: + EXAMPLES: - :meth:`is_proper` + The identity operator is always cross-positive:: - ALGORITHM: + sage: set_random_seed() + sage: K = random_cone(max_ambient_dim=8) + sage: L = identity_matrix(K.lattice_dim()) + sage: is_cross_positive_on(L,K) + True - The codimension formula from the second reference is used. We find - all pairs `(x,s)` in the complementarity set of `K` such that `x` - and `s` are rays of our cone. It is known that these vectors are - sufficient to apply the codimension formula. Once we have all such - pairs, we "brute force" the codimension formula by finding all - linearly-independent `xs^{T}`. + The "zero" operator is always cross-positive:: - REFERENCES: + sage: K = random_cone(max_ambient_dim=8) + sage: R = K.lattice().vector_space().base_ring() + sage: L = zero_matrix(R, K.lattice_dim()) + sage: is_cross_positive_on(L,K) + True - 1. M.S. Gowda and J. Tao. On the bilinearity rank of a proper cone - and Lyapunov-like transformations, Mathematical Programming, 147 - (2014) 155-170. + TESTS: - 2. G. Rudolf, N. Noyan, D. Papp, and F. Alizadeh, Bilinear - optimality constraints for the cone of positive polynomials, - Mathematical Programming, Series B, 129 (2011) 5-31. + Everything in ``K.cross_positive_operators_gens()`` should be + cross-positive on ``K``:: - EXAMPLES: + sage: K = random_cone(max_ambient_dim=5) + sage: all([ is_cross_positive_on(L,K) # long time + ....: for L in K.cross_positive_operators_gens() ]) # long time + True + sage: all([ is_cross_positive_on(L.change_ring(SR),K) # long time + ....: for L in K.cross_positive_operators_gens() ]) # long time + True + + """ + if not is_Cone(K): + raise TypeError('K must be a Cone') + if not L.base_ring().is_exact() and not L.base_ring() is SR: + raise ValueError('base ring of operator L is neither SR nor exact') - The nonnegative orthant in `\mathbb{R}^{n}` always has rank `n`:: + return all([ s*(L*x) >= 0 + for (x,s) in K.discrete_complementarity_set() ]) - sage: positives = Cone([(1,)]) - sage: lyapunov_rank(positives) - 1 - sage: quadrant = Cone([(1,0), (0,1)]) - sage: lyapunov_rank(quadrant) - 2 - sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)]) - sage: lyapunov_rank(octant) - 3 +def is_Z_on(L,K): + r""" + Determine whether or not ``L`` is a Z-operator on ``K``. + + We say that ``L`` is a Z-operator on a closed convex cone``K`` if + `\left\langle L\left\lparenx\right\rparen,s\right\rangle \le 0` for + all pairs `\left\langle x,s \right\rangle` in the complementarity + set of ``K``. It is known that this property need only be checked + for generators of ``K`` and its dual. - The `L^{3}_{1}` cone is known to have a Lyapunov rank of one:: + A matrix is a Z-operator on ``K`` if and only if its negation is a + cross-positive operator on ``K``. - sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)]) - sage: lyapunov_rank(L31) - 1 + To reliably check whether or not ``L`` is a Z operator, its base + ring must be either exact (for example, the rationals) or ``SR``. An + exact ring is more reliable, but in some cases a matrix whose + entries contain symbolic constants like ``e`` and ``pi`` will work. - Likewise for the `L^{3}_{\infty}` cone:: + INPUT: - sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)]) - sage: lyapunov_rank(L3infty) - 1 + - ``L`` -- A matrix over either an exact ring or ``SR``. - The Lyapunov rank should be additive on a product of cones:: + - ``K`` -- A polyhedral closed convex cone. - sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)]) - sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)]) - sage: K = L31.cartesian_product(octant) - sage: lyapunov_rank(K) == lyapunov_rank(L31) + lyapunov_rank(octant) - True + OUTPUT: + + If the base ring of ``L`` is exact, then ``True`` will be returned if + and only if ``L`` is a Z-operator on ``K``. - Two isomorphic cones should have the same Lyapunov rank. The cone - ``K`` in the following example is isomorphic to the nonnegative - octant in `\mathbb{R}^{3}`:: + If the base ring of ``L`` is ``SR``, then the situation is more + complicated: - sage: K = Cone([(1,2,3), (-1,1,0), (1,0,6)]) - sage: lyapunov_rank(K) - 3 + - ``True`` will be returned if it can be proven that ``L`` + is a Z-operator on ``K``. + - ``False`` will be returned if it can be proven that ``L`` + is not a Z-operator on ``K``. + - ``False`` will also be returned if we can't decide; specifically + if we arrive at a symbolic inequality that cannot be resolved. - The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K`` - itself:: + EXAMPLES: - sage: K = Cone([(2,2,4), (-1,9,0), (2,0,6)]) - sage: lyapunov_rank(K) == lyapunov_rank(K.dual()) + The identity operator is always a Z-operator:: + + sage: set_random_seed() + sage: K = random_cone(max_ambient_dim=8) + sage: L = identity_matrix(K.lattice_dim()) + sage: is_Z_on(L,K) + True + + The "zero" operator is always a Z-operator:: + + sage: K = random_cone(max_ambient_dim=8) + sage: R = K.lattice().vector_space().base_ring() + sage: L = zero_matrix(R, K.lattice_dim()) + sage: is_Z_on(L,K) True TESTS: - The Lyapunov rank should be additive on a product of cones:: + Everything in ``K.Z_operators_gens()`` should be a Z-operator + on ``K``:: - sage: K1 = random_cone(0,10,0,10) - sage: K2 = random_cone(0,10,0,10) - sage: K = K1.cartesian_product(K2) - sage: lyapunov_rank(K) == lyapunov_rank(K1) + lyapunov_rank(K2) + sage: K = random_cone(max_ambient_dim=5) + sage: all([ is_Z_on(L,K) # long time + ....: for L in K.Z_operators_gens() ]) # long time True + sage: all([ is_Z_on(L.change_ring(SR),K) # long time + ....: for L in K.Z_operators_gens() ]) # long time + True + + """ + return is_cross_positive_on(-L,K) + + +def is_lyapunov_like_on(L,K): + r""" + Determine whether or not ``L`` is Lyapunov-like on ``K``. + + We say that ``L`` is Lyapunov-like on a closed convex cone ``K`` if + `\left\langle L\left\lparenx\right\rparen,s\right\rangle = 0` for + all pairs `\left\langle x,s \right\rangle` in the complementarity + set of ``K``. This property need only be checked for generators of + ``K`` and its dual. + + To reliably check whether or not ``L`` is Lyapunov-like, its base + ring must be either exact (for example, the rationals) or ``SR``. An + exact ring is more reliable, but in some cases a matrix whose + entries contain symbolic constants like ``e`` and ``pi`` will work. + + INPUT: - The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K`` - itself:: + - ``L`` -- A matrix over either an exact ring or ``SR``. - sage: K = random_cone(0,10,0,10) - sage: lyapunov_rank(K) == lyapunov_rank(K.dual()) + - ``K`` -- A polyhedral closed convex cone. + + OUTPUT: + + If the base ring of ``L`` is exact, then ``True`` will be returned if + and only if ``L`` is Lyapunov-like on ``K``. + + If the base ring of ``L`` is ``SR``, then the situation is more + complicated: + + - ``True`` will be returned if it can be proven that ``L`` + is Lyapunov-like on ``K``. + - ``False`` will be returned if it can be proven that ``L`` + is not Lyapunov-like on ``K``. + - ``False`` will also be returned if we can't decide; specifically + if we arrive at a symbolic inequality that cannot be resolved. + + EXAMPLES: + + Diagonal matrices are Lyapunov-like operators on the nonnegative + orthant:: + + sage: K = Cone([(1,0,0),(0,1,0),(0,0,1)]) + sage: L = diagonal_matrix(random_vector(QQ,3)) + sage: is_lyapunov_like_on(L,K) True - """ - V = K.lattice().vector_space() + TESTS: - xs = [V(x) for x in K.rays()] - ss = [V(s) for s in K.dual().rays()] + The identity operator is always Lyapunov-like:: - # WARNING: This isn't really C(K), it only contains the pairs - # (x,s) in C(K) where x,s are extreme in their respective cones. - C_of_K = [(x,s) for x in xs for s in ss if x.inner_product(s) == 0] + sage: set_random_seed() + sage: K = random_cone(max_ambient_dim=8) + sage: L = identity_matrix(K.lattice_dim()) + sage: is_lyapunov_like_on(L,K) + True - matrices = [x.column() * s.row() for (x,s) in C_of_K] + The "zero" operator is always Lyapunov-like:: - # Sage doesn't think matrices are vectors, so we have to convert - # our matrices to vectors explicitly before we can figure out how - # many are linearly-indepenedent. - # - # The space W has the same base ring as V, but dimension - # dim(V)^2. So it has the same dimension as the space of linear - # transformations on V. In other words, it's just the right size - # to create an isomorphism between it and our matrices. - W = VectorSpace(V.base_ring(), V.dimension()**2) + sage: K = random_cone(max_ambient_dim=8) + sage: R = K.lattice().vector_space().base_ring() + sage: L = zero_matrix(R, K.lattice_dim()) + sage: is_lyapunov_like_on(L,K) + True - def phi(m): - r""" - Convert a matrix to a vector isomorphically. - """ - return W(m.list()) + Everything in ``K.lyapunov_like_basis()`` should be Lyapunov-like + on ``K``:: - vectors = [phi(m) for m in matrices] + sage: K = random_cone(max_ambient_dim=5) + sage: all([ is_lyapunov_like_on(L,K) # long time + ....: for L in K.lyapunov_like_basis() ]) # long time + True + sage: all([ is_lyapunov_like_on(L.change_ring(SR),K) # long time + ....: for L in K.lyapunov_like_basis() ]) # long time + True - return (W.dimension() - W.span(vectors).rank()) + """ + if not is_Cone(K): + raise TypeError('K must be a Cone') + if not L.base_ring().is_exact() and not L.base_ring() is SR: + raise ValueError('base ring of operator L is neither SR nor exact') + + return all([ s*(L*x) == 0 + for (x,s) in K.discrete_complementarity_set() ]) + + +def LL_cone(K): + gens = K.lyapunov_like_basis() + L = ToricLattice(K.lattice_dim()**2) + return Cone([ g.list() for g in gens ], lattice=L, check=False) + +def Sigma_cone(K): + gens = K.cross_positive_operators_gens() + L = ToricLattice(K.lattice_dim()**2) + return Cone([ g.list() for g in gens ], lattice=L, check=False) + +def Z_cone(K): + gens = K.Z_operators_gens() + L = ToricLattice(K.lattice_dim()**2) + return Cone([ g.list() for g in gens ], lattice=L, check=False) + +def pi_cone(K1, K2=None): + if K2 is None: + K2 = K1 + gens = K1.positive_operators_gens(K2) + L = ToricLattice(K1.lattice_dim()*K2.lattice_dim()) + return Cone([ g.list() for g in gens ], lattice=L, check=False)