X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=mjo%2Fcone%2Fcone.py;h=01ef5f9d5cc56b087344237b9b5c915423ac0921;hb=344df04250bdc5982eab521276fc08b509fd239f;hp=2e3dc8afb78cd1f4e1d5f368eb5f75733d127e82;hpb=433b47ffeac9305beb9101afa302c7e924b60744;p=sage.d.git diff --git a/mjo/cone/cone.py b/mjo/cone/cone.py index 2e3dc8a..01ef5f9 100644 --- a/mjo/cone/cone.py +++ b/mjo/cone/cone.py @@ -1,263 +1,326 @@ -# 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 * - -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 ``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``. + + INPUT: + + - ``L`` -- A linear transformation or matrix. - 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. + - ``K`` -- A polyhedral closed convex cone. OUTPUT: - A list of pairs `(x,s)` such that, + ``True`` if it can be proven that ``L`` is positive on ``K``, + and ``False`` otherwise. - * `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. + .. WARNING:: + + If this function returns ``True``, then ``L`` is positive + on ``K``. However, if ``False`` is returned, that could mean one + of two things. The first is that ``L`` is definitely not + positive on ``K``. The second is more of an "I don't know" + answer, returned (for example) if we cannot prove that an inner + product is nonnegative. EXAMPLES: - The discrete complementarity set of the nonnegative orthant consists - of pairs of standard basis vectors:: + Positive operators on the nonnegative orthant are nonnegative + matrices:: - 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 is always positive in a nontrivial space:: - 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(min_ambient_dim=1, 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) - [] + As is the "zero" transformation:: - TESTS: + sage: K = random_cone(min_ambient_dim=1, 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(min_ambient_dim=1, max_ambient_dim=6) + sage: all([ is_positive_on(L,K) + ....: for L in K.positive_operators_gens() ]) + True + sage: all([ is_positive_on(L.change_ring(SR),K) + ....: for L in K.positive_operators_gens() ]) True """ - V = K.lattice().vector_space() + if L.base_ring().is_exact(): + # This could potentially be extended to other types of ``K``... + return all([ L*x in K for x in K ]) + elif L.base_ring() is SR: + # 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() ]) + else: + # The only inexact ring that we're willing to work with is SR, + # since it can still be exact when working with symbolic + # constants like pi and e. + raise ValueError('base ring of operator L is neither SR nor exact') + + +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 ``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] + INPUT: + - ``L`` -- A linear transformation or matrix. -def LL(K): - r""" - Compute the space `\mathbf{LL}` of all Lyapunov-like transformations - on this cone. + - ``K`` -- A polyhedral closed convex cone. OUTPUT: - A ``MatrixSpace`` object `M` such that every matrix `L \in M` is - Lyapunov-like on this cone. + ``True`` if it can be proven that ``L`` is cross-positive on ``K``, + and ``False`` otherwise. - """ - V = K.lattice().vector_space() + .. WARNING:: - C_of_K = discrete_complementarity_set(K) + If this function returns ``True``, then ``L`` is cross-positive + on ``K``. However, if ``False`` is returned, that could mean one + of two things. The first is that ``L`` is definitely not + cross-positive on ``K``. The second is more of an "I don't know" + answer, returned (for example) if we cannot prove that an inner + product is nonnegative. - 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 is always cross-positive in a nontrivial space:: - # Turn our matrices into long vectors... - vectors = [ W(m.list()) for m in tensor_products ] + sage: set_random_seed() + sage: K = random_cone(min_ambient_dim=1, 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() + As is the "zero" transformation:: - # Now construct an ambient MatrixSpace in which to stick our - # transformations. - M = MatrixSpace(V.base_ring(), V.dimension()) + sage: K = random_cone(min_ambient_dim=1, 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(min_ambient_dim=1, max_ambient_dim=6) + sage: all([ is_cross_positive_on(L,K) + ....: for L in K.cross_positive_operators_gens() ]) + True + sage: all([ is_cross_positive_on(L.change_ring(SR),K) + ....: for L in K.cross_positive_operators_gens() ]) + True + """ + if L.base_ring().is_exact() or L.base_ring() is SR: + return all([ s*(L*x) >= 0 + for (x,s) in K.discrete_complementarity_set() ]) + else: + # The only inexact ring that we're willing to work with is SR, + # since it can still be exact when working with symbolic + # constants like pi and e. + raise ValueError('base ring of operator L is neither SR nor exact') -def lyapunov_rank(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: +def is_Z_on(L,K): + r""" + Determine whether or not ``L`` is a Z-operator on ``K``. - 1. The dimension of the Lie algebra of the automorphism group of the - cone. + We say that ``L`` is a Z-operator on ``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. - 2. The dimension of the linear space of all Lyapunov-like - transformations on the cone. + A matrix is a Z-operator on ``K`` if and only if its negation is a + cross-positive operator on ``K``. INPUT: - A closed, convex polyhedral cone. + - ``L`` -- A linear transformation or matrix. + + - ``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). + ``True`` if it can be proven that ``L`` is a Z-operator on ``K``, + and ``False`` otherwise. - .. note:: + .. WARNING:: - In the references, the cones are always assumed to be proper. We - do not impose this restriction. + If this function returns ``True``, then ``L`` is a Z-operator + on ``K``. However, if ``False`` is returned, that could mean one + of two things. The first is that ``L`` is definitely not + a Z-operator on ``K``. The second is more of an "I don't know" + answer, returned (for example) if we cannot prove that an inner + product is nonnegative. - .. seealso:: + EXAMPLES: - :meth:`is_proper` + The identity is always a Z-operator in a nontrivial space:: - ALGORITHM: + sage: set_random_seed() + sage: K = random_cone(min_ambient_dim=1, 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}`. + As is the "zero" transformation:: - REFERENCES: + sage: K = random_cone(min_ambient_dim=1, 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(min_ambient_dim=1, max_ambient_dim=6) + sage: all([ is_Z_on(L,K) + ....: for L in K.Z_operators_gens() ]) + True + sage: all([ is_Z_on(L.change_ring(SR),K) + ....: for L in K.Z_operators_gens() ]) + 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 ``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.]_:: + INPUT: - sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)]) - sage: lyapunov_rank(L3infty) - 1 + - ``L`` -- A linear transformation or matrix. - The Lyapunov rank should be additive on a product of cones - [Rudolf et al.]_:: + - ``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: + + ``True`` if it can be proven that ``L`` is Lyapunov-like on ``K``, + and ``False`` otherwise. - 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}`:: + .. WARNING:: - sage: K = Cone([(1,2,3), (-1,1,0), (1,0,6)]) - sage: lyapunov_rank(K) - 3 + If this function returns ``True``, then ``L`` is Lyapunov-like + on ``K``. However, if ``False`` is returned, that could mean one + of two things. The first is that ``L`` is definitely not + Lyapunov-like on ``K``. The second is more of an "I don't know" + answer, returned (for example) if we cannot prove that an inner + product is zero. - The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K`` - itself [Rudolf et al.]_:: + EXAMPLES: - sage: K = Cone([(2,2,4), (-1,9,0), (2,0,6)]) - sage: lyapunov_rank(K) == lyapunov_rank(K.dual()) + Lyapunov-like operators on the nonnegative orthant are diagonal + matrices:: + + 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 is always Lyapunov-like in a nontrivial space:: - 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(min_ambient_dim=1, 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.]_:: + As is the "zero" transformation:: - sage: K = random_cone(max_dim=10, max_rays=10) - sage: lyapunov_rank(K) == lyapunov_rank(K.dual()) + sage: K = random_cone(min_ambient_dim=1, 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(min_ambient_dim=1, max_ambient_dim=6) + sage: all([ is_lyapunov_like_on(L,K) + ....: for L in K.lyapunov_like_basis() ]) + True + sage: all([ is_lyapunov_like_on(L.change_ring(SR),K) + ....: for L in K.lyapunov_like_basis() ]) True - sage: b == n-1 - False """ - return len(LL(K)) + if L.base_ring().is_exact() or L.base_ring() is SR: + # The "fast method" of creating a vector space based on a + # ``lyapunov_like_basis`` is actually slower than this. + return all([ s*(L*x) == 0 + for (x,s) in K.discrete_complementarity_set() ]) + else: + # The only inexact ring that we're willing to work with is SR, + # since it can still be exact when working with symbolic + # constants like pi and e. + raise ValueError('base ring of operator L is neither SR nor exact') + +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)