mjo/cone/cone.py

   1 from sage.all import *
   2
   3 def is_cross_positive(L,K):
   4     r"""
   5     Determine whether or not ``L`` is cross-positive on ``K``.
   6
   7     We say that ``L`` is cross-positive on ``K`` if `\left\langle
   8     L\left\lparenx\right\rparen,s\right\rangle >= 0` for all pairs
   9     `\left\langle x,s \right\rangle` in the complementarity set of
  10     ``K``. It is known that this property need only be
  11     checked for generators of ``K`` and its dual.
  12
  13     INPUT:
  14
  15     - ``L`` -- A linear transformation or matrix.
  16
  17     - ``K`` -- A polyhedral closed convex cone.
  18
  19     OUTPUT:
  20
  21     ``True`` if it can be proven that ``L`` is cross-positive on ``K``,
  22     and ``False`` otherwise.
  23
  24     .. WARNING::
  25
  26         If this function returns ``True``, then ``L`` is cross-positive
  27         on ``K``. However, if ``False`` is returned, that could mean one
  28         of two things. The first is that ``L`` is definitely not
  29         cross-positive on ``K``. The second is more of an "I don't know"
  30         answer, returned (for example) if we cannot prove that an inner
  31         product is nonnegative.
  32
  33     EXAMPLES:
  34
  35     The identity is always cross-positive in a nontrivial space::
  36
  37         sage: set_random_seed()
  38         sage: K = random_cone(min_ambient_dim=1, max_ambient_dim=8)
  39         sage: L = identity_matrix(K.lattice_dim())
  40         sage: is_cross_positive(L,K)
  41         True
  42
  43     As is the "zero" transformation::
  44
  45         sage: K = random_cone(min_ambient_dim=1, max_ambient_dim=8)
  46         sage: R = K.lattice().vector_space().base_ring()
  47         sage: L = zero_matrix(R, K.lattice_dim())
  48         sage: is_cross_positive(L,K)
  49         True
  50
  51         Everything in ``K.cross_positive_operator_gens()`` should be
  52         cross-positive on ``K``::
  53
  54         sage: K = random_cone(min_ambient_dim=1, max_ambient_dim=6)
  55         sage: all([ is_cross_positive(L,K)
  56         ....:       for L in K.cross_positive_operator_gens() ])
  57         True
  58
  59     """
  60     if L.base_ring().is_exact() or L.base_ring() is SR:
  61         return all([ s*(L*x) >= 0
  62                      for (x,s) in K.discrete_complementarity_set() ])
  63     else:
  64         # The only inexact ring that we're willing to work with is SR,
  65         # since it can still be exact when working with symbolic
  66         # constants like pi and e.
  67         raise ValueError('base ring of operator L is neither SR nor exact')
  68
  69
  70 def is_lyapunov_like(L,K):
  71     r"""
  72     Determine whether or not ``L`` is Lyapunov-like on ``K``.
  73
  74     We say that ``L`` is Lyapunov-like on ``K`` if `\left\langle
  75     L\left\lparenx\right\rparen,s\right\rangle = 0` for all pairs
  76     `\left\langle x,s \right\rangle` in the complementarity set of
  77     ``K``. It is known [Orlitzky]_ that this property need only be
  78     checked for generators of ``K`` and its dual.
  79
  80     There are faster ways of checking this property. For example, we
  81     could compute a `lyapunov_like_basis` of the cone, and then test
  82     whether or not the given matrix is contained in the span of that
  83     basis. The value of this function is that it works on symbolic
  84     matrices.
  85
  86     INPUT:
  87
  88     - ``L`` -- A linear transformation or matrix.
  89
  90     - ``K`` -- A polyhedral closed convex cone.
  91
  92     OUTPUT:
  93
  94     ``True`` if it can be proven that ``L`` is Lyapunov-like on ``K``,
  95     and ``False`` otherwise.
  96
  97     .. WARNING::
  98
  99         If this function returns ``True``, then ``L`` is Lyapunov-like
 100         on ``K``. However, if ``False`` is returned, that could mean one
 101         of two things. The first is that ``L`` is definitely not
 102         Lyapunov-like on ``K``. The second is more of an "I don't know"
 103         answer, returned (for example) if we cannot prove that an inner
 104         product is zero.
 105
 106     REFERENCES:
 107
 108     M. Orlitzky. The Lyapunov rank of an improper cone.
 109     http://www.optimization-online.org/DB_HTML/2015/10/5135.html
 110
 111     EXAMPLES:
 112
 113     The identity is always Lyapunov-like in a nontrivial space::
 114
 115         sage: set_random_seed()
 116         sage: K = random_cone(min_ambient_dim=1, max_ambient_dim=8)
 117         sage: L = identity_matrix(K.lattice_dim())
 118         sage: is_lyapunov_like(L,K)
 119         True
 120
 121     As is the "zero" transformation::
 122
 123         sage: K = random_cone(min_ambient_dim=1, max_ambient_dim=8)
 124         sage: R = K.lattice().vector_space().base_ring()
 125         sage: L = zero_matrix(R, K.lattice_dim())
 126         sage: is_lyapunov_like(L,K)
 127         True
 128
 129         Everything in ``K.lyapunov_like_basis()`` should be Lyapunov-like
 130         on ``K``::
 131
 132         sage: K = random_cone(min_ambient_dim=1, max_ambient_dim=6)
 133         sage: all([ is_lyapunov_like(L,K) for L in K.lyapunov_like_basis() ])
 134         True
 135
 136     """
 137     if L.base_ring().is_exact() or L.base_ring() is SR:
 138         V = VectorSpace(K.lattice().base_field(), K.lattice_dim()**2)
 139         LL_of_K = V.span([ V(m.list()) for m in K.lyapunov_like_basis() ])
 140         return V(L.list()) in LL_of_K
 141     else:
 142         # The only inexact ring that we're willing to work with is SR,
 143         # since it can still be exact when working with symbolic
 144         # constants like pi and e.
 145         raise ValueError('base ring of operator L is neither SR nor exact')
 146
 147 def LL_cone(K):
 148     gens = K.lyapunov_like_basis()
 149     L = ToricLattice(K.lattice_dim()**2)
 150     return Cone([ g.list() for g in gens ], lattice=L, check=False)
 151
 152 def Sigma_cone(K):
 153     gens = K.cross_positive_operator_gens()
 154     L = ToricLattice(K.lattice_dim()**2)
 155     return Cone([ g.list() for g in gens ], lattice=L, check=False)
 156
 157 def Z_cone(K):
 158     gens = K.Z_operator_gens()
 159     L = ToricLattice(K.lattice_dim()**2)
 160     return Cone([ g.list() for g in gens ], lattice=L, check=False)
 161
 162 def pi_cone(K1, K2=None):
 163     if K2 is None:
 164         K2 = K1
 165     gens = K1.positive_operator_gens(K2)
 166     L = ToricLattice(K1.lattice_dim()*K2.lattice_dim())
 167     return Cone([ g.list() for g in gens ], lattice=L, check=False)