mjo/cone/cone.py

   1 from sage.all import *
   2 from sage.geometry.cone import is_Cone
   3
   4 def is_positive_on(L,K):
   5     r"""
   6     Determine whether or not ``L`` is positive on ``K``.
   7
   8     We say that ``L`` is positive on a closed convex cone ``K`` if
   9     `L\left\lparen x \right\rparen` belongs to ``K`` for all `x` in
  10     ``K``. This property need only be checked for generators of ``K``.
  11
  12     To reliably check whether or not ``L`` is positive, its base ring
  13     must be either exact (for example, the rationals) or ``SR``. An
  14     exact ring is more reliable, but in some cases a matrix whose
  15     entries contain symbolic constants like ``e`` and ``pi`` will work.
  16
  17     INPUT:
  18
  19     - ``L`` -- A matrix over either an exact ring or ``SR``.
  20
  21     - ``K`` -- A polyhedral closed convex cone.
  22
  23     OUTPUT:
  24
  25     If the base ring of ``L`` is exact, then ``True`` will be returned if
  26     and only if ``L`` is positive on ``K``.
  27
  28     If the base ring of ``L`` is ``SR``, then the situation is more
  29     complicated:
  30
  31     - ``True`` will be returned if it can be proven that ``L``
  32       is positive on ``K``.
  33     - ``False`` will be returned if it can be proven that ``L``
  34       is not positive on ``K``.
  35     - ``False`` will also be returned if we can't decide; specifically
  36       if we arrive at a symbolic inequality that cannot be resolved.
  37
  38     .. SEEALSO::
  39
  40           :func:`is_cross_positive_on`,
  41           :func:`is_Z_operator_on`,
  42           :func:`is_lyapunov_like_on`
  43
  44     EXAMPLES:
  45
  46     Nonnegative matrices are positive operators on the nonnegative
  47     orthant::
  48
  49         sage: K = Cone([(1,0,0),(0,1,0),(0,0,1)])
  50         sage: L = random_matrix(QQ,3).apply_map(abs)
  51         sage: is_positive_on(L,K)
  52         True
  53
  54     Your matrix can be over any exact ring, but you may get unexpected
  55     answers with weirder rings. For example, any rational matrix is
  56     positive on the plane, but if your matrix contains polynomial
  57     variables, the answer will be negative::
  58
  59         sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
  60         sage: K.is_full_space()
  61         True
  62         sage: L = matrix(QQ[x], [[x,0],[0,1]])
  63         sage: is_positive_on(L,K)
  64         False
  65
  66     The previous example is "unexpected" because it depends on how we
  67     check whether or not ``L`` is positive. For exact base rings, we
  68     check whether or not ``L*z`` belongs to ``K`` for each ``z in K``.
  69     If ``K`` is closed, then an equally-valid test would be to check
  70     whether the inner product of ``L*z`` and ``s`` is nonnegative for
  71     every ``z`` in ``K`` and ``s`` in ``K.dual()``. In fact, that is
  72     what we do over inexact rings. In the previous example, that test
  73     would return an affirmative answer::
  74
  75         sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
  76         sage: L = matrix(QQ[x], [[x,0],[0,1]])
  77         sage: all([ (L*z).inner_product(s) for z in K for s in K.dual() ])
  78         True
  79         sage: is_positive_on(L.change_ring(SR), K)
  80         True
  81
  82     TESTS:
  83
  84     The identity operator is always positive::
  85
  86         sage: set_random_seed()
  87         sage: K = random_cone(max_ambient_dim=8)
  88         sage: L = identity_matrix(K.lattice_dim())
  89         sage: is_positive_on(L,K)
  90         True
  91
  92     The "zero" operator is always positive::
  93
  94         sage: K = random_cone(max_ambient_dim=8)
  95         sage: R = K.lattice().vector_space().base_ring()
  96         sage: L = zero_matrix(R, K.lattice_dim())
  97         sage: is_positive_on(L,K)
  98         True
  99
 100     Everything in ``K.positive_operators_gens()`` should be
 101     positive on ``K``::
 102
 103         sage: K = random_cone(max_ambient_dim=5)
 104         sage: all([ is_positive_on(L,K)                     # long time
 105         ....:       for L in K.positive_operators_gens() ]) # long time
 106         True
 107         sage: all([ is_positive_on(L.change_ring(SR),K)     # long time
 108         ....:       for L in K.positive_operators_gens() ]) # long time
 109         True
 110
 111     Technically we could test this, but for now only closed convex cones
 112     are supported as our ``K`` argument::
 113
 114         sage: K = [ vector([1,2,3]), vector([5,-1,7]) ]
 115         sage: L = identity_matrix(3)
 116         sage: is_positive_on(L,K)
 117         Traceback (most recent call last):
 118         ...
 119         TypeError: K must be a Cone.
 120
 121     We can't give reliable answers over inexact rings::
 122
 123         sage: K = Cone([(1,2,3), (4,5,6)])
 124         sage: L = identity_matrix(RR,3)
 125         sage: is_positive_on(L,K)
 126         Traceback (most recent call last):
 127         ...
 128         ValueError: The base ring of L is neither SR nor exact.
 129
 130     """
 131
 132     if not is_Cone(K):
 133         raise TypeError('K must be a Cone.')
 134     if not L.base_ring().is_exact() and not L.base_ring() is SR:
 135         raise ValueError('The base ring of L is neither SR nor exact.')
 136
 137     if L.base_ring().is_exact():
 138         # This should be way faster than computing the dual and
 139         # checking a bunch of inequalities, but it doesn't work if
 140         # ``L*x`` is symbolic. For example, ``e in Cone([(1,)])``
 141         # is true, but returns ``False``.
 142         return all([ L*x in K for x in K ])
 143     else:
 144         # Fall back to inequality-checking when the entries of ``L``
 145         # might be symbolic.
 146         return all([ s*(L*x) >= 0 for x in K for s in K.dual() ])
 147
 148
 149 def is_cross_positive_on(L,K):
 150     r"""
 151     Determine whether or not ``L`` is cross-positive on ``K``.
 152
 153     We say that ``L`` is cross-positive on a closed convex cone``K`` if
 154     `\left\langle L\left\lparenx\right\rparen,s\right\rangle \ge 0` for
 155     all pairs `\left\langle x,s \right\rangle` in the complementarity
 156     set of ``K``. This property need only be checked for generators of
 157     ``K`` and its dual.
 158
 159     To reliably check whether or not ``L`` is cross-positive, its base
 160     ring must be either exact (for example, the rationals) or ``SR``. An
 161     exact ring is more reliable, but in some cases a matrix whose
 162     entries contain symbolic constants like ``e`` and ``pi`` will work.
 163
 164     INPUT:
 165
 166     - ``L`` -- A matrix over either an exact ring or ``SR``.
 167
 168     - ``K`` -- A polyhedral closed convex cone.
 169
 170     OUTPUT:
 171
 172     If the base ring of ``L`` is exact, then ``True`` will be returned if
 173     and only if ``L`` is cross-positive on ``K``.
 174
 175     If the base ring of ``L`` is ``SR``, then the situation is more
 176     complicated:
 177
 178     - ``True`` will be returned if it can be proven that ``L``
 179       is cross-positive on ``K``.
 180     - ``False`` will be returned if it can be proven that ``L``
 181       is not cross-positive on ``K``.
 182     - ``False`` will also be returned if we can't decide; specifically
 183       if we arrive at a symbolic inequality that cannot be resolved.
 184
 185     .. SEEALSO::
 186
 187           :func:`is_positive_on`,
 188           :func:`is_Z_operator_on`,
 189           :func:`is_lyapunov_like_on`
 190
 191     EXAMPLES:
 192
 193     The identity operator is always cross-positive::
 194
 195         sage: set_random_seed()
 196         sage: K = random_cone(max_ambient_dim=8)
 197         sage: L = identity_matrix(K.lattice_dim())
 198         sage: is_cross_positive_on(L,K)
 199         True
 200
 201     The "zero" operator is always cross-positive::
 202
 203         sage: K = random_cone(max_ambient_dim=8)
 204         sage: R = K.lattice().vector_space().base_ring()
 205         sage: L = zero_matrix(R, K.lattice_dim())
 206         sage: is_cross_positive_on(L,K)
 207         True
 208
 209     TESTS:
 210
 211     Everything in ``K.cross_positive_operators_gens()`` should be
 212     cross-positive on ``K``::
 213
 214         sage: K = random_cone(max_ambient_dim=5)
 215         sage: all([ is_cross_positive_on(L,K)                     # long time
 216         ....:       for L in K.cross_positive_operators_gens() ]) # long time
 217         True
 218         sage: all([ is_cross_positive_on(L.change_ring(SR),K)     # long time
 219         ....:       for L in K.cross_positive_operators_gens() ]) # long time
 220         True
 221
 222     Technically we could test this, but for now only closed convex cones
 223     are supported as our ``K`` argument::
 224
 225         sage: L = identity_matrix(3)
 226         sage: K = [ vector([8,2,-8]), vector([5,-5,7]) ]
 227         sage: is_cross_positive_on(L,K)
 228         Traceback (most recent call last):
 229         ...
 230         TypeError: K must be a Cone.
 231
 232     We can't give reliable answers over inexact rings::
 233
 234         sage: K = Cone([(1,2,3), (4,5,6)])
 235         sage: L = identity_matrix(RR,3)
 236         sage: is_cross_positive_on(L,K)
 237         Traceback (most recent call last):
 238         ...
 239         ValueError: The base ring of L is neither SR nor exact.
 240
 241     """
 242     if not is_Cone(K):
 243         raise TypeError('K must be a Cone.')
 244     if not L.base_ring().is_exact() and not L.base_ring() is SR:
 245         raise ValueError('The base ring of L is neither SR nor exact.')
 246
 247     return all([ s*(L*x) >= 0
 248                  for (x,s) in K.discrete_complementarity_set() ])
 249
 250 def is_Z_operator_on(L,K):
 251     r"""
 252     Determine whether or not ``L`` is a Z-operator on ``K``.
 253
 254     We say that ``L`` is a Z-operator on a closed convex cone``K`` if
 255     `\left\langle L\left\lparenx\right\rparen,s\right\rangle \le 0` for
 256     all pairs `\left\langle x,s \right\rangle` in the complementarity
 257     set of ``K``. It is known that this property need only be checked
 258     for generators of ``K`` and its dual.
 259
 260     A matrix is a Z-operator on ``K`` if and only if its negation is a
 261     cross-positive operator on ``K``.
 262
 263     To reliably check whether or not ``L`` is a Z operator, its base
 264     ring must be either exact (for example, the rationals) or ``SR``. An
 265     exact ring is more reliable, but in some cases a matrix whose
 266     entries contain symbolic constants like ``e`` and ``pi`` will work.
 267
 268     INPUT:
 269
 270     - ``L`` -- A matrix over either an exact ring or ``SR``.
 271
 272     - ``K`` -- A polyhedral closed convex cone.
 273
 274     OUTPUT:
 275
 276     If the base ring of ``L`` is exact, then ``True`` will be returned if
 277     and only if ``L`` is a Z-operator on ``K``.
 278
 279     If the base ring of ``L`` is ``SR``, then the situation is more
 280     complicated:
 281
 282     - ``True`` will be returned if it can be proven that ``L``
 283       is a Z-operator on ``K``.
 284     - ``False`` will be returned if it can be proven that ``L``
 285       is not a Z-operator on ``K``.
 286     - ``False`` will also be returned if we can't decide; specifically
 287       if we arrive at a symbolic inequality that cannot be resolved.
 288
 289     .. SEEALSO::
 290
 291           :func:`is_positive_on`,
 292           :func:`is_cross_positive_on`,
 293           :func:`is_lyapunov_like_on`
 294
 295     EXAMPLES:
 296
 297     The identity operator is always a Z-operator::
 298
 299         sage: set_random_seed()
 300         sage: K = random_cone(max_ambient_dim=8)
 301         sage: L = identity_matrix(K.lattice_dim())
 302         sage: is_Z_operator_on(L,K)
 303         True
 304
 305     The "zero" operator is always a Z-operator::
 306
 307         sage: K = random_cone(max_ambient_dim=8)
 308         sage: R = K.lattice().vector_space().base_ring()
 309         sage: L = zero_matrix(R, K.lattice_dim())
 310         sage: is_Z_operator_on(L,K)
 311         True
 312
 313     TESTS:
 314
 315     Everything in ``K.Z_operators_gens()`` should be a Z-operator
 316     on ``K``::
 317
 318         sage: K = random_cone(max_ambient_dim=5)
 319         sage: all([ is_Z_operator_on(L,K)            # long time
 320         ....:       for L in K.Z_operators_gens() ]) # long time
 321         True
 322         sage: all([ is_Z_operator_on(L.change_ring(SR),K) # long time
 323         ....:       for L in K.Z_operators_gens() ])      # long time
 324         True
 325
 326     Technically we could test this, but for now only closed convex cones
 327     are supported as our ``K`` argument::
 328
 329         sage: L = identity_matrix(3)
 330         sage: K = [ vector([-4,20,3]), vector([1,-5,2]) ]
 331         sage: is_Z_operator_on(L,K)
 332         Traceback (most recent call last):
 333         ...
 334         TypeError: K must be a Cone.
 335
 336
 337     We can't give reliable answers over inexact rings::
 338
 339         sage: K = Cone([(1,2,3), (4,5,6)])
 340         sage: L = identity_matrix(RR,3)
 341         sage: is_Z_operator_on(L,K)
 342         Traceback (most recent call last):
 343         ...
 344         ValueError: The base ring of L is neither SR nor exact.
 345
 346     """
 347     return is_cross_positive_on(-L,K)
 348
 349
 350 def is_lyapunov_like_on(L,K):
 351     r"""
 352     Determine whether or not ``L`` is Lyapunov-like on ``K``.
 353
 354     We say that ``L`` is Lyapunov-like on a closed convex cone ``K`` if
 355     `\left\langle L\left\lparenx\right\rparen,s\right\rangle = 0` for
 356     all pairs `\left\langle x,s \right\rangle` in the complementarity
 357     set of ``K``. This property need only be checked for generators of
 358     ``K`` and its dual.
 359
 360     An operator is Lyapunov-like on ``K`` if and only if both the
 361     operator itself and its negation are cross-positive on ``K``.
 362
 363     To reliably check whether or not ``L`` is Lyapunov-like, its base
 364     ring must be either exact (for example, the rationals) or ``SR``. An
 365     exact ring is more reliable, but in some cases a matrix whose
 366     entries contain symbolic constants like ``e`` and ``pi`` will work.
 367
 368     INPUT:
 369
 370     - ``L`` -- A matrix over either an exact ring or ``SR``.
 371
 372     - ``K`` -- A polyhedral closed convex cone.
 373
 374     OUTPUT:
 375
 376     If the base ring of ``L`` is exact, then ``True`` will be returned if
 377     and only if ``L`` is Lyapunov-like on ``K``.
 378
 379     If the base ring of ``L`` is ``SR``, then the situation is more
 380     complicated:
 381
 382     - ``True`` will be returned if it can be proven that ``L``
 383       is Lyapunov-like on ``K``.
 384     - ``False`` will be returned if it can be proven that ``L``
 385       is not Lyapunov-like on ``K``.
 386     - ``False`` will also be returned if we can't decide; specifically
 387       if we arrive at a symbolic inequality that cannot be resolved.
 388
 389     .. SEEALSO::
 390
 391           :func:`is_positive_on`,
 392           :func:`is_cross_positive_on`,
 393           :func:`is_Z_operator_on`
 394
 395     EXAMPLES:
 396
 397     Diagonal matrices are Lyapunov-like operators on the nonnegative
 398     orthant::
 399
 400         sage: K = Cone([(1,0,0),(0,1,0),(0,0,1)])
 401         sage: L = diagonal_matrix(random_vector(QQ,3))
 402         sage: is_lyapunov_like_on(L,K)
 403         True
 404
 405     TESTS:
 406
 407     The identity operator is always Lyapunov-like::
 408
 409         sage: set_random_seed()
 410         sage: K = random_cone(max_ambient_dim=8)
 411         sage: L = identity_matrix(K.lattice_dim())
 412         sage: is_lyapunov_like_on(L,K)
 413         True
 414
 415     The "zero" operator is always Lyapunov-like::
 416
 417         sage: K = random_cone(max_ambient_dim=8)
 418         sage: R = K.lattice().vector_space().base_ring()
 419         sage: L = zero_matrix(R, K.lattice_dim())
 420         sage: is_lyapunov_like_on(L,K)
 421         True
 422
 423     Everything in ``K.lyapunov_like_basis()`` should be Lyapunov-like
 424     on ``K``::
 425
 426         sage: K = random_cone(max_ambient_dim=5)
 427         sage: all([ is_lyapunov_like_on(L,K)            # long time
 428         ....:       for L in K.lyapunov_like_basis() ]) # long time
 429         True
 430         sage: all([ is_lyapunov_like_on(L.change_ring(SR),K) # long time
 431         ....:       for L in K.lyapunov_like_basis() ])      # long time
 432         True
 433
 434     Technically we could test this, but for now only closed convex cones
 435     are supported as our ``K`` argument::
 436
 437         sage: L = identity_matrix(3)
 438         sage: K = [ vector([2,2,-1]), vector([5,4,-3]) ]
 439         sage: is_lyapunov_like_on(L,K)
 440         Traceback (most recent call last):
 441         ...
 442         TypeError: K must be a Cone.
 443
 444     We can't give reliable answers over inexact rings::
 445
 446         sage: K = Cone([(1,2,3), (4,5,6)])
 447         sage: L = identity_matrix(RR,3)
 448         sage: is_lyapunov_like_on(L,K)
 449         Traceback (most recent call last):
 450         ...
 451         ValueError: The base ring of L is neither SR nor exact.
 452
 453     An operator is Lyapunov-like on a cone if and only if both the
 454     operator and its negation are cross-positive on the cone::
 455
 456         sage: K = random_cone(max_ambient_dim=5)
 457         sage: R = K.lattice().vector_space().base_ring()
 458         sage: L = random_matrix(R, K.lattice_dim())
 459         sage: actual = is_lyapunov_like_on(L,K)          # long time
 460         sage: expected = (is_cross_positive_on(L,K) and  # long time
 461         ....:             is_cross_positive_on(-L,K))    # long time
 462         sage: actual == expected                         # long time
 463         True
 464
 465     """
 466     if not is_Cone(K):
 467         raise TypeError('K must be a Cone.')
 468     if not L.base_ring().is_exact() and not L.base_ring() is SR:
 469         raise ValueError('The base ring of L is neither SR nor exact.')
 470
 471     # Even though ``discrete_complementarity_set`` is a cached method
 472     # of cones, this is faster than calling ``is_cross_positive_on``
 473     # twice: doing so checks twice as many inequalities as the number
 474     # of equalities that we're about to check.
 475     return all([ s*(L*x) == 0
 476                  for (x,s) in K.discrete_complementarity_set() ])
 477
 478
 479 def LL_cone(K):
 480     gens = K.lyapunov_like_basis()
 481     L = ToricLattice(K.lattice_dim()**2)
 482     return Cone([ g.list() for g in gens ], lattice=L, check=False)
 483
 484 def Sigma_cone(K):
 485     gens = K.cross_positive_operators_gens()
 486     L = ToricLattice(K.lattice_dim()**2)
 487     return Cone([ g.list() for g in gens ], lattice=L, check=False)
 488
 489 def Z_cone(K):
 490     gens = K.Z_operators_gens()
 491     L = ToricLattice(K.lattice_dim()**2)
 492     return Cone([ g.list() for g in gens ], lattice=L, check=False)
 493
 494 def pi_cone(K1, K2=None):
 495     if K2 is None:
 496         K2 = K1
 497     gens = K1.positive_operators_gens(K2)
 498     L = ToricLattice(K1.lattice_dim()*K2.lattice_dim())
 499     return Cone([ g.list() for g in gens ], lattice=L, check=False)