X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=mjo%2Fcone%2Fcone.py;h=cbbfe9b692f3bdf36f68f202deedd467c684cef0;hb=115ecffcebcd0b7f86358aa38ddb52f1115e966c;hp=2e3dc8afb78cd1f4e1d5f368eb5f75733d127e82;hpb=433b47ffeac9305beb9101afa302c7e924b60744;p=sage.d.git diff --git a/mjo/cone/cone.py b/mjo/cone/cone.py index 2e3dc8a..cbbfe9b 100644 --- a/mjo/cone/cone.py +++ b/mjo/cone/cone.py @@ -1,263 +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 discrete_complementarity_set(K): +def is_positive_on(L,K): r""" - Compute the discrete complementarity set of this 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``. + + 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: - The complementarity set of this cone is the set of all orthogonal - pairs `(x,s)` such that `x` is in this cone, and `s` is in its - dual. The discrete complementarity set restricts `x` and `s` to be - generators of their respective cones. + - ``L`` -- A matrix over either an exact ring or ``SR``. + + - ``K`` -- A polyhedral closed convex cone. OUTPUT: - A list of pairs `(x,s)` such that, + 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: - * `x` is in this cone. - * `x` is a generator of this cone. - * `s` is in this cone's dual. - * `s` is a generator of this cone's dual. - * `x` and `s` are orthogonal. + - ``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: - The discrete complementarity set of the nonnegative orthant consists - of pairs of standard basis vectors:: + Nonnegative matrices are positive operators on the nonnegative + orthant:: - sage: K = Cone([(1,0),(0,1)]) - sage: discrete_complementarity_set(K) - [((1, 0), (0, 1)), ((0, 1), (1, 0))] + 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 - If the cone consists of a single ray, the second components of the - discrete complementarity set should generate the orthogonal - complement of that ray:: + TESTS: - sage: K = Cone([(1,0)]) - sage: discrete_complementarity_set(K) - [((1, 0), (0, 1)), ((1, 0), (0, -1))] - sage: K = Cone([(1,0,0)]) - sage: discrete_complementarity_set(K) - [((1, 0, 0), (0, 1, 0)), - ((1, 0, 0), (0, -1, 0)), - ((1, 0, 0), (0, 0, 1)), - ((1, 0, 0), (0, 0, -1))] + The identity operator is always positive:: - When the cone is the entire space, its dual is the trivial cone, so - the discrete complementarity set is empty:: + 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 - sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)]) - sage: discrete_complementarity_set(K) - [] + The "zero" operator is always positive:: - TESTS: + 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 complementarity set of the dual can be obtained by switching the - components of the complementarity set of the original cone:: + Everything in ``K.positive_operators_gens()`` should be + positive on ``K``:: - sage: K1 = random_cone(max_dim=10, max_rays=10) - sage: K2 = K1.dual() - sage: expected = [(x,s) for (s,x) in discrete_complementarity_set(K2)] - sage: actual = discrete_complementarity_set(K1) - sage: actual == expected + 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 """ - V = K.lattice().vector_space() + 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``. - # Convert the rays to vectors so that we can compute inner - # products. - xs = [V(x) for x in K.rays()] - ss = [V(s) for s in K.dual().rays()] + 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. - return [(x,s) for x in xs for s in ss if x.inner_product(s) == 0] + 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: -def LL(K): - r""" - Compute the space `\mathbf{LL}` of all Lyapunov-like transformations - on this cone. + - ``L`` -- A matrix over either an exact ring or ``SR``. + + - ``K`` -- A polyhedral closed convex cone. OUTPUT: - A ``MatrixSpace`` object `M` such that every matrix `L \in M` is - Lyapunov-like on this cone. + If the base ring of ``L`` is exact, then ``True`` will be returned if + and only if ``L`` is cross-positive on ``K``. - """ - V = K.lattice().vector_space() + If the base ring of ``L`` is ``SR``, then the situation is more + complicated: - C_of_K = discrete_complementarity_set(K) + - ``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. - tensor_products = [s.tensor_product(x) for (x,s) in C_of_K] + EXAMPLES: - # 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) + The identity operator is always cross-positive:: - # Turn our matrices into long vectors... - vectors = [ W(m.list()) for m in tensor_products ] + 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 - # Vector space representation of Lyapunov-like matrices - # (i.e. vec(L) where L is Luapunov-like). - LL_vector = W.span(vectors).complement() + The "zero" operator is always cross-positive:: - # Now construct an ambient MatrixSpace in which to stick our - # transformations. - M = MatrixSpace(V.base_ring(), V.dimension()) + 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 - matrix_basis = [ M(v.list()) for v in LL_vector.basis() ] + TESTS: - return matrix_basis + 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 + + """ + 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 lyapunov_rank(K): +def is_Z_on(L,K): r""" - Compute the Lyapunov (or bilinearity) rank of this cone. + Determine whether or not ``L`` is a Z-operator on ``K``. - The Lyapunov rank of a cone can be thought of in (mainly) two ways: + 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. - 1. The dimension of the Lie algebra of the automorphism group of the - cone. + A matrix is a Z-operator on ``K`` if and only if its negation is a + cross-positive operator on ``K``. - 2. The dimension of the linear space of all Lyapunov-like - transformations on the cone. + 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: - 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 a Z-operator 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 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. - .. seealso:: + EXAMPLES: - :meth:`is_proper` + The identity operator is always a Z-operator:: - ALGORITHM: + 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 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 a Z-operator:: - 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_Z_on(L,K) + True - .. [Gowda/Tao] 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: - .. [Rudolf et al.] 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.Z_operators_gens()`` should be a Z-operator + on ``K``:: - EXAMPLES: + 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 - The nonnegative orthant in `\mathbb{R}^{n}` always has rank `n` - [Rudolf et al.]_:: + """ + return is_cross_positive_on(-L,K) - 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 - The `L^{3}_{1}` cone is known to have a Lyapunov rank of one - [Rudolf et al.]_:: +def is_lyapunov_like_on(L,K): + r""" + Determine whether or not ``L`` is Lyapunov-like on ``K``. - sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)]) - sage: lyapunov_rank(L31) - 1 + 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. - Likewise for the `L^{3}_{\infty}` cone [Rudolf et al.]_:: + 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. - sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)]) - sage: lyapunov_rank(L3infty) - 1 + INPUT: - The Lyapunov rank should be additive on a product of cones - [Rudolf et al.]_:: + - ``L`` -- A matrix over either an exact ring or ``SR``. - 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 + - ``K`` -- A polyhedral closed convex cone. - Two isomorphic cones should have the same Lyapunov rank [Rudolf et al.]_. - The cone ``K`` in the following example is isomorphic to the nonnegative - octant in `\mathbb{R}^{3}`:: + OUTPUT: - sage: K = Cone([(1,2,3), (-1,1,0), (1,0,6)]) - sage: lyapunov_rank(K) - 3 + If the base ring of ``L`` is exact, then ``True`` will be returned if + and only if ``L`` is Lyapunov-like on ``K``. - The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K`` - itself [Rudolf et al.]_:: + If the base ring of ``L`` is ``SR``, then the situation is more + complicated: - sage: K = Cone([(2,2,4), (-1,9,0), (2,0,6)]) - sage: lyapunov_rank(K) == lyapunov_rank(K.dual()) + - ``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 Lyapunov rank should be additive on a product of cones - [Rudolf et al.]_:: + The identity operator is always Lyapunov-like:: - sage: K1 = random_cone(max_dim=10, max_rays=10) - sage: K2 = random_cone(max_dim=10, max_rays=10) - sage: K = K1.cartesian_product(K2) - sage: lyapunov_rank(K) == lyapunov_rank(K1) + lyapunov_rank(K2) + 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 dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K`` - itself [Rudolf et al.]_:: + The "zero" operator is always Lyapunov-like:: - sage: K = random_cone(max_dim=10, max_rays=10) - sage: lyapunov_rank(K) == lyapunov_rank(K.dual()) + 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 - The Lyapunov rank of a proper polyhedral cone in `n` dimensions can - be any number between `1` and `n` inclusive, excluding `n-1` - [Gowda/Tao]_. By accident, the `n-1` restriction will hold for the - trivial cone in a trivial space as well. However, in zero dimensions, - the Lyapunov rank of the trivial cone will be zero:: + Everything in ``K.lyapunov_like_basis()`` should be Lyapunov-like + on ``K``:: - sage: K = random_cone(max_dim=10, strictly_convex=True, solid=True) - sage: b = lyapunov_rank(K) - sage: n = K.lattice_dim() - sage: (n == 0 or 1 <= b) and b <= n + 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 - sage: b == n-1 - False """ - return len(LL(K)) + 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)