mjo/cone/cone.py

   1 from sage.all import *
   2
   3 def is_lyapunov_like(L,K):
   4     r"""
   5     Determine whether or not ``L`` is Lyapunov-like on ``K``.
   6
   7     We say that ``L`` is Lyapunov-like 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 [Orlitzky]_ 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 Lyapunov-like on ``K``,
  22     and ``False`` otherwise.
  23
  24     .. WARNING::
  25
  26         If this function returns ``True``, then ``L`` is Lyapunov-like
  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         Lyapunov-like 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 zero.
  32
  33     REFERENCES:
  34
  35     M. Orlitzky. The Lyapunov rank of an improper cone.
  36     http://www.optimization-online.org/DB_HTML/2015/10/5135.html
  37
  38     EXAMPLES:
  39
  40     The identity is always Lyapunov-like in a nontrivial space::
  41
  42         sage: set_random_seed()
  43         sage: K = random_cone(min_ambient_dim=1, max_ambient_dim=8)
  44         sage: L = identity_matrix(K.lattice_dim())
  45         sage: is_lyapunov_like(L,K)
  46         True
  47
  48     As is the "zero" transformation::
  49
  50         sage: K = random_cone(min_ambient_dim=1, max_ambient_dim=8)
  51         sage: R = K.lattice().vector_space().base_ring()
  52         sage: L = zero_matrix(R, K.lattice_dim())
  53         sage: is_lyapunov_like(L,K)
  54         True
  55
  56         Everything in ``K.lyapunov_like_basis()`` should be Lyapunov-like
  57         on ``K``::
  58
  59         sage: K = random_cone(min_ambient_dim=1, max_ambient_dim=6)
  60         sage: all([ is_lyapunov_like(L,K) for L in K.lyapunov_like_basis() ])
  61         True
  62
  63     """
  64     return all([(L*x).inner_product(s) == 0
  65                 for (x,s) in K.discrete_complementarity_set()])
  66
  67
  68 def motzkin_decomposition(K):
  69     r"""
  70     Return the pair of components in the Motzkin decomposition of this cone.
  71
  72     Every convex cone is the direct sum of a strictly convex cone and a
  73     linear subspace [Stoer-Witzgall]_. Return a pair ``(P,S)`` of cones
  74     such that ``P`` is strictly convex, ``S`` is a subspace, and ``K``
  75     is the direct sum of ``P`` and ``S``.
  76
  77     OUTPUT:
  78
  79     An ordered pair ``(P,S)`` of closed convex polyhedral cones where
  80     ``P`` is strictly convex, ``S`` is a subspace, and ``K`` is the
  81     direct sum of ``P`` and ``S``.
  82
  83     REFERENCES:
  84
  85     .. [Stoer-Witzgall] J. Stoer and C. Witzgall. Convexity and
  86        Optimization in Finite Dimensions I. Springer-Verlag, New
  87        York, 1970.
  88
  89     EXAMPLES:
  90
  91     The nonnegative orthant is strictly convex, so it is its own
  92     strictly convex component and its subspace component is trivial::
  93
  94         sage: K = Cone([(1,0,0),(0,1,0),(0,0,1)])
  95         sage: (P,S) = motzkin_decomposition(K)
  96         sage: K.is_equivalent(P)
  97         True
  98         sage: S.is_trivial()
  99         True
 100
 101     Likewise, full spaces are their own subspace components::
 102
 103         sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
 104         sage: K.is_full_space()
 105         True
 106         sage: (P,S) = motzkin_decomposition(K)
 107         sage: K.is_equivalent(S)
 108         True
 109         sage: P.is_trivial()
 110         True
 111
 112     TESTS:
 113
 114     A random point in the cone should belong to either the strictly
 115     convex component or the subspace component. If the point is nonzero,
 116     it cannot be in both::
 117
 118         sage: set_random_seed()
 119         sage: K = random_cone(max_ambient_dim=8)
 120         sage: (P,S) = motzkin_decomposition(K)
 121         sage: x = K.random_element(ring=QQ)
 122         sage: P.contains(x) or S.contains(x)
 123         True
 124         sage: x.is_zero() or (P.contains(x) != S.contains(x))
 125         True
 126
 127     The strictly convex component should always be strictly convex, and
 128     the subspace component should always be a subspace::
 129
 130         sage: set_random_seed()
 131         sage: K = random_cone(max_ambient_dim=8)
 132         sage: (P,S) = motzkin_decomposition(K)
 133         sage: P.is_strictly_convex()
 134         True
 135         sage: S.lineality() == S.dim()
 136         True
 137
 138     The generators of the components are obtained from orthogonal
 139     projections of the original generators [Stoer-Witzgall]_::
 140
 141         sage: set_random_seed()
 142         sage: K = random_cone(max_ambient_dim=8)
 143         sage: (P,S) = motzkin_decomposition(K)
 144         sage: A = S.linear_subspace().complement().matrix()
 145         sage: proj_S_perp = A.transpose() * (A*A.transpose()).inverse() * A
 146         sage: expected_P = Cone([ proj_S_perp*g for g in K ], K.lattice())
 147         sage: P.is_equivalent(expected_P)
 148         True
 149         sage: A = S.linear_subspace().matrix()
 150         sage: proj_S = A.transpose() * (A*A.transpose()).inverse() * A
 151         sage: expected_S = Cone([ proj_S*g for g in K ], K.lattice())
 152         sage: S.is_equivalent(expected_S)
 153         True
 154     """
 155     # The lines() method only returns one generator per line. For a true
 156     # line, we also need a generator pointing in the opposite direction.
 157     S_gens = [ direction*gen for direction in [1,-1] for gen in K.lines() ]
 158     S = Cone(S_gens, K.lattice())
 159
 160     # Since ``S`` is a subspace, the rays of its dual generate its
 161     # orthogonal complement.
 162     S_perp = Cone(S.dual(), K.lattice())
 163     P = K.intersection(S_perp)
 164
 165     return (P,S)
 166
 167
 168 def positive_operator_gens(K):
 169     r"""
 170     Compute generators of the cone of positive operators on this cone.
 171
 172     OUTPUT:
 173
 174     A list of `n`-by-``n`` matrices where ``n == K.lattice_dim()``.
 175     Each matrix ``P`` in the list should have the property that ``P*x``
 176     is an element of ``K`` whenever ``x`` is an element of
 177     ``K``. Moreover, any nonnegative linear combination of these
 178     matrices shares the same property.
 179
 180     EXAMPLES:
 181
 182     Positive operators on the nonnegative orthant are nonnegative matrices::
 183
 184         sage: K = Cone([(1,)])
 185         sage: positive_operator_gens(K)
 186         [[1]]
 187
 188         sage: K = Cone([(1,0),(0,1)])
 189         sage: positive_operator_gens(K)
 190         [
 191         [1 0]  [0 1]  [0 0]  [0 0]
 192         [0 0], [0 0], [1 0], [0 1]
 193         ]
 194
 195     The trivial cone in a trivial space has no positive operators::
 196
 197         sage: K = Cone([], ToricLattice(0))
 198         sage: positive_operator_gens(K)
 199         []
 200
 201     Every operator is positive on the trivial cone::
 202
 203         sage: K = Cone([(0,)])
 204         sage: positive_operator_gens(K)
 205         [[1], [-1]]
 206
 207         sage: K = Cone([(0,0)])
 208         sage: K.is_trivial()
 209         True
 210         sage: positive_operator_gens(K)
 211         [
 212         [1 0]  [-1  0]  [0 1]  [ 0 -1]  [0 0]  [ 0  0]  [0 0]  [ 0  0]
 213         [0 0], [ 0  0], [0 0], [ 0  0], [1 0], [-1  0], [0 1], [ 0 -1]
 214         ]
 215
 216     Every operator is positive on the ambient vector space::
 217
 218         sage: K = Cone([(1,),(-1,)])
 219         sage: K.is_full_space()
 220         True
 221         sage: positive_operator_gens(K)
 222         [[1], [-1]]
 223
 224         sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
 225         sage: K.is_full_space()
 226         True
 227         sage: positive_operator_gens(K)
 228         [
 229         [1 0]  [-1  0]  [0 1]  [ 0 -1]  [0 0]  [ 0  0]  [0 0]  [ 0  0]
 230         [0 0], [ 0  0], [0 0], [ 0  0], [1 0], [-1  0], [0 1], [ 0 -1]
 231         ]
 232
 233     A non-obvious application is to find the positive operators on the
 234     right half-plane::
 235
 236         sage: K = Cone([(1,0),(0,1),(0,-1)])
 237         sage: positive_operator_gens(K)
 238         [
 239         [1 0]  [0 0]  [ 0  0]  [0 0]  [ 0  0]
 240         [0 0], [1 0], [-1  0], [0 1], [ 0 -1]
 241         ]
 242
 243     TESTS:
 244
 245     Each positive operator generator should send the generators of the
 246     cone into the cone::
 247
 248         sage: set_random_seed()
 249         sage: K = random_cone(max_ambient_dim=5)
 250         sage: pi_of_K = positive_operator_gens(K)
 251         sage: all([ K.contains(P*x) for P in pi_of_K for x in K ])
 252         True
 253
 254     Each positive operator generator should send a random element of the
 255     cone into the cone::
 256
 257         sage: set_random_seed()
 258         sage: K = random_cone(max_ambient_dim=5)
 259         sage: pi_of_K = positive_operator_gens(K)
 260         sage: all([ K.contains(P*K.random_element(QQ)) for P in pi_of_K ])
 261         True
 262
 263     A random element of the positive operator cone should send the
 264     generators of the cone into the cone::
 265
 266         sage: set_random_seed()
 267         sage: K = random_cone(max_ambient_dim=5)
 268         sage: pi_of_K = positive_operator_gens(K)
 269         sage: L = ToricLattice(K.lattice_dim()**2)
 270         sage: pi_cone = Cone([ g.list() for g in pi_of_K ], lattice=L)
 271         sage: P = matrix(K.lattice_dim(), pi_cone.random_element(QQ).list())
 272         sage: all([ K.contains(P*x) for x in K ])
 273         True
 274
 275     A random element of the positive operator cone should send a random
 276     element of the cone into the cone::
 277
 278         sage: set_random_seed()
 279         sage: K = random_cone(max_ambient_dim=5)
 280         sage: pi_of_K = positive_operator_gens(K)
 281         sage: L = ToricLattice(K.lattice_dim()**2)
 282         sage: pi_cone = Cone([ g.list() for g in pi_of_K ], lattice=L)
 283         sage: P = matrix(K.lattice_dim(), pi_cone.random_element(QQ).list())
 284         sage: K.contains(P*K.random_element(ring=QQ))
 285         True
 286
 287     The lineality space of the dual of the cone of positive operators
 288     can be computed from the lineality spaces of the cone and its dual::
 289
 290         sage: set_random_seed()
 291         sage: K = random_cone(max_ambient_dim=5)
 292         sage: pi_of_K = positive_operator_gens(K)
 293         sage: L = ToricLattice(K.lattice_dim()**2)
 294         sage: pi_cone = Cone([ g.list() for g in pi_of_K ], lattice=L)
 295         sage: actual = pi_cone.dual().linear_subspace()
 296         sage: U1 = [ vector((s.tensor_product(x)).list())
 297         ....:        for x in K.lines()
 298         ....:        for s in K.dual() ]
 299         sage: U2 = [ vector((s.tensor_product(x)).list())
 300         ....:        for x in K
 301         ....:        for s in K.dual().lines() ]
 302         sage: expected = pi_cone.lattice().vector_space().span(U1 + U2)
 303         sage: actual == expected
 304         True
 305
 306     The lineality of the dual of the cone of positive operators
 307     is known from its lineality space::
 308
 309         sage: set_random_seed()
 310         sage: K = random_cone(max_ambient_dim=5)
 311         sage: n = K.lattice_dim()
 312         sage: m = K.dim()
 313         sage: l = K.lineality()
 314         sage: pi_of_K = positive_operator_gens(K)
 315         sage: L = ToricLattice(n**2)
 316         sage: pi_cone = Cone([p.list() for p in pi_of_K], lattice=L)
 317         sage: actual = pi_cone.dual().lineality()
 318         sage: expected = l*(m - l) + m*(n - m)
 319         sage: actual == expected
 320         True
 321
 322     The dimension of the cone of positive operators is given by the
 323     corollary in my paper::
 324
 325         sage: set_random_seed()
 326         sage: K = random_cone(max_ambient_dim=5)
 327         sage: n = K.lattice_dim()
 328         sage: m = K.dim()
 329         sage: l = K.lineality()
 330         sage: pi_of_K = positive_operator_gens(K)
 331         sage: L = ToricLattice(n**2)
 332         sage: actual = Cone([p.list() for p in pi_of_K], lattice=L).dim()
 333         sage: expected = n**2 - l*(m - l) - (n - m)*m
 334         sage: actual == expected
 335         True
 336
 337     The trivial cone, full space, and half-plane all give rise to the
 338     expected dimensions::
 339
 340         sage: n = ZZ.random_element().abs()
 341         sage: K = Cone([[0] * n], ToricLattice(n))
 342         sage: K.is_trivial()
 343         True
 344         sage: L = ToricLattice(n^2)
 345         sage: pi_of_K = positive_operator_gens(K)
 346         sage: actual = Cone([p.list() for p in pi_of_K], lattice=L).dim()
 347         sage: actual == n^2
 348         True
 349         sage: K = K.dual()
 350         sage: K.is_full_space()
 351         True
 352         sage: pi_of_K = positive_operator_gens(K)
 353         sage: actual = Cone([p.list() for p in pi_of_K], lattice=L).dim()
 354         sage: actual == n^2
 355         True
 356         sage: K = Cone([(1,0),(0,1),(0,-1)])
 357         sage: pi_of_K = positive_operator_gens(K)
 358         sage: actual = Cone([p.list() for p in pi_of_K]).dim()
 359         sage: actual == 3
 360         True
 361
 362     The lineality of the cone of positive operators follows from the
 363     description of its generators::
 364
 365         sage: set_random_seed()
 366         sage: K = random_cone(max_ambient_dim=5)
 367         sage: n = K.lattice_dim()
 368         sage: pi_of_K = positive_operator_gens(K)
 369         sage: L = ToricLattice(n**2)
 370         sage: actual = Cone([p.list() for p in pi_of_K], lattice=L).lineality()
 371         sage: expected = n**2 - K.dim()*K.dual().dim()
 372         sage: actual == expected
 373         True
 374
 375     The trivial cone, full space, and half-plane all give rise to the
 376     expected linealities::
 377
 378         sage: n = ZZ.random_element().abs()
 379         sage: K = Cone([[0] * n], ToricLattice(n))
 380         sage: K.is_trivial()
 381         True
 382         sage: L = ToricLattice(n^2)
 383         sage: pi_of_K = positive_operator_gens(K)
 384         sage: actual = Cone([p.list() for p in pi_of_K], lattice=L).lineality()
 385         sage: actual == n^2
 386         True
 387         sage: K = K.dual()
 388         sage: K.is_full_space()
 389         True
 390         sage: pi_of_K = positive_operator_gens(K)
 391         sage: actual = Cone([p.list() for p in pi_of_K], lattice=L).lineality()
 392         sage: actual == n^2
 393         True
 394         sage: K = Cone([(1,0),(0,1),(0,-1)])
 395         sage: pi_of_K = positive_operator_gens(K)
 396         sage: actual = Cone([p.list() for p in pi_of_K]).lineality()
 397         sage: actual == 2
 398         True
 399
 400     A cone is proper if and only if its cone of positive operators
 401     is proper::
 402
 403         sage: set_random_seed()
 404         sage: K = random_cone(max_ambient_dim=5)
 405         sage: pi_of_K = positive_operator_gens(K)
 406         sage: L = ToricLattice(K.lattice_dim()**2)
 407         sage: pi_cone = Cone([p.list() for p in pi_of_K], lattice=L)
 408         sage: K.is_proper() == pi_cone.is_proper()
 409         True
 410     """
 411     # Matrices are not vectors in Sage, so we have to convert them
 412     # to vectors explicitly before we can find a basis. We need these
 413     # two values to construct the appropriate "long vector" space.
 414     F = K.lattice().base_field()
 415     n = K.lattice_dim()
 416
 417     tensor_products = [ s.tensor_product(x) for x in K for s in K.dual() ]
 418
 419     # Convert those tensor products to long vectors.
 420     W = VectorSpace(F, n**2)
 421     vectors = [ W(tp.list()) for tp in tensor_products ]
 422
 423     # Create the *dual* cone of the positive operators, expressed as
 424     # long vectors. WARNING: takes forever unless we pass check=False
 425     # to Cone().
 426     pi_dual = Cone(vectors, ToricLattice(W.dimension()))
 427
 428     # Now compute the desired cone from its dual...
 429     pi_cone = pi_dual.dual()
 430
 431     # And finally convert its rays back to matrix representations.
 432     M = MatrixSpace(F, n)
 433     return [ M(v.list()) for v in pi_cone.rays() ]
 434
 435
 436 def Z_transformation_gens(K):
 437     r"""
 438     Compute generators of the cone of Z-transformations on this cone.
 439
 440     OUTPUT:
 441
 442     A list of `n`-by-``n`` matrices where ``n == K.lattice_dim()``.
 443     Each matrix ``L`` in the list should have the property that
 444     ``(L*x).inner_product(s) <= 0`` whenever ``(x,s)`` is an element the
 445     discrete complementarity set of ``K``. Moreover, any nonnegative
 446     linear combination of these matrices shares the same property.
 447
 448     EXAMPLES:
 449
 450     Z-transformations on the nonnegative orthant are just Z-matrices.
 451     That is, matrices whose off-diagonal elements are nonnegative::
 452
 453         sage: K = Cone([(1,0),(0,1)])
 454         sage: Z_transformation_gens(K)
 455         [
 456         [ 0 -1]  [ 0  0]  [-1  0]  [1 0]  [ 0  0]  [0 0]
 457         [ 0  0], [-1  0], [ 0  0], [0 0], [ 0 -1], [0 1]
 458         ]
 459         sage: K = Cone([(1,0,0,0),(0,1,0,0),(0,0,1,0),(0,0,0,1)])
 460         sage: all([ z[i][j] <= 0 for z in Z_transformation_gens(K)
 461         ....:                    for i in range(z.nrows())
 462         ....:                    for j in range(z.ncols())
 463         ....:                    if i != j ])
 464         True
 465
 466     The trivial cone in a trivial space has no Z-transformations::
 467
 468         sage: K = Cone([], ToricLattice(0))
 469         sage: Z_transformation_gens(K)
 470         []
 471
 472     Z-transformations on a subspace are Lyapunov-like and vice-versa::
 473
 474         sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
 475         sage: K.is_full_space()
 476         True
 477         sage: lls = span([ vector(l.list()) for l in K.lyapunov_like_basis() ])
 478         sage: zs  = span([ vector(z.list()) for z in Z_transformation_gens(K) ])
 479         sage: zs == lls
 480         True
 481
 482     TESTS:
 483
 484     The Z-property is possessed by every Z-transformation::
 485
 486         sage: set_random_seed()
 487         sage: K = random_cone(max_ambient_dim=6)
 488         sage: Z_of_K = Z_transformation_gens(K)
 489         sage: dcs = K.discrete_complementarity_set()
 490         sage: all([(z*x).inner_product(s) <= 0 for z in Z_of_K
 491         ....:                                  for (x,s) in dcs])
 492         True
 493
 494     The lineality space of Z is LL::
 495
 496         sage: set_random_seed()
 497         sage: K = random_cone(min_ambient_dim=1, max_ambient_dim=6)
 498         sage: lls = span([ vector(l.list()) for l in K.lyapunov_like_basis() ])
 499         sage: z_cone  = Cone([ z.list() for z in Z_transformation_gens(K) ])
 500         sage: z_cone.linear_subspace() == lls
 501         True
 502
 503     And thus, the lineality of Z is the Lyapunov rank::
 504
 505         sage: set_random_seed()
 506         sage: K = random_cone(max_ambient_dim=6)
 507         sage: Z_of_K = Z_transformation_gens(K)
 508         sage: L = ToricLattice(K.lattice_dim()**2)
 509         sage: z_cone  = Cone([ z.list() for z in Z_of_K ], lattice=L)
 510         sage: z_cone.lineality() == K.lyapunov_rank()
 511         True
 512
 513     The lineality spaces of pi-star and Z-star are equal:
 514
 515         sage: set_random_seed()
 516         sage: K = random_cone(max_ambient_dim=5)
 517         sage: pi_of_K = positive_operator_gens(K)
 518         sage: Z_of_K = Z_transformation_gens(K)
 519         sage: L = ToricLattice(K.lattice_dim()**2)
 520         sage: pi_star = Cone([p.list() for p in pi_of_K], lattice=L).dual()
 521         sage: z_star  = Cone([ z.list() for z in Z_of_K], lattice=L).dual()
 522         sage: pi_star.linear_subspace() == z_star.linear_subspace()
 523         True
 524     """
 525     # Matrices are not vectors in Sage, so we have to convert them
 526     # to vectors explicitly before we can find a basis. We need these
 527     # two values to construct the appropriate "long vector" space.
 528     F = K.lattice().base_field()
 529     n = K.lattice_dim()
 530
 531     # These tensor products contain generators for the dual cone of
 532     # the cross-positive transformations.
 533     tensor_products = [ s.tensor_product(x)
 534                         for (x,s) in K.discrete_complementarity_set() ]
 535
 536     # Turn our matrices into long vectors...
 537     W = VectorSpace(F, n**2)
 538     vectors = [ W(m.list()) for m in tensor_products ]
 539
 540     # Create the *dual* cone of the cross-positive operators,
 541     # expressed as long vectors. WARNING: takes forever unless we pass
 542     # check=False to Cone().
 543     Sigma_dual = Cone(vectors, lattice=ToricLattice(W.dimension()))
 544
 545     # Now compute the desired cone from its dual...
 546     Sigma_cone = Sigma_dual.dual()
 547
 548     # And finally convert its rays back to matrix representations.
 549     # But first, make them negative, so we get Z-transformations and
 550     # not cross-positive ones.
 551     M = MatrixSpace(F, n)
 552     return [ -M(v.list()) for v in Sigma_cone.rays() ]
 553
 554
 555 def Z_cone(K):
 556     gens = Z_transformation_gens(K)
 557     L = None
 558     if len(gens) == 0:
 559         L = ToricLattice(0)
 560     return Cone([ g.list() for g in gens ], lattice=L)
 561
 562 def pi_cone(K):
 563     gens = positive_operator_gens(K)
 564     L = None
 565     if len(gens) == 0:
 566         L = ToricLattice(0)
 567     return Cone([ g.list() for g in gens ], lattice=L)