dunshire/cones.py

   1 """
   2 Class definitions for all of the symmetric cones (and their superclass,
   3 :class:`SymmetricCone`) supported by CVXOPT.
   4 """
   5
   6 from cvxopt import matrix
   7
   8 from .matrices import eigenvalues, norm
   9 from . import options
  10
  11 class SymmetricCone:
  12     """
  13     An instance of a symmetric (self-dual and homogeneous) cone.
  14
  15     There are three types of symmetric cones supported by CVXOPT:
  16
  17       1. The nonnegative orthant in the real n-space.
  18       2. The Lorentz "ice cream" cone, or the second-order cone.
  19       3. The cone of symmetric positive-semidefinite matrices.
  20
  21     This class is intended to encompass them all.
  22
  23     When constructing a single symmetric cone (i.e. not a
  24     :class:`CartesianProduct` of them), the only information that we
  25     need is its dimension. We take that dimension as a parameter, and
  26     store it for later.
  27
  28     Parameters
  29     ----------
  30
  31     dimension : int
  32         The dimension of this cone.
  33
  34     Raises
  35     ------
  36
  37     ValueError
  38         If you try to create a cone with dimension zero or less.
  39
  40     """
  41     def __init__(self, dimension):
  42         """
  43         A generic constructor for symmetric cones.
  44         """
  45         if dimension <= 0:
  46             raise ValueError('cones must have dimension greater than zero')
  47
  48         self._dimension = dimension
  49
  50
  51     def __contains__(self, point):
  52         """
  53         Return whether or not ``point`` belongs to this cone.
  54
  55         Parameters
  56         ----------
  57
  58         point : matrix
  59             The point to test for membership in this cone.
  60
  61         Raises
  62         ------
  63
  64         NotImplementedError
  65             Always, this method must be implemented in subclasses.
  66
  67         Examples
  68         --------
  69
  70             >>> K = SymmetricCone(5)
  71             >>> matrix([1,2]) in K
  72             Traceback (most recent call last):
  73             ...
  74             NotImplementedError
  75
  76         """
  77         raise NotImplementedError
  78
  79
  80     def ball_radius(self, point):
  81         """
  82         Return the radius of a ball around ``point`` in this cone.
  83
  84         Since a radius cannot be negative, the ``point`` must be
  85         contained in this cone; otherwise, an error is raised.
  86
  87         Parameters
  88         ----------
  89
  90         point : matrix
  91             A point contained in this cone.
  92
  93         Returns
  94         -------
  95
  96         float
  97             A radius ``r`` such that the ball of radius ``r`` centered at
  98             ``point`` is contained entirely within this cone.
  99
 100         Raises
 101         ------
 102
 103         NotImplementedError
 104             Always, this method must be implemented in subclasses.
 105
 106         Examples
 107         --------
 108
 109             >>> K = SymmetricCone(5)
 110             >>> K.ball_radius(matrix([1,1,1,1,1]))
 111             Traceback (most recent call last):
 112             ...
 113             NotImplementedError
 114
 115         """
 116         raise NotImplementedError
 117
 118
 119     def dimension(self):
 120         """
 121         Return the dimension of this symmetric cone.
 122
 123         The dimension of this symmetric cone is recorded during its
 124         creation. This method simply returns the recorded value, and
 125         should not need to be overridden in subclasses. We prefer to do
 126         any special computation in ``__init__()`` and record the result
 127         in ``self._dimension``.
 128
 129         Returns
 130         -------
 131
 132         int
 133             The stored dimension (from when this cone was constructed)
 134             of this cone.
 135
 136         Examples
 137         --------
 138
 139             >>> K = SymmetricCone(5)
 140             >>> K.dimension()
 141             5
 142
 143         """
 144         return self._dimension
 145
 146
 147 class NonnegativeOrthant(SymmetricCone):
 148     """
 149     The nonnegative orthant in the given number of dimensions.
 150
 151     Examples
 152     --------
 153
 154         >>> K = NonnegativeOrthant(3)
 155         >>> print(K)
 156         Nonnegative orthant in the real 3-space
 157
 158     """
 159     def __str__(self):
 160         """
 161         Output a human-readable description of myself.
 162         """
 163         tpl = 'Nonnegative orthant in the real {:d}-space'
 164         return tpl.format(self.dimension())
 165
 166
 167     def __contains__(self, point):
 168         """
 169         Return whether or not ``point`` belongs to this cone.
 170
 171         Since this test is expected to work on points whose components
 172         are floating point numbers, it doesn't make any sense to
 173         distinguish between strict and non-strict containment -- the
 174         test uses a tolerance parameter.
 175
 176         Parameters
 177         ----------
 178
 179         point : matrix
 180             A :class:`cvxopt.base.matrix` having dimensions ``(n,1)``
 181             where ``n`` is the :meth:`dimension` of this cone.
 182
 183         Returns
 184         -------
 185
 186         bool
 187
 188            ``True`` if ``point`` belongs to this cone, ``False`` otherwise.
 189
 190         Raises
 191         ------
 192
 193         TypeError
 194             If ``point`` is not a :class:`cvxopt.base.matrix`.
 195
 196         TypeError
 197             If ``point`` has the wrong dimensions.
 198
 199         Examples
 200         --------
 201
 202         All of these coordinates are positive enough:
 203
 204             >>> K = NonnegativeOrthant(3)
 205             >>> matrix([1,2,3]) in K
 206             True
 207
 208         The one negative coordinate pushes this point outside of ``K``:
 209
 210             >>> K = NonnegativeOrthant(3)
 211             >>> matrix([1,-0.1,3]) in K
 212             False
 213
 214         A boundary point is considered inside of ``K``:
 215             >>> K = NonnegativeOrthant(3)
 216             >>> matrix([1,0,3]) in K
 217             True
 218
 219         Junk arguments don't work:
 220
 221             >>> K = NonnegativeOrthant(3)
 222             >>> [1,2,3] in K
 223             Traceback (most recent call last):
 224             ...
 225             TypeError: the given point is not a cvxopt.base.matrix
 226
 227             >>> K = NonnegativeOrthant(3)
 228             >>> matrix([1,2]) in K
 229             Traceback (most recent call last):
 230             ...
 231             TypeError: the given point has the wrong dimensions
 232
 233         """
 234         if not isinstance(point, matrix):
 235             raise TypeError('the given point is not a cvxopt.base.matrix')
 236         if not point.size == (self.dimension(), 1):
 237             raise TypeError('the given point has the wrong dimensions')
 238
 239         return all([x > -options.ABS_TOL for x in point])
 240
 241
 242     def ball_radius(self, point):
 243         """
 244         Return the radius of a ball around ``point`` in this cone.
 245
 246         Since a radius cannot be negative, the ``point`` must be
 247         contained in this cone; otherwise, an error is raised.
 248
 249         The minimum distance from ``point`` to the complement of this
 250         cone will always occur at its projection onto that set. It is
 251         not hard to see that the projection is directly along one of the
 252         coordinates, and so the minimum distance from ``point`` to the
 253         boundary of this cone is the smallest coordinate of
 254         ``point``. Therefore any radius less than that will work; we
 255         divide it in half somewhat arbitrarily.
 256
 257         Parameters
 258         ----------
 259
 260         point : matrix
 261             A point contained in this cone.
 262
 263         Returns
 264         -------
 265
 266         float
 267             A radius ``r`` such that the ball of radius ``r`` centered at
 268             ``point`` is contained entirely within this cone.
 269
 270         Raises
 271         ------
 272
 273         TypeError
 274             If ``point`` is not a :class:`cvxopt.base.matrix`.
 275
 276         TypeError
 277             If ``point`` has the wrong dimensions.
 278
 279         ValueError
 280             if ``point`` is not contained in this cone.
 281
 282         Examples
 283         --------
 284
 285             >>> K = NonnegativeOrthant(5)
 286             >>> K.ball_radius(matrix([1,2,3,4,5]))
 287             0.5
 288
 289         """
 290         if not isinstance(point, matrix):
 291             raise TypeError('the given point is not a cvxopt.base.matrix')
 292         if not point.size == (self.dimension(), 1):
 293             raise TypeError('the given point has the wrong dimensions')
 294         if not point in self:
 295             raise ValueError('the given point does not lie in the cone')
 296
 297         return min(list(point)) / 2.0
 298
 299
 300
 301 class IceCream(SymmetricCone):
 302     """
 303     The Lorentz "ice cream" cone in the given number of dimensions.
 304
 305     Examples
 306     --------
 307
 308         >>> K = IceCream(3)
 309         >>> print(K)
 310         Lorentz "ice cream" cone in the real 3-space
 311
 312     """
 313     def __str__(self):
 314         """
 315         Output a human-readable description of myself.
 316         """
 317         tpl = 'Lorentz "ice cream" cone in the real {:d}-space'
 318         return tpl.format(self.dimension())
 319
 320
 321     def __contains__(self, point):
 322         """
 323         Return whether or not ``point`` belongs to this cone.
 324
 325         Since this test is expected to work on points whose components
 326         are floating point numbers, it doesn't make any sense to
 327         distinguish between strict and non-strict containment -- the
 328         test uses a tolerance parameter.
 329
 330         Parameters
 331         ----------
 332
 333         point : matrix
 334             A :class:`cvxopt.base.matrix` having dimensions ``(n,1)``
 335             where ``n`` is the :meth:`dimension` of this cone.
 336
 337         Returns
 338         -------
 339
 340         bool
 341
 342            ``True`` if ``point`` belongs to this cone, ``False`` otherwise.
 343
 344         Raises
 345         ------
 346
 347         TypeError
 348             If ``point`` is not a :class:`cvxopt.base.matrix`.
 349
 350         TypeError
 351             If ``point`` has the wrong dimensions.
 352
 353         Examples
 354         --------
 355
 356         This point lies well within the ice cream cone:
 357
 358             >>> K = IceCream(3)
 359             >>> matrix([1,0.5,0.5]) in K
 360             True
 361
 362         This one lies on its boundary:
 363
 364             >>> K = IceCream(3)
 365             >>> matrix([1,0,1]) in K
 366             True
 367
 368         This point lies entirely outside of the ice cream cone:
 369
 370             >>> K = IceCream(3)
 371             >>> matrix([1,1,1]) in K
 372             False
 373
 374         Junk arguments don't work:
 375
 376             >>> K = IceCream(3)
 377             >>> [1,2,3] in K
 378             Traceback (most recent call last):
 379             ...
 380             TypeError: the given point is not a cvxopt.base.matrix
 381
 382             >>> K = IceCream(3)
 383             >>> matrix([1,2]) in K
 384             Traceback (most recent call last):
 385             ...
 386             TypeError: the given point has the wrong dimensions
 387
 388         """
 389         if not isinstance(point, matrix):
 390             raise TypeError('the given point is not a cvxopt.base.matrix')
 391         if not point.size == (self.dimension(), 1):
 392             raise TypeError('the given point has the wrong dimensions')
 393
 394         height = point[0]
 395         if self.dimension() == 1:
 396             # In one dimension, the ice cream cone is the nonnegative
 397             # orthant.
 398             return height > -options.ABS_TOL
 399         else:
 400             radius = point[1:]
 401             return norm(radius) < (height + options.ABS_TOL)
 402
 403
 404     def ball_radius(self, point):
 405         r"""
 406         Return the radius of a ball around ``point`` in this cone.
 407
 408         Since a radius cannot be negative, the ``point`` must be
 409         contained in this cone; otherwise, an error is raised.
 410
 411         The minimum distance from ``point`` to the complement of this
 412         cone will always occur at its projection onto that set. It is
 413         not hard to see that the projection is at a "down and out" angle
 414         of :math:`\pi/4` towards the outside of the cone. If one draws a
 415         right triangle involving the ``point`` and that projection, it
 416         becomes clear that the distance between ``point`` and its
 417         projection is a factor of :math:`1/\sqrt(2)` times the "horizontal"
 418         distance from ``point`` to boundary of this cone. For simplicity
 419         we take :math:`1/2` instead.
 420
 421         Parameters
 422         ----------
 423
 424         point : matrix
 425             A point contained in this cone.
 426
 427         Returns
 428         -------
 429
 430         float
 431             A radius ``r`` such that the ball of radius ``r`` centered at
 432             ``point`` is contained entirely within this cone.
 433
 434         Raises
 435         ------
 436
 437         TypeError
 438             If ``point`` is not a :class:`cvxopt.base.matrix`.
 439
 440         TypeError
 441             If ``point`` has the wrong dimensions.
 442
 443         ValueError
 444             if ``point`` is not contained in this cone.
 445
 446         Examples
 447         --------
 448
 449         The height of ``x`` below is one (its first coordinate), and so
 450         the radius of the circle obtained from a height-one cross
 451         section is also one. Note that the last two coordinates of ``x``
 452         are half of the way to the boundary of the cone, and in the
 453         direction of a 30-60-90 triangle. If one follows those
 454         coordinates, they hit at :math:`\left(1, \frac{\sqrt(3)}{2},
 455         \frac{1}{2}\right)` having unit norm. Thus the "horizontal"
 456         distance to the boundary of the cone is :math:`1 - \left\lVert x
 457         \right\rVert`, which simplifies to :math:`1/2`. And rather than
 458         involve a square root, we divide by two for a final safe radius
 459         of :math:`1/4`.
 460
 461             >>> from math import sqrt
 462             >>> K = IceCream(3)
 463             >>> x = matrix([1, sqrt(3)/4.0, 1/4.0])
 464             >>> K.ball_radius(x)
 465             0.25
 466
 467         """
 468         if not isinstance(point, matrix):
 469             raise TypeError('the given point is not a cvxopt.base.matrix')
 470         if not point.size == (self.dimension(), 1):
 471             raise TypeError('the given point has the wrong dimensions')
 472         if not point in self:
 473             raise ValueError('the given point does not lie in the cone')
 474
 475         height = point[0]
 476         radius = norm(point[1:])
 477         return (height - radius) / 2.0
 478
 479
 480 class SymmetricPSD(SymmetricCone):
 481     r"""
 482     The cone of real symmetric positive-semidefinite matrices.
 483
 484     This cone has a dimension ``n`` associated with it, but we let ``n``
 485     refer to the dimension of the domain of our matrices and not the
 486     dimension of the (much larger) space in which the matrices
 487     themselves live. In other words, our ``n`` is the ``n`` that appears
 488     in the usual notation :math:`S^{n}` for symmetric matrices.
 489
 490     As a result, the cone ``SymmetricPSD(n)`` lives in a space of dimension
 491     :math:`\left(n^{2} + n\right)/2)`.
 492
 493     Examples
 494     --------
 495
 496         >>> K = SymmetricPSD(3)
 497         >>> print(K)
 498         Cone of symmetric positive-semidefinite matrices on the real 3-space
 499         >>> K.dimension()
 500         3
 501
 502     """
 503     def __str__(self):
 504         """
 505         Output a human-readable description of myself.
 506         """
 507         tpl = 'Cone of symmetric positive-semidefinite matrices ' \
 508               'on the real {:d}-space'
 509         return tpl.format(self.dimension())
 510
 511
 512     def __contains__(self, point):
 513         """
 514         Return whether or not ``point`` belongs to this cone.
 515
 516         Since this test is expected to work on points whose components
 517         are floating point numbers, it doesn't make any sense to
 518         distinguish between strict and non-strict containment -- the
 519         test uses a tolerance parameter.
 520
 521         Parameters
 522         ----------
 523
 524         point : matrix
 525             A :class:`cvxopt.base.matrix` having dimensions ``(n,n)``
 526             where ``n`` is the :meth:`dimension` of this cone.
 527
 528         Returns
 529         -------
 530
 531         bool
 532
 533            ``True`` if ``point`` belongs to this cone, ``False`` otherwise.
 534
 535         Raises
 536         ------
 537
 538         TypeError
 539             If ``point`` is not a :class:`cvxopt.base.matrix`.
 540
 541         TypeError
 542             If ``point`` has the wrong dimensions.
 543
 544         Examples
 545         --------
 546
 547         These all lie in the interior of the Symmetric PSD cone:
 548
 549             >>> K = SymmetricPSD(2)
 550             >>> matrix([[1,0],[0,1]]) in K
 551             True
 552
 553             >>> K = SymmetricPSD(3)
 554             >>> matrix([[2,-1,0],[-1,2,-1],[0,-1,2]]) in K
 555             True
 556
 557             >>> K = SymmetricPSD(5)
 558             >>> A = matrix([[5,4,3,2,1],
 559             ...             [4,5,4,3,2],
 560             ...             [3,4,5,4,3],
 561             ...             [2,3,4,5,4],
 562             ...             [1,2,3,4,5]])
 563             >>> A in K
 564             True
 565
 566         Boundary points lie in the cone as well:
 567
 568             >>> K = SymmetricPSD(2)
 569             >>> matrix([[0,0],[0,0]]) in K
 570             True
 571
 572             >>> K = SymmetricPSD(5)
 573             >>> A = matrix([[1,0,0,0,0],
 574             ...             [0,1,0,0,0],
 575             ...             [0,0,0,0,0],
 576             ...             [0,0,0,1,0],
 577             ...             [0,0,0,0,1]])
 578             >>> A in K
 579             True
 580
 581         However, this matrix has a negative eigenvalue:
 582
 583            >>> K = SymmetricPSD(2)
 584            >>> A = matrix([[ 1, -2],
 585            ...             [-2,  1]])
 586            >>> A in K
 587            False
 588
 589         An asymmetric cone with positive eigenvalues is not in the cone:
 590
 591            >>> K = SymmetricPSD(2)
 592            >>> A = matrix([[10, 2],
 593            ...             [4,  8]])
 594            >>> A in K
 595            False
 596
 597         Junk arguments don't work:
 598
 599             >>> K = SymmetricPSD(2)
 600             >>> [[1,2],[2,3]] in K
 601             Traceback (most recent call last):
 602             ...
 603             TypeError: the given point is not a cvxopt.base.matrix
 604
 605             >>> K = SymmetricPSD(3)
 606             >>> matrix([[1,2],[3,4]]) in K
 607             Traceback (most recent call last):
 608             ...
 609             TypeError: the given point has the wrong dimensions
 610
 611         """
 612         if not isinstance(point, matrix):
 613             raise TypeError('the given point is not a cvxopt.base.matrix')
 614         if not point.size == (self.dimension(), self.dimension()):
 615             raise TypeError('the given point has the wrong dimensions')
 616         if not point.typecode == 'd':
 617             point = matrix(point, (self.dimension(), self.dimension()), 'd')
 618         if not norm(point - point.trans()) < options.ABS_TOL:
 619             # It's not symmetric.
 620             return False
 621         return all([e > -options.ABS_TOL for e in eigenvalues(point)])
 622
 623
 624
 625 class CartesianProduct(SymmetricCone):
 626     """
 627     A cartesian product of symmetric cones, which is itself a symmetric
 628     cone.
 629
 630     Examples
 631     --------
 632
 633         >>> K = CartesianProduct(NonnegativeOrthant(3), IceCream(2))
 634         >>> print(K)
 635         Cartesian product of dimension 5 with 2 factors:
 636           * Nonnegative orthant in the real 3-space
 637           * Lorentz "ice cream" cone in the real 2-space
 638
 639     """
 640     def __init__(self, *factors):
 641         my_dimension = sum([f.dimension() for f in factors])
 642         super().__init__(my_dimension)
 643         self._factors = factors
 644
 645
 646     def __str__(self):
 647         """
 648         Output a human-readable description of myself.
 649         """
 650         tpl = 'Cartesian product of dimension {:d} with {:d} factors:'
 651         tpl += '\n  * {!s}' * len(self.factors())
 652         format_args = [self.dimension(), len(self.factors())]
 653         format_args += list(self.factors())
 654         return tpl.format(*format_args)
 655
 656
 657     def __contains__(self, point):
 658         """
 659         Return whether or not ``point`` belongs to this cone.
 660
 661         The ``point`` is expected to be a tuple of points which will be
 662         tested for membership in this cone's factors. If each point in
 663         the tuple belongs to its corresponding factor, then the whole
 664         point belongs to this cone. Otherwise, it doesn't.
 665
 666         Parameters
 667         ----------
 668
 669         point : tuple of matrix
 670             A tuple of :class:`cvxopt.base.matrix` corresponding to the
 671             :meth:`factors` of this cartesian product.
 672
 673         Returns
 674         -------
 675
 676         bool
 677
 678            ``True`` if ``point`` belongs to this cone, ``False`` otherwise.
 679
 680         Raises
 681         ------
 682
 683         TypeError
 684             If ``point`` is not a tuple of :class:`cvxopt.base.matrix`.
 685
 686         TypeError
 687             If any element of ``point`` has the wrong dimensions.
 688
 689         Examples
 690         --------
 691
 692         The result depends on how containment is defined for our factors:
 693
 694             >>> K = CartesianProduct(NonnegativeOrthant(3), IceCream(3))
 695             >>> (matrix([1,2,3]), matrix([1,0.5,0.5])) in K
 696             True
 697
 698             >>> K = CartesianProduct(NonnegativeOrthant(3), IceCream(3))
 699             >>> (matrix([0,0,0]), matrix([1,0,1])) in K
 700             True
 701
 702             >>> K = CartesianProduct(NonnegativeOrthant(3), IceCream(3))
 703             >>> (matrix([1,1,1]), matrix([1,1,1])) in K
 704             False
 705
 706             >>> K = CartesianProduct(NonnegativeOrthant(3), IceCream(3))
 707             >>> (matrix([1,-1,1]), matrix([1,0,1])) in K
 708             False
 709
 710         Junk arguments don't work:
 711
 712             >>> K = CartesianProduct(NonnegativeOrthant(3), IceCream(3))
 713             >>> [[1,2,3],[4,5,6]] in K
 714             Traceback (most recent call last):
 715             ...
 716             TypeError: the given point is not a cvxopt.base.matrix
 717
 718             >>> K = CartesianProduct(NonnegativeOrthant(3), IceCream(3))
 719             >>> (matrix([1,2]), matrix([3,4,5,6])) in K
 720             Traceback (most recent call last):
 721             ...
 722             TypeError: the given point has the wrong dimensions
 723
 724         """
 725         return all([p in f for (p, f) in zip(point, self.factors())])
 726
 727
 728
 729     def factors(self):
 730         """
 731         Return a tuple containing the factors (in order) of this
 732         cartesian product.
 733
 734         Returns
 735         -------
 736
 737         tuple of :class:`SymmetricCone`.
 738             The factors of this cartesian product.
 739
 740         Examples
 741         --------
 742
 743             >>> K = CartesianProduct(NonnegativeOrthant(3), IceCream(2))
 744             >>> len(K.factors())
 745             2
 746
 747         """
 748         return self._factors
 749
 750
 751     def cvxopt_dims(self):
 752         """
 753         Return a dictionary of dimensions corresponding to the
 754         factors of this cartesian product. The format of this dictionary
 755         is described in the `CVXOPT user's guide
 756         <http://cvxopt.org/userguide/coneprog.html#linear-cone-programs>`_.
 757
 758         Returns
 759         -------
 760
 761         dict
 762             A dimension dictionary suitable to feed to CVXOPT.
 763
 764         Examples
 765         --------
 766
 767             >>> K = CartesianProduct(NonnegativeOrthant(3),
 768             ...                      IceCream(2),
 769             ...                      IceCream(3))
 770             >>> d = K.cvxopt_dims()
 771             >>> (d['l'], d['q'], d['s'])
 772             (3, [2, 3], [])
 773
 774         """
 775         dims = {'l':0, 'q':[], 's':[]}
 776         dims['l'] += sum([K.dimension()
 777                           for K in self.factors()
 778                           if isinstance(K, NonnegativeOrthant)])
 779         dims['q'] = [K.dimension()
 780                      for K in self.factors()
 781                      if isinstance(K, IceCream)]
 782         dims['s'] = [K.dimension()
 783                      for K in self.factors()
 784                      if isinstance(K, SymmetricPSD)]
 785         return dims
 786
 787