X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=mjo%2Fcone%2Fcone.py;h=4221566a52c1210fa1045b9df9305f04a45924a6;hb=3f2f1e2cc17d8796878a6971050f7aed7a37f737;hp=a7c16a520cbe19242398e176635ad3250e5c7671;hpb=9aea9a588808849ed2f01d0a2887fa9dffca9f36;p=sage.d.git diff --git a/mjo/cone/cone.py b/mjo/cone/cone.py index a7c16a5..4221566 100644 --- a/mjo/cone/cone.py +++ b/mjo/cone/cone.py @@ -1,122 +1,206 @@ from sage.all import * +from sage.geometry.cone import is_Cone def is_positive_on(L,K): r""" 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``. + 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: - - ``L`` -- A linear transformation or matrix. + - ``L`` -- A matrix over either an exact ring or ``SR``. - ``K`` -- A polyhedral closed convex cone. OUTPUT: - ``True`` if it can be proven that ``L`` is positive on ``K``, - and ``False`` otherwise. + 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. - .. WARNING:: + .. SEEALSO:: - 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. + :func:`is_cross_positive_on`, + :func:`is_Z_operator_on`, + :func:`is_lyapunov_like_on` EXAMPLES: - The identity is always positive in a nontrivial space:: + 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 + + Your matrix can be over any exact ring, but you may get unexpected + answers with weirder rings. For example, any rational matrix is + positive on the plane, but if your matrix contains polynomial + variables, the answer will be negative:: + + sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)]) + sage: K.is_full_space() + True + sage: L = matrix(QQ[x], [[x,0],[0,1]]) + sage: is_positive_on(L,K) + False + + The previous example is "unexpected" because it depends on how we + check whether or not ``L`` is positive. For exact base rings, we + check whether or not ``L*z`` belongs to ``K`` for each ``z in K``. + If ``K`` is closed, then an equally-valid test would be to check + whether the inner product of ``L*z`` and ``s`` is nonnegative for + every ``z`` in ``K`` and ``s`` in ``K.dual()``. In fact, that is + what we do over inexact rings. In the previous example, that test + would return an affirmative answer:: + + sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)]) + sage: L = matrix(QQ[x], [[x,0],[0,1]]) + sage: all([ (L*z).inner_product(s) for z in K for s in K.dual() ]) + True + sage: is_positive_on(L.change_ring(SR), K) + True + + TESTS: + + The identity operator is always positive:: sage: set_random_seed() - sage: K = random_cone(min_ambient_dim=1, max_ambient_dim=8) + sage: K = random_cone(max_ambient_dim=8) sage: L = identity_matrix(K.lattice_dim()) sage: is_positive_on(L,K) True - As is the "zero" transformation:: + The "zero" operator is always positive:: - sage: K = random_cone(min_ambient_dim=1, max_ambient_dim=8) + 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 - TESTS: - Everything in ``K.positive_operators_gens()`` should be positive on ``K``:: - 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() ]) + 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) - ....: for L in K.positive_operators_gens() ]) + sage: all([ is_positive_on(L.change_ring(SR),K) # long time + ....: for L in K.positive_operators_gens() ]) # long time True + Technically we could test this, but for now only closed convex cones + are supported as our ``K`` argument:: + + sage: K = [ vector([1,2,3]), vector([5,-1,7]) ] + sage: L = identity_matrix(3) + sage: is_positive_on(L,K) + Traceback (most recent call last): + ... + TypeError: K must be a Cone. + + We can't give reliable answers over inexact rings:: + + sage: K = Cone([(1,2,3), (4,5,6)]) + sage: L = identity_matrix(RR,3) + sage: is_positive_on(L,K) + Traceback (most recent call last): + ... + ValueError: The base ring of L is neither SR nor exact. + """ + + 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('The base ring of L is neither SR nor exact.') + if L.base_ring().is_exact(): - # This could potentially be extended to other types of ``K``... + # 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 ]) - elif L.base_ring() is SR: + 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 ]) - 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') + 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 ``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 + 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: - - ``L`` -- A linear transformation or matrix. + - ``L`` -- A matrix over either an exact ring or ``SR``. - ``K`` -- A polyhedral closed convex cone. OUTPUT: - ``True`` if it can be proven that ``L`` is cross-positive on ``K``, - and ``False`` otherwise. + If the base ring of ``L`` is exact, then ``True`` will be returned if + and only if ``L`` is cross-positive on ``K``. + + If the base ring of ``L`` is ``SR``, then the situation is more + complicated: - .. WARNING:: + - ``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. - 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. + .. SEEALSO:: + + :func:`is_positive_on`, + :func:`is_Z_operator_on`, + :func:`is_lyapunov_like_on` EXAMPLES: - The identity is always cross-positive in a nontrivial space:: + The identity operator is always cross-positive:: sage: set_random_seed() - sage: K = random_cone(min_ambient_dim=1, max_ambient_dim=8) + sage: K = random_cone(max_ambient_dim=8) sage: L = identity_matrix(K.lattice_dim()) sage: is_cross_positive_on(L,K) True - As is the "zero" transformation:: + The "zero" operator is always cross-positive:: - sage: K = random_cone(min_ambient_dim=1, max_ambient_dim=8) + 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) @@ -127,74 +211,103 @@ def is_cross_positive_on(L,K): 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() ]) + 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) - ....: for L in K.cross_positive_operators_gens() ]) + sage: all([ is_cross_positive_on(L.change_ring(SR),K) # long time + ....: for L in K.cross_positive_operators_gens() ]) # long time True + Technically we could test this, but for now only closed convex cones + are supported as our ``K`` argument:: + + sage: L = identity_matrix(3) + sage: K = [ vector([8,2,-8]), vector([5,-5,7]) ] + sage: is_cross_positive_on(L,K) + Traceback (most recent call last): + ... + TypeError: K must be a Cone. + + We can't give reliable answers over inexact rings:: + + sage: K = Cone([(1,2,3), (4,5,6)]) + sage: L = identity_matrix(RR,3) + sage: is_cross_positive_on(L,K) + Traceback (most recent call last): + ... + ValueError: The base ring of L is neither SR nor exact. + """ - 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') + 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('The base ring of L is neither SR nor exact.') + return all([ s*(L*x) >= 0 + for (x,s) in K.discrete_complementarity_set() ]) -def is_Z_on(L,K): +def is_Z_operator_on(L,K): r""" Determine whether or not ``L`` is a Z-operator on ``K``. - 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. + 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. A matrix is a Z-operator on ``K`` if and only if its negation is a cross-positive operator on ``K``. + 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 linear transformation or matrix. + - ``L`` -- A matrix over either an exact ring or ``SR``. - ``K`` -- A polyhedral closed convex cone. OUTPUT: - ``True`` if it can be proven that ``L`` is a Z-operator on ``K``, - and ``False`` otherwise. + 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. - .. WARNING:: + .. SEEALSO:: - 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. + :func:`is_positive_on`, + :func:`is_cross_positive_on`, + :func:`is_lyapunov_like_on` EXAMPLES: - The identity is always a Z-operator in a nontrivial space:: + The identity operator is always a Z-operator:: sage: set_random_seed() - sage: K = random_cone(min_ambient_dim=1, max_ambient_dim=8) + sage: K = random_cone(max_ambient_dim=8) sage: L = identity_matrix(K.lattice_dim()) - sage: is_Z_on(L,K) + sage: is_Z_operator_on(L,K) True - As is the "zero" transformation:: + The "zero" operator is always a Z-operator:: - sage: K = random_cone(min_ambient_dim=1, max_ambient_dim=8) + 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) + sage: is_Z_operator_on(L,K) True TESTS: @@ -202,14 +315,34 @@ def is_Z_on(L,K): Everything in ``K.Z_operators_gens()`` should be a Z-operator on ``K``:: - 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() ]) + sage: K = random_cone(max_ambient_dim=5) + sage: all([ is_Z_operator_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) - ....: for L in K.Z_operators_gens() ]) + sage: all([ is_Z_operator_on(L.change_ring(SR),K) # long time + ....: for L in K.Z_operators_gens() ]) # long time True + Technically we could test this, but for now only closed convex cones + are supported as our ``K`` argument:: + + sage: L = identity_matrix(3) + sage: K = [ vector([-4,20,3]), vector([1,-5,2]) ] + sage: is_Z_operator_on(L,K) + Traceback (most recent call last): + ... + TypeError: K must be a Cone. + + + We can't give reliable answers over inexact rings:: + + sage: K = Cone([(1,2,3), (4,5,6)]) + sage: L = identity_matrix(RR,3) + sage: is_Z_operator_on(L,K) + Traceback (most recent call last): + ... + ValueError: The base ring of L is neither SR nor exact. + """ return is_cross_positive_on(-L,K) @@ -218,74 +351,130 @@ 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 ``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 + 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. + An operator is Lyapunov-like on ``K`` if and only if both the + operator itself and its negation are cross-positive on ``K``. + + 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: - - ``L`` -- A linear transformation or matrix. + - ``L`` -- A matrix over either an exact ring or ``SR``. - ``K`` -- A polyhedral closed convex cone. OUTPUT: - ``True`` if it can be proven that ``L`` is Lyapunov-like on ``K``, - and ``False`` otherwise. + 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: - .. WARNING:: + - ``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. - 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. + .. SEEALSO:: + + :func:`is_positive_on`, + :func:`is_cross_positive_on`, + :func:`is_Z_operator_on` EXAMPLES: - The identity is always Lyapunov-like in a nontrivial space:: + 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(min_ambient_dim=1, max_ambient_dim=8) + sage: K = random_cone(max_ambient_dim=8) sage: L = identity_matrix(K.lattice_dim()) sage: is_lyapunov_like_on(L,K) True - As is the "zero" transformation:: + The "zero" operator is always Lyapunov-like:: - sage: K = random_cone(min_ambient_dim=1, max_ambient_dim=8) + 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 - TESTS: - Everything in ``K.lyapunov_like_basis()`` should be Lyapunov-like on ``K``:: - 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() ]) + 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) - ....: for L in K.lyapunov_like_basis() ]) + sage: all([ is_lyapunov_like_on(L.change_ring(SR),K) # long time + ....: for L in K.lyapunov_like_basis() ]) # long time + True + + Technically we could test this, but for now only closed convex cones + are supported as our ``K`` argument:: + + sage: L = identity_matrix(3) + sage: K = [ vector([2,2,-1]), vector([5,4,-3]) ] + sage: is_lyapunov_like_on(L,K) + Traceback (most recent call last): + ... + TypeError: K must be a Cone. + + We can't give reliable answers over inexact rings:: + + sage: K = Cone([(1,2,3), (4,5,6)]) + sage: L = identity_matrix(RR,3) + sage: is_lyapunov_like_on(L,K) + Traceback (most recent call last): + ... + ValueError: The base ring of L is neither SR nor exact. + + An operator is Lyapunov-like on a cone if and only if both the + operator and its negation are cross-positive on the cone:: + + sage: K = random_cone(max_ambient_dim=5) + sage: R = K.lattice().vector_space().base_ring() + sage: L = random_matrix(R, K.lattice_dim()) + sage: actual = is_lyapunov_like_on(L,K) # long time + sage: expected = (is_cross_positive_on(L,K) and # long time + ....: is_cross_positive_on(-L,K)) # long time + sage: actual == expected # long time True """ - 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') + 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('The base ring of L is neither SR nor exact.') + + # Even though ``discrete_complementarity_set`` is a cached method + # of cones, this is faster than calling ``is_cross_positive_on`` + # twice: doing so checks twice as many inequalities as the number + # of equalities that we're about to check. + return all([ s*(L*x) == 0 + for (x,s) in K.discrete_complementarity_set() ]) + def LL_cone(K): gens = K.lyapunov_like_basis()