]> gitweb.michael.orlitzky.com - dunshire.git/blob - dunshire/cones.py
Fix the test script to appease pylint.
[dunshire.git] / 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