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         """
 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 ``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 ``1/sqrt(2)`` times the "horizontal"
 418         distance from ``point`` to boundary of this cone. For simplicity
 419         we take ``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`` is one (its first coordinate), and so the
 450         radius of the circle obtained from a height-one cross section is
 451         also one. Note that the last two coordinates of ``x`` are half
 452         of the way to the boundary of the cone, and in the direction of
 453         a 30-60-90 triangle. If one follows those coordinates, they hit
 454         at ``(1, sqrt(3)/2, 1/2)`` having unit norm. Thus the
 455         "horizontal" distance to the boundary of the cone is ``(1 -
 456         norm(x)``, which simplifies to ``1/2``. And rather than involve
 457         a square root, we divide by two for a final safe radius of
 458         ``1/4``.
 459
 460             >>> from math import sqrt
 461             >>> K = IceCream(3)
 462             >>> x = matrix([1, sqrt(3)/4.0, 1/4.0])
 463             >>> K.ball_radius(x)
 464             0.25
 465
 466         """
 467         if not isinstance(point, matrix):
 468             raise TypeError('the given point is not a cvxopt.base.matrix')
 469         if not point.size == (self.dimension(), 1):
 470             raise TypeError('the given point has the wrong dimensions')
 471         if not point in self:
 472             raise ValueError('the given point does not lie in the cone')
 473
 474         height = point[0]
 475         radius = norm(point[1:])
 476         return (height - radius) / 2.0
 477
 478
 479 class SymmetricPSD(SymmetricCone):
 480     r"""
 481     The cone of real symmetric positive-semidefinite matrices.
 482
 483     This cone has a dimension ``n`` associated with it, but we let ``n``
 484     refer to the dimension of the domain of our matrices and not the
 485     dimension of the (much larger) space in which the matrices
 486     themselves live. In other words, our ``n`` is the ``n`` that appears
 487     in the usual notation :math:`S^{n}` for symmetric matrices.
 488
 489     As a result, the cone ``SymmetricPSD(n)`` lives in a space of dimension
 490     :math:`\left(n^{2} + n\right)/2)`.
 491
 492     Examples
 493     --------
 494
 495         >>> K = SymmetricPSD(3)
 496         >>> print(K)
 497         Cone of symmetric positive-semidefinite matrices on the real 3-space
 498         >>> K.dimension()
 499         3
 500
 501     """
 502     def __str__(self):
 503         """
 504         Output a human-readable description of myself.
 505         """
 506         tpl = 'Cone of symmetric positive-semidefinite matrices ' \
 507               'on the real {:d}-space'
 508         return tpl.format(self.dimension())
 509
 510
 511     def __contains__(self, point):
 512         """
 513         Return whether or not ``point`` belongs to this cone.
 514
 515         Since this test is expected to work on points whose components
 516         are floating point numbers, it doesn't make any sense to
 517         distinguish between strict and non-strict containment -- the
 518         test uses a tolerance parameter.
 519
 520         Parameters
 521         ----------
 522
 523         point : matrix
 524             A :class:`cvxopt.base.matrix` having dimensions ``(n,n)``
 525             where ``n`` is the :meth:`dimension` of this cone.
 526
 527         Returns
 528         -------
 529
 530         bool
 531
 532            ``True`` if ``point`` belongs to this cone, ``False`` otherwise.
 533
 534         Raises
 535         ------
 536
 537         TypeError
 538             If ``point`` is not a :class:`cvxopt.base.matrix`.
 539
 540         TypeError
 541             If ``point`` has the wrong dimensions.
 542
 543         Examples
 544         --------
 545
 546         These all lie in the interior of the Symmetric PSD cone:
 547
 548             >>> K = SymmetricPSD(2)
 549             >>> matrix([[1,0],[0,1]]) in K
 550             True
 551
 552             >>> K = SymmetricPSD(3)
 553             >>> matrix([[2,-1,0],[-1,2,-1],[0,-1,2]]) in K
 554             True
 555
 556             >>> K = SymmetricPSD(5)
 557             >>> A = matrix([[5,4,3,2,1],
 558             ...             [4,5,4,3,2],
 559             ...             [3,4,5,4,3],
 560             ...             [2,3,4,5,4],
 561             ...             [1,2,3,4,5]])
 562             >>> A in K
 563             True
 564
 565         Boundary points lie in the cone as well:
 566
 567             >>> K = SymmetricPSD(2)
 568             >>> matrix([[0,0],[0,0]]) in K
 569             True
 570
 571             >>> K = SymmetricPSD(5)
 572             >>> A = matrix([[1,0,0,0,0],
 573             ...             [0,1,0,0,0],
 574             ...             [0,0,0,0,0],
 575             ...             [0,0,0,1,0],
 576             ...             [0,0,0,0,1]])
 577             >>> A in K
 578             True
 579
 580         However, this matrix has a negative eigenvalue:
 581
 582            >>> K = SymmetricPSD(2)
 583            >>> A = matrix([[ 1, -2],
 584            ...             [-2,  1]])
 585            >>> A in K
 586            False
 587
 588         An asymmetric cone with positive eigenvalues is not in the cone:
 589
 590            >>> K = SymmetricPSD(2)
 591            >>> A = matrix([[10, 2],
 592            ...             [4,  8]])
 593            >>> A in K
 594            False
 595
 596         Junk arguments don't work:
 597
 598             >>> K = SymmetricPSD(2)
 599             >>> [[1,2],[2,3]] in K
 600             Traceback (most recent call last):
 601             ...
 602             TypeError: the given point is not a cvxopt.base.matrix
 603
 604             >>> K = SymmetricPSD(3)
 605             >>> matrix([[1,2],[3,4]]) in K
 606             Traceback (most recent call last):
 607             ...
 608             TypeError: the given point has the wrong dimensions
 609
 610         """
 611         if not isinstance(point, matrix):
 612             raise TypeError('the given point is not a cvxopt.base.matrix')
 613         if not point.size == (self.dimension(), self.dimension()):
 614             raise TypeError('the given point has the wrong dimensions')
 615         if not point.typecode == 'd':
 616             point = matrix(point, (self.dimension(), self.dimension()), 'd')
 617         if not norm(point - point.trans()) < options.ABS_TOL:
 618             # It's not symmetric.
 619             return False
 620         return all([e > -options.ABS_TOL for e in eigenvalues(point)])
 621
 622
 623
 624 class CartesianProduct(SymmetricCone):
 625     """
 626     A cartesian product of symmetric cones, which is itself a symmetric
 627     cone.
 628
 629     Examples
 630     --------
 631
 632         >>> K = CartesianProduct(NonnegativeOrthant(3), IceCream(2))
 633         >>> print(K)
 634         Cartesian product of dimension 5 with 2 factors:
 635           * Nonnegative orthant in the real 3-space
 636           * Lorentz "ice cream" cone in the real 2-space
 637
 638     """
 639     def __init__(self, *factors):
 640         my_dimension = sum([f.dimension() for f in factors])
 641         super().__init__(my_dimension)
 642         self._factors = factors
 643
 644
 645     def __str__(self):
 646         """
 647         Output a human-readable description of myself.
 648         """
 649         tpl = 'Cartesian product of dimension {:d} with {:d} factors:'
 650         tpl += '\n  * {!s}' * len(self.factors())
 651         format_args = [self.dimension(), len(self.factors())]
 652         format_args += list(self.factors())
 653         return tpl.format(*format_args)
 654
 655
 656     def __contains__(self, point):
 657         """
 658         Return whether or not ``point`` belongs to this cone.
 659
 660         The ``point`` is expected to be a tuple of points which will be
 661         tested for membership in this cone's factors. If each point in
 662         the tuple belongs to its corresponding factor, then the whole
 663         point belongs to this cone. Otherwise, it doesn't.
 664
 665         Parameters
 666         ----------
 667
 668         point : tuple of matrix
 669             A tuple of :class:`cvxopt.base.matrix` corresponding to the
 670             :meth:`factors` of this cartesian product.
 671
 672         Returns
 673         -------
 674
 675         bool
 676
 677            ``True`` if ``point`` belongs to this cone, ``False`` otherwise.
 678
 679         Raises
 680         ------
 681
 682         TypeError
 683             If ``point`` is not a tuple of :class:`cvxopt.base.matrix`.
 684
 685         TypeError
 686             If any element of ``point`` has the wrong dimensions.
 687
 688         Examples
 689         --------
 690
 691         The result depends on how containment is defined for our factors:
 692
 693             >>> K = CartesianProduct(NonnegativeOrthant(3), IceCream(3))
 694             >>> (matrix([1,2,3]), matrix([1,0.5,0.5])) in K
 695             True
 696
 697             >>> K = CartesianProduct(NonnegativeOrthant(3), IceCream(3))
 698             >>> (matrix([0,0,0]), matrix([1,0,1])) in K
 699             True
 700
 701             >>> K = CartesianProduct(NonnegativeOrthant(3), IceCream(3))
 702             >>> (matrix([1,1,1]), matrix([1,1,1])) in K
 703             False
 704
 705             >>> K = CartesianProduct(NonnegativeOrthant(3), IceCream(3))
 706             >>> (matrix([1,-1,1]), matrix([1,0,1])) in K
 707             False
 708
 709         Junk arguments don't work:
 710
 711             >>> K = CartesianProduct(NonnegativeOrthant(3), IceCream(3))
 712             >>> [[1,2,3],[4,5,6]] in K
 713             Traceback (most recent call last):
 714             ...
 715             TypeError: the given point is not a cvxopt.base.matrix
 716
 717             >>> K = CartesianProduct(NonnegativeOrthant(3), IceCream(3))
 718             >>> (matrix([1,2]), matrix([3,4,5,6])) in K
 719             Traceback (most recent call last):
 720             ...
 721             TypeError: the given point has the wrong dimensions
 722
 723         """
 724         return all([p in f for (p, f) in zip(point, self.factors())])
 725
 726
 727
 728     def factors(self):
 729         """
 730         Return a tuple containing the factors (in order) of this
 731         cartesian product.
 732
 733         Returns
 734         -------
 735
 736         tuple of :class:`SymmetricCone`.
 737             The factors of this cartesian product.
 738
 739         Examples
 740         --------
 741
 742             >>> K = CartesianProduct(NonnegativeOrthant(3), IceCream(2))
 743             >>> len(K.factors())
 744             2
 745
 746         """
 747         return self._factors
 748
 749
 750     def cvxopt_dims(self):
 751         """
 752         Return a dictionary of dimensions corresponding to the factors
 753         of this cartesian product. The format of this dictionary is
 754         described in the CVXOPT user's guide:
 755
 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