X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=mjo%2Fcone%2Fcone.py;h=cbbfe9b692f3bdf36f68f202deedd467c684cef0;hb=115ecffcebcd0b7f86358aa38ddb52f1115e966c;hp=23043381bca83acd6d7f17479ef4769980b9960f;hpb=012175f3a2591586099b4e955bb440958f2d63d7;p=sage.d.git diff --git a/mjo/cone/cone.py b/mjo/cone/cone.py index 2304338..cbbfe9b 100644 --- a/mjo/cone/cone.py +++ b/mjo/cone/cone.py @@ -1,141 +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 lyapunov_rank(K): +def is_positive_on(L,K): r""" - Compute the Lyapunov (or bilinearity) rank of this cone. - - The Lyapunov rank of a cone can be thought of in (mainly) two ways: + Determine whether or not ``L`` is positive on ``K``. - 1. The dimension of the Lie algebra of the automorphism group of the - cone. + 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``. - 2. The dimension of the linear space of all Lyapunov-like - transformations on the cone. + 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: - 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 positive 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 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: + + The identity operator is always positive:: + + 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:: + + 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 + + Everything in ``K.positive_operators_gens()`` should be + positive on ``K``:: + + 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 + + """ + 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``. - .. note:: + 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. - In the references, the cones are always assumed to be proper. We - do not impose this restriction. + 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. - .. seealso:: + INPUT: - :meth:`is_proper` + - ``L`` -- A matrix over either an exact ring or ``SR``. - ALGORITHM: + - ``K`` -- A polyhedral closed convex cone. - 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}`. + OUTPUT: - REFERENCES: + If the base ring of ``L`` is exact, then ``True`` will be returned if + and only if ``L`` is cross-positive on ``K``. - 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. + If the base ring of ``L`` is ``SR``, then the situation is more + complicated: - 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. + - ``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. EXAMPLES: - The nonnegative orthant in `\mathbb{R}^{n}` always has rank `n`:: + The identity operator is always cross-positive:: + + 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 "zero" operator is always cross-positive:: - 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 + 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 + + TESTS: + + Everything in ``K.cross_positive_operators_gens()`` should be + cross-positive on ``K``:: + + 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 - The `L^{3}_{1}` cone is known to have a Lyapunov rank of one:: + """ + 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') - sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)]) - sage: lyapunov_rank(L31) - 1 + return all([ s*(L*x) >= 0 + for (x,s) in K.discrete_complementarity_set() ]) - Likewise for the `L^{3}_{\infty}` cone:: +def is_Z_on(L,K): + r""" + Determine whether or not ``L`` is a Z-operator on ``K``. - sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)]) - sage: lyapunov_rank(L3infty) - 1 + 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 Lyapunov rank should be additive on a product of cones:: + 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: 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) + 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. + + INPUT: + + - ``L`` -- A matrix over either an exact ring or ``SR``. + + - ``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 a Z-operator 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 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. + + EXAMPLES: + + 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 - 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}`:: + 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 - sage: K = Cone([(1,2,3), (-1,1,0), (1,0,6)]) - sage: lyapunov_rank(K) - 3 + TESTS: - The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K`` - itself:: + Everything in ``K.Z_operators_gens()`` should be a Z-operator + on ``K``:: - sage: K = Cone([(2,2,4), (-1,9,0), (2,0,6)]) - sage: lyapunov_rank(K) == lyapunov_rank(K.dual()) + 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 """ - V = K.lattice().vector_space() + return is_cross_positive_on(-L,K) + - xs = [V(x) for x in K.rays()] - ss = [V(s) for s in K.dual().rays()] +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. - # 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] + 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. - matrices = [x.column() * s.row() for (x,s) in C_of_K] + INPUT: - # 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) + - ``L`` -- A matrix over either an exact ring or ``SR``. - def phi(m): - r""" - Convert a matrix to a vector isomorphically. - """ - return W(m.list()) + - ``K`` -- A polyhedral closed convex cone. + + OUTPUT: - vectors = [phi(m) for m in matrices] + If the base ring of ``L`` is exact, then ``True`` will be returned if + and only if ``L`` is Lyapunov-like on ``K``. - return (W.dimension() - W.span(vectors).rank()) + 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 + + TESTS: + + The identity operator is always Lyapunov-like:: + + 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 + + The "zero" operator is always Lyapunov-like:: + + 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 + + Everything in ``K.lyapunov_like_basis()`` should be Lyapunov-like + on ``K``:: + + 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 + + """ + 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)