]> gitweb.michael.orlitzky.com - dunshire.git/blob - dunshire/games.py
e9ae21dd9fca7803ce0f67d5ba60f3d3555419c7
[dunshire.git] / dunshire / games.py
1 """
2 Symmetric linear games and their solutions.
3
4 This module contains the main :class:`SymmetricLinearGame` class that
5 knows how to solve a linear game.
6 """
7 from cvxopt import matrix, printing, solvers
8 from .cones import CartesianProduct
9 from .errors import GameUnsolvableException, PoorScalingException
10 from .matrices import (append_col, append_row, condition_number, identity,
11 inner_product, norm, specnorm)
12 from .options import ABS_TOL, FLOAT_FORMAT, DEBUG_FLOAT_FORMAT
13
14 printing.options['dformat'] = FLOAT_FORMAT
15
16
17 class Solution:
18 """
19 A representation of the solution of a linear game. It should contain
20 the value of the game, and both players' strategies.
21
22 Examples
23 --------
24
25 >>> print(Solution(10, matrix([1,2]), matrix([3,4])))
26 Game value: 10.000...
27 Player 1 optimal:
28 [ 1]
29 [ 2]
30 Player 2 optimal:
31 [ 3]
32 [ 4]
33
34 """
35 def __init__(self, game_value, p1_optimal, p2_optimal):
36 """
37 Create a new Solution object from a game value and two optimal
38 strategies for the players.
39 """
40 self._game_value = game_value
41 self._player1_optimal = p1_optimal
42 self._player2_optimal = p2_optimal
43
44 def __str__(self):
45 """
46 Return a string describing the solution of a linear game.
47
48 The three data that are described are,
49
50 * The value of the game.
51 * The optimal strategy of player one.
52 * The optimal strategy of player two.
53
54 The two optimal strategy vectors are indented by two spaces.
55 """
56 tpl = 'Game value: {:.7f}\n' \
57 'Player 1 optimal:{:s}\n' \
58 'Player 2 optimal:{:s}'
59
60 p1_str = '\n{!s}'.format(self.player1_optimal())
61 p1_str = '\n '.join(p1_str.splitlines())
62 p2_str = '\n{!s}'.format(self.player2_optimal())
63 p2_str = '\n '.join(p2_str.splitlines())
64
65 return tpl.format(self.game_value(), p1_str, p2_str)
66
67
68 def game_value(self):
69 """
70 Return the game value for this solution.
71
72 Examples
73 --------
74
75 >>> s = Solution(10, matrix([1,2]), matrix([3,4]))
76 >>> s.game_value()
77 10
78
79 """
80 return self._game_value
81
82
83 def player1_optimal(self):
84 """
85 Return player one's optimal strategy in this solution.
86
87 Examples
88 --------
89
90 >>> s = Solution(10, matrix([1,2]), matrix([3,4]))
91 >>> print(s.player1_optimal())
92 [ 1]
93 [ 2]
94 <BLANKLINE>
95
96 """
97 return self._player1_optimal
98
99
100 def player2_optimal(self):
101 """
102 Return player two's optimal strategy in this solution.
103
104 Examples
105 --------
106
107 >>> s = Solution(10, matrix([1,2]), matrix([3,4]))
108 >>> print(s.player2_optimal())
109 [ 3]
110 [ 4]
111 <BLANKLINE>
112
113 """
114 return self._player2_optimal
115
116
117 class SymmetricLinearGame:
118 r"""
119 A representation of a symmetric linear game.
120
121 The data for a symmetric linear game are,
122
123 * A "payoff" operator ``L``.
124 * A symmetric cone ``K``.
125 * Two points ``e1`` and ``e2`` in the interior of ``K``.
126
127 The ambient space is assumed to be the span of ``K``.
128
129 With those data understood, the game is played as follows. Players
130 one and two choose points :math:`x` and :math:`y` respectively, from
131 their respective strategy sets,
132
133 .. math::
134 \begin{aligned}
135 \Delta_{1}
136 &=
137 \left\{
138 x \in K \ \middle|\ \left\langle x, e_{2} \right\rangle = 1
139 \right\}\\
140 \Delta_{2}
141 &=
142 \left\{
143 y \in K \ \middle|\ \left\langle y, e_{1} \right\rangle = 1
144 \right\}.
145 \end{aligned}
146
147 Afterwards, a "payout" is computed as :math:`\left\langle
148 L\left(x\right), y \right\rangle` and is paid to player one out of
149 player two's pocket. The game is therefore zero sum, and we suppose
150 that player one would like to guarantee himself the largest minimum
151 payout possible. That is, player one wishes to,
152
153 .. math::
154 \begin{aligned}
155 \text{maximize }
156 &\underset{y \in \Delta_{2}}{\min}\left(
157 \left\langle L\left(x\right), y \right\rangle
158 \right)\\
159 \text{subject to } & x \in \Delta_{1}.
160 \end{aligned}
161
162 Player two has the simultaneous goal to,
163
164 .. math::
165 \begin{aligned}
166 \text{minimize }
167 &\underset{x \in \Delta_{1}}{\max}\left(
168 \left\langle L\left(x\right), y \right\rangle
169 \right)\\
170 \text{subject to } & y \in \Delta_{2}.
171 \end{aligned}
172
173 These goals obviously conflict (the game is zero sum), but an
174 existence theorem guarantees at least one optimal min-max solution
175 from which neither player would like to deviate. This class is
176 able to find such a solution.
177
178 Parameters
179 ----------
180
181 L : list of list of float
182 A matrix represented as a list of **rows**. This representation
183 agrees with (for example) `SageMath <http://www.sagemath.org/>`_
184 and `NumPy <http://www.numpy.org/>`_, but not with CVXOPT (whose
185 matrix constructor accepts a list of columns). In reality, ``L``
186 can be any iterable type of the correct length; however, you
187 should be extremely wary of the way we interpret anything other
188 than a list of rows.
189
190 K : dunshire.cones.SymmetricCone
191 The symmetric cone instance over which the game is played.
192
193 e1 : iterable float
194 The interior point of ``K`` belonging to player one; it
195 can be of any iterable type having the correct length.
196
197 e2 : iterable float
198 The interior point of ``K`` belonging to player two; it
199 can be of any enumerable type having the correct length.
200
201 Raises
202 ------
203
204 ValueError
205 If either ``e1`` or ``e2`` lie outside of the cone ``K``.
206
207 Examples
208 --------
209
210 >>> from dunshire import *
211 >>> K = NonnegativeOrthant(3)
212 >>> L = [[1,-5,-15],[-1,2,-3],[-12,-15,1]]
213 >>> e1 = [1,1,1]
214 >>> e2 = [1,2,3]
215 >>> SLG = SymmetricLinearGame(L, K, e1, e2)
216 >>> print(SLG)
217 The linear game (L, K, e1, e2) where
218 L = [ 1 -5 -15]
219 [ -1 2 -3]
220 [-12 -15 1],
221 K = Nonnegative orthant in the real 3-space,
222 e1 = [ 1]
223 [ 1]
224 [ 1],
225 e2 = [ 1]
226 [ 2]
227 [ 3]
228
229 Lists can (and probably should) be used for every argument::
230
231 >>> from dunshire import *
232 >>> K = NonnegativeOrthant(2)
233 >>> L = [[1,0],[0,1]]
234 >>> e1 = [1,1]
235 >>> e2 = [1,1]
236 >>> G = SymmetricLinearGame(L, K, e1, e2)
237 >>> print(G)
238 The linear game (L, K, e1, e2) where
239 L = [ 1 0]
240 [ 0 1],
241 K = Nonnegative orthant in the real 2-space,
242 e1 = [ 1]
243 [ 1],
244 e2 = [ 1]
245 [ 1]
246
247 The points ``e1`` and ``e2`` can also be passed as some other
248 enumerable type (of the correct length) without much harm, since
249 there is no row/column ambiguity::
250
251 >>> import cvxopt
252 >>> import numpy
253 >>> from dunshire import *
254 >>> K = NonnegativeOrthant(2)
255 >>> L = [[1,0],[0,1]]
256 >>> e1 = cvxopt.matrix([1,1])
257 >>> e2 = numpy.matrix([1,1])
258 >>> G = SymmetricLinearGame(L, K, e1, e2)
259 >>> print(G)
260 The linear game (L, K, e1, e2) where
261 L = [ 1 0]
262 [ 0 1],
263 K = Nonnegative orthant in the real 2-space,
264 e1 = [ 1]
265 [ 1],
266 e2 = [ 1]
267 [ 1]
268
269 However, ``L`` will always be intepreted as a list of rows, even
270 if it is passed as a :class:`cvxopt.base.matrix` which is
271 otherwise indexed by columns::
272
273 >>> import cvxopt
274 >>> from dunshire import *
275 >>> K = NonnegativeOrthant(2)
276 >>> L = [[1,2],[3,4]]
277 >>> e1 = [1,1]
278 >>> e2 = e1
279 >>> G = SymmetricLinearGame(L, K, e1, e2)
280 >>> print(G)
281 The linear game (L, K, e1, e2) where
282 L = [ 1 2]
283 [ 3 4],
284 K = Nonnegative orthant in the real 2-space,
285 e1 = [ 1]
286 [ 1],
287 e2 = [ 1]
288 [ 1]
289 >>> L = cvxopt.matrix(L)
290 >>> print(L)
291 [ 1 3]
292 [ 2 4]
293 <BLANKLINE>
294 >>> G = SymmetricLinearGame(L, K, e1, e2)
295 >>> print(G)
296 The linear game (L, K, e1, e2) where
297 L = [ 1 2]
298 [ 3 4],
299 K = Nonnegative orthant in the real 2-space,
300 e1 = [ 1]
301 [ 1],
302 e2 = [ 1]
303 [ 1]
304
305 """
306 def __init__(self, L, K, e1, e2):
307 """
308 Create a new SymmetricLinearGame object.
309 """
310 self._K = K
311 self._e1 = matrix(e1, (K.dimension(), 1))
312 self._e2 = matrix(e2, (K.dimension(), 1))
313
314 # Our input ``L`` is indexed by rows but CVXOPT matrices are
315 # indexed by columns, so we need to transpose the input before
316 # feeding it to CVXOPT.
317 self._L = matrix(L, (K.dimension(), K.dimension())).trans()
318
319 if not self._e1 in K:
320 raise ValueError('the point e1 must lie in the interior of K')
321
322 if not self._e2 in K:
323 raise ValueError('the point e2 must lie in the interior of K')
324
325 # Initial value of cached method.
326 self._L_specnorm_value = None
327
328
329 def __str__(self):
330 """
331 Return a string representation of this game.
332 """
333 tpl = 'The linear game (L, K, e1, e2) where\n' \
334 ' L = {:s},\n' \
335 ' K = {!s},\n' \
336 ' e1 = {:s},\n' \
337 ' e2 = {:s}'
338 indented_L = '\n '.join(str(self.L()).splitlines())
339 indented_e1 = '\n '.join(str(self.e1()).splitlines())
340 indented_e2 = '\n '.join(str(self.e2()).splitlines())
341
342 return tpl.format(indented_L,
343 str(self.K()),
344 indented_e1,
345 indented_e2)
346
347
348 def L(self):
349 """
350 Return the matrix ``L`` passed to the constructor.
351
352 Returns
353 -------
354
355 matrix
356 The matrix that defines this game's :meth:`payoff` operator.
357
358 Examples
359 --------
360
361 >>> from dunshire import *
362 >>> K = NonnegativeOrthant(3)
363 >>> L = [[1,-5,-15],[-1,2,-3],[-12,-15,1]]
364 >>> e1 = [1,1,1]
365 >>> e2 = [1,2,3]
366 >>> SLG = SymmetricLinearGame(L, K, e1, e2)
367 >>> print(SLG.L())
368 [ 1 -5 -15]
369 [ -1 2 -3]
370 [-12 -15 1]
371 <BLANKLINE>
372
373 """
374 return self._L
375
376
377 def K(self):
378 """
379 Return the cone over which this game is played.
380
381 Returns
382 -------
383
384 SymmetricCone
385 The :class:`SymmetricCone` over which this game is played.
386
387 Examples
388 --------
389
390 >>> from dunshire import *
391 >>> K = NonnegativeOrthant(3)
392 >>> L = [[1,-5,-15],[-1,2,-3],[-12,-15,1]]
393 >>> e1 = [1,1,1]
394 >>> e2 = [1,2,3]
395 >>> SLG = SymmetricLinearGame(L, K, e1, e2)
396 >>> print(SLG.K())
397 Nonnegative orthant in the real 3-space
398
399 """
400 return self._K
401
402
403 def e1(self):
404 """
405 Return player one's interior point.
406
407 Returns
408 -------
409
410 matrix
411 The point interior to :meth:`K` affiliated with player one.
412
413 Examples
414 --------
415
416 >>> from dunshire import *
417 >>> K = NonnegativeOrthant(3)
418 >>> L = [[1,-5,-15],[-1,2,-3],[-12,-15,1]]
419 >>> e1 = [1,1,1]
420 >>> e2 = [1,2,3]
421 >>> SLG = SymmetricLinearGame(L, K, e1, e2)
422 >>> print(SLG.e1())
423 [ 1]
424 [ 1]
425 [ 1]
426 <BLANKLINE>
427
428 """
429 return self._e1
430
431
432 def e2(self):
433 """
434 Return player two's interior point.
435
436 Returns
437 -------
438
439 matrix
440 The point interior to :meth:`K` affiliated with player one.
441
442 Examples
443 --------
444
445 >>> from dunshire import *
446 >>> K = NonnegativeOrthant(3)
447 >>> L = [[1,-5,-15],[-1,2,-3],[-12,-15,1]]
448 >>> e1 = [1,1,1]
449 >>> e2 = [1,2,3]
450 >>> SLG = SymmetricLinearGame(L, K, e1, e2)
451 >>> print(SLG.e2())
452 [ 1]
453 [ 2]
454 [ 3]
455 <BLANKLINE>
456
457 """
458 return self._e2
459
460
461 def payoff(self, strategy1, strategy2):
462 r"""
463 Return the payoff associated with ``strategy1`` and ``strategy2``.
464
465 The payoff operator takes pairs of strategies to a real
466 number. For example, if player one's strategy is :math:`x` and
467 player two's strategy is :math:`y`, then the associated payoff
468 is :math:`\left\langle L\left(x\right),y \right\rangle \in
469 \mathbb{R}`. Here, :math:`L` denotes the same linear operator as
470 :meth:`L`. This method computes the payoff given the two
471 players' strategies.
472
473 Parameters
474 ----------
475
476 strategy1 : matrix
477 Player one's strategy.
478
479 strategy2 : matrix
480 Player two's strategy.
481
482 Returns
483 -------
484
485 float
486 The payoff for the game when player one plays ``strategy1``
487 and player two plays ``strategy2``.
488
489 Examples
490 --------
491
492 The value of the game should be the payoff at the optimal
493 strategies::
494
495 >>> from dunshire import *
496 >>> from dunshire.options import ABS_TOL
497 >>> K = NonnegativeOrthant(3)
498 >>> L = [[1,-5,-15],[-1,2,-3],[-12,-15,1]]
499 >>> e1 = [1,1,1]
500 >>> e2 = [1,1,1]
501 >>> SLG = SymmetricLinearGame(L, K, e1, e2)
502 >>> soln = SLG.solution()
503 >>> x_bar = soln.player1_optimal()
504 >>> y_bar = soln.player2_optimal()
505 >>> abs(SLG.payoff(x_bar, y_bar) - soln.game_value()) < ABS_TOL
506 True
507
508 """
509 return inner_product(self.L()*strategy1, strategy2)
510
511
512 def dimension(self):
513 """
514 Return the dimension of this game.
515
516 The dimension of a game is not needed for the theory, but it is
517 useful for the implementation. We define the dimension of a game
518 to be the dimension of its underlying cone. Or what is the same,
519 the dimension of the space from which the strategies are chosen.
520
521 Returns
522 -------
523
524 int
525 The dimension of the cone :meth:`K`, or of the space where
526 this game is played.
527
528 Examples
529 --------
530
531 The dimension of a game over the nonnegative quadrant in the
532 plane should be two (the dimension of the plane)::
533
534 >>> from dunshire import *
535 >>> K = NonnegativeOrthant(2)
536 >>> L = [[1,-5],[-1,2]]
537 >>> e1 = [1,1]
538 >>> e2 = [1,4]
539 >>> SLG = SymmetricLinearGame(L, K, e1, e2)
540 >>> SLG.dimension()
541 2
542
543 """
544 return self.K().dimension()
545
546
547 def _zero(self):
548 """
549 Return a column of zeros that fits ``K``.
550
551 This is used in our CVXOPT construction.
552
553 .. warning::
554
555 It is not safe to cache any of the matrices passed to
556 CVXOPT, because it can clobber them.
557
558 Returns
559 -------
560
561 matrix
562 A ``self.dimension()``-by-``1`` column vector of zeros.
563
564 Examples
565 --------
566
567 >>> from dunshire import *
568 >>> K = NonnegativeOrthant(3)
569 >>> L = identity(3)
570 >>> e1 = [1,1,1]
571 >>> e2 = e1
572 >>> SLG = SymmetricLinearGame(L, K, e1, e2)
573 >>> print(SLG._zero())
574 [0.0000000]
575 [0.0000000]
576 [0.0000000]
577 <BLANKLINE>
578
579 """
580 return matrix(0, (self.dimension(), 1), tc='d')
581
582
583 def A(self):
584 r"""
585 Return the matrix ``A`` used in our CVXOPT construction.
586
587 This matrix :math:`A` appears on the right-hand side of :math:`Ax
588 = b` in the statement of the CVXOPT conelp program.
589
590 .. warning::
591
592 It is not safe to cache any of the matrices passed to
593 CVXOPT, because it can clobber them.
594
595 Returns
596 -------
597
598 matrix
599 A ``1``-by-``(1 + self.dimension())`` row vector. Its first
600 entry is zero, and the rest are the entries of :meth:`e2`.
601
602 Examples
603 --------
604
605 >>> from dunshire import *
606 >>> K = NonnegativeOrthant(3)
607 >>> L = [[1,1,1],[1,1,1],[1,1,1]]
608 >>> e1 = [1,1,1]
609 >>> e2 = [1,2,3]
610 >>> SLG = SymmetricLinearGame(L, K, e1, e2)
611 >>> print(SLG.A())
612 [0.0000000 1.0000000 2.0000000 3.0000000]
613 <BLANKLINE>
614
615 """
616 return matrix([0, self.e2()], (1, self.dimension() + 1), 'd')
617
618
619
620 def G(self):
621 r"""
622 Return the matrix ``G`` used in our CVXOPT construction.
623
624 Thus matrix :math:`G` appears on the left-hand side of :math:`Gx
625 + s = h` in the statement of the CVXOPT conelp program.
626
627 .. warning::
628
629 It is not safe to cache any of the matrices passed to
630 CVXOPT, because it can clobber them.
631
632 Returns
633 -------
634
635 matrix
636 A ``2*self.dimension()``-by-``(1 + self.dimension())`` matrix.
637
638 Examples
639 --------
640
641 >>> from dunshire import *
642 >>> K = NonnegativeOrthant(3)
643 >>> L = [[4,5,6],[7,8,9],[10,11,12]]
644 >>> e1 = [1,2,3]
645 >>> e2 = [1,1,1]
646 >>> SLG = SymmetricLinearGame(L, K, e1, e2)
647 >>> print(SLG.G())
648 [ 0.0000000 -1.0000000 0.0000000 0.0000000]
649 [ 0.0000000 0.0000000 -1.0000000 0.0000000]
650 [ 0.0000000 0.0000000 0.0000000 -1.0000000]
651 [ 1.0000000 -4.0000000 -5.0000000 -6.0000000]
652 [ 2.0000000 -7.0000000 -8.0000000 -9.0000000]
653 [ 3.0000000 -10.0000000 -11.0000000 -12.0000000]
654 <BLANKLINE>
655
656 """
657 identity_matrix = identity(self.dimension())
658 return append_row(append_col(self._zero(), -identity_matrix),
659 append_col(self.e1(), -self.L()))
660
661
662 def c(self):
663 r"""
664 Return the vector ``c`` used in our CVXOPT construction.
665
666 The column vector :math:`c` appears in the objective function
667 value :math:`\left\langle c,x \right\rangle` in the statement of
668 the CVXOPT conelp program.
669
670 .. warning::
671
672 It is not safe to cache any of the matrices passed to
673 CVXOPT, because it can clobber them.
674
675 Returns
676 -------
677
678 matrix
679 A :meth:`dimension`-by-``1`` column vector.
680
681 Examples
682 --------
683
684 >>> from dunshire import *
685 >>> K = NonnegativeOrthant(3)
686 >>> L = [[4,5,6],[7,8,9],[10,11,12]]
687 >>> e1 = [1,2,3]
688 >>> e2 = [1,1,1]
689 >>> SLG = SymmetricLinearGame(L, K, e1, e2)
690 >>> print(SLG.c())
691 [-1.0000000]
692 [ 0.0000000]
693 [ 0.0000000]
694 [ 0.0000000]
695 <BLANKLINE>
696
697 """
698 return matrix([-1, self._zero()])
699
700
701 def C(self):
702 """
703 Return the cone ``C`` used in our CVXOPT construction.
704
705 This is the cone over which the conelp program takes place.
706
707 Returns
708 -------
709
710 CartesianProduct
711 The cartesian product of ``K`` with itself.
712
713 Examples
714 --------
715
716 >>> from dunshire import *
717 >>> K = NonnegativeOrthant(3)
718 >>> L = [[4,5,6],[7,8,9],[10,11,12]]
719 >>> e1 = [1,2,3]
720 >>> e2 = [1,1,1]
721 >>> SLG = SymmetricLinearGame(L, K, e1, e2)
722 >>> print(SLG.C())
723 Cartesian product of dimension 6 with 2 factors:
724 * Nonnegative orthant in the real 3-space
725 * Nonnegative orthant in the real 3-space
726
727 """
728 return CartesianProduct(self._K, self._K)
729
730 def h(self):
731 r"""
732 Return the ``h`` vector used in our CVXOPT construction.
733
734 The :math:`h` vector appears on the right-hand side of :math:`Gx + s
735 = h` in the statement of the CVXOPT conelp program.
736
737 .. warning::
738
739 It is not safe to cache any of the matrices passed to
740 CVXOPT, because it can clobber them.
741
742 Returns
743 -------
744
745 matrix
746 A ``2*self.dimension()``-by-``1`` column vector of zeros.
747
748 Examples
749 --------
750
751 >>> from dunshire import *
752 >>> K = NonnegativeOrthant(3)
753 >>> L = [[4,5,6],[7,8,9],[10,11,12]]
754 >>> e1 = [1,2,3]
755 >>> e2 = [1,1,1]
756 >>> SLG = SymmetricLinearGame(L, K, e1, e2)
757 >>> print(SLG.h())
758 [0.0000000]
759 [0.0000000]
760 [0.0000000]
761 [0.0000000]
762 [0.0000000]
763 [0.0000000]
764 <BLANKLINE>
765
766 """
767
768 return matrix([self._zero(), self._zero()])
769
770
771 @staticmethod
772 def b():
773 r"""
774 Return the ``b`` vector used in our CVXOPT construction.
775
776 The vector ``b`` appears on the right-hand side of :math:`Ax =
777 b` in the statement of the CVXOPT conelp program.
778
779 This method is static because the dimensions and entries of
780 ``b`` are known beforehand, and don't depend on any other
781 properties of the game.
782
783 .. warning::
784
785 It is not safe to cache any of the matrices passed to
786 CVXOPT, because it can clobber them.
787
788 Returns
789 -------
790
791 matrix
792 A ``1``-by-``1`` matrix containing a single entry ``1``.
793
794 Examples
795 --------
796
797 >>> from dunshire import *
798 >>> K = NonnegativeOrthant(3)
799 >>> L = [[4,5,6],[7,8,9],[10,11,12]]
800 >>> e1 = [1,2,3]
801 >>> e2 = [1,1,1]
802 >>> SLG = SymmetricLinearGame(L, K, e1, e2)
803 >>> print(SLG.b())
804 [1.0000000]
805 <BLANKLINE>
806
807 """
808 return matrix([1], tc='d')
809
810
811 def player1_start(self):
812 """
813 Return a feasible starting point for player one.
814
815 This starting point is for the CVXOPT formulation and not for
816 the original game. The basic premise is that if you scale
817 :meth:`e2` by the reciprocal of its squared norm, then you get a
818 point in :meth:`K` that makes a unit inner product with
819 :meth:`e2`. We then get to choose the primal objective function
820 value such that the constraint involving :meth:`L` is satisfied.
821
822 Returns
823 -------
824
825 dict
826 A dictionary with two keys, ``'x'`` and ``'s'``, which
827 contain the vectors of the same name in the CVXOPT primal
828 problem formulation.
829
830 The vector ``x`` consists of the primal objective function
831 value concatenated with the strategy (for player one) that
832 achieves it. The vector ``s`` is essentially a dummy
833 variable, and is computed from the equality constraing in
834 the CVXOPT primal problem.
835
836 """
837 p = self.e2() / (norm(self.e2()) ** 2)
838 dist = self.K().ball_radius(self.e1())
839 nu = - self._L_specnorm()/(dist*norm(self.e2()))
840 x = matrix([nu, p], (self.dimension() + 1, 1))
841 s = - self.G()*x
842
843 return {'x': x, 's': s}
844
845
846 def player2_start(self):
847 """
848 Return a feasible starting point for player two.
849
850 This starting point is for the CVXOPT formulation and not for
851 the original game. The basic premise is that if you scale
852 :meth:`e1` by the reciprocal of its squared norm, then you get a
853 point in :meth:`K` that makes a unit inner product with
854 :meth:`e1`. We then get to choose the dual objective function
855 value such that the constraint involving :meth:`L` is satisfied.
856
857 Returns
858 -------
859
860 dict
861 A dictionary with two keys, ``'y'`` and ``'z'``, which
862 contain the vectors of the same name in the CVXOPT dual
863 problem formulation.
864
865 The ``1``-by-``1`` vector ``y`` consists of the dual
866 objective function value. The last :meth:`dimension` entries
867 of the vector ``z`` contain the strategy (for player two)
868 that achieves it. The remaining entries of ``z`` are
869 essentially dummy variables, computed from the equality
870 constraint in the CVXOPT dual problem.
871
872 """
873 q = self.e1() / (norm(self.e1()) ** 2)
874 dist = self.K().ball_radius(self.e2())
875 omega = self._L_specnorm()/(dist*norm(self.e1()))
876 y = matrix([omega])
877 z2 = q
878 z1 = y*self.e2() - self.L().trans()*z2
879 z = matrix([z1, z2], (self.dimension()*2, 1))
880
881 return {'y': y, 'z': z}
882
883
884 def _L_specnorm(self):
885 """
886 Compute the spectral norm of :meth:`L` and cache it.
887
888 The spectral norm of the matrix :meth:`L` is used in a few
889 places. Since it can be expensive to compute, we want to cache
890 its value. That is not possible in :func:`specnorm`, which lies
891 outside of a class, so this is the place to do it.
892
893 Returns
894 -------
895
896 float
897 A nonnegative real number; the largest singular value of
898 the matrix :meth:`L`.
899
900 Examples
901 --------
902
903 >>> from dunshire import *
904 >>> from dunshire.matrices import specnorm
905 >>> L = [[1,2],[3,4]]
906 >>> K = NonnegativeOrthant(2)
907 >>> e1 = [1,1]
908 >>> e2 = e1
909 >>> SLG = SymmetricLinearGame(L,K,e1,e2)
910 >>> specnorm(SLG.L()) == SLG._L_specnorm()
911 True
912
913 """
914 if self._L_specnorm_value is None:
915 self._L_specnorm_value = specnorm(self.L())
916 return self._L_specnorm_value
917
918
919 def tolerance_scale(self, solution):
920 r"""
921 Return a scaling factor that should be applied to :const:`ABS_TOL`
922 for this game.
923
924 When performing certain comparisons, the default tolerance
925 :const:`ABS_TOL` may not be appropriate. For example, if we expect
926 ``x`` and ``y`` to be within :const:`ABS_TOL` of each other,
927 than the inner product of ``L*x`` and ``y`` can be as far apart
928 as the spectral norm of ``L`` times the sum of the norms of
929 ``x`` and ``y``. Such a comparison is made in :meth:`solution`,
930 and in many of our unit tests.
931
932 The returned scaling factor found from the inner product
933 mentioned above is
934
935 .. math::
936
937 \left\lVert L \right\rVert_{2}
938 \left( \left\lVert \bar{x} \right\rVert
939 + \left\lVert \bar{y} \right\rVert
940 \right),
941
942 where :math:`\bar{x}` and :math:`\bar{y}` are optimal solutions
943 for players one and two respectively. This scaling factor is not
944 formally justified, but attempting anything smaller leads to
945 test failures.
946
947 .. warning::
948
949 Optimal solutions are not unique, so the scaling factor
950 obtained from ``solution`` may not work when comparing other
951 solutions.
952
953 Parameters
954 ----------
955
956 solution : Solution
957 A solution of this game, used to obtain the norms of the
958 optimal strategies.
959
960 Returns
961 -------
962
963 float
964 A scaling factor to be multiplied by :const:`ABS_TOL` when
965 making comparisons involving solutions of this game.
966
967 Examples
968 --------
969
970 The spectral norm of ``L`` in this case is around ``5.464``, and
971 the optimal strategies both have norm one, so we expect the
972 tolerance scale to be somewhere around ``2 * 5.464``, or
973 ``10.929``::
974
975 >>> from dunshire import *
976 >>> L = [[1,2],[3,4]]
977 >>> K = NonnegativeOrthant(2)
978 >>> e1 = [1,1]
979 >>> e2 = e1
980 >>> SLG = SymmetricLinearGame(L,K,e1,e2)
981 >>> SLG.tolerance_scale(SLG.solution())
982 10.929...
983
984 """
985 norm_p1_opt = norm(solution.player1_optimal())
986 norm_p2_opt = norm(solution.player2_optimal())
987 scale = self._L_specnorm()*(norm_p1_opt + norm_p2_opt)
988
989 # Don't return anything smaller than 1... we can't go below
990 # out "minimum tolerance."
991 return max(1, scale)
992
993
994 def solution(self):
995 """
996 Solve this linear game and return a :class:`Solution`.
997
998 Returns
999 -------
1000
1001 :class:`Solution`
1002 A :class:`Solution` object describing the game's value and
1003 the optimal strategies of both players.
1004
1005 Raises
1006 ------
1007 GameUnsolvableException
1008 If the game could not be solved (if an optimal solution to its
1009 associated cone program was not found).
1010
1011 PoorScalingException
1012 If the game could not be solved because CVXOPT crashed while
1013 trying to take the square root of a negative number.
1014
1015 Examples
1016 --------
1017
1018 This example is computed in Gowda and Ravindran in the section
1019 "The value of a Z-transformation"::
1020
1021 >>> from dunshire import *
1022 >>> K = NonnegativeOrthant(3)
1023 >>> L = [[1,-5,-15],[-1,2,-3],[-12,-15,1]]
1024 >>> e1 = [1,1,1]
1025 >>> e2 = [1,1,1]
1026 >>> SLG = SymmetricLinearGame(L, K, e1, e2)
1027 >>> print(SLG.solution())
1028 Game value: -6.172...
1029 Player 1 optimal:
1030 [0.551...]
1031 [0.000...]
1032 [0.448...]
1033 Player 2 optimal:
1034 [0.448...]
1035 [0.000...]
1036 [0.551...]
1037
1038 The value of the following game can be computed using the fact
1039 that the identity is invertible::
1040
1041 >>> from dunshire import *
1042 >>> K = NonnegativeOrthant(3)
1043 >>> L = [[1,0,0],[0,1,0],[0,0,1]]
1044 >>> e1 = [1,2,3]
1045 >>> e2 = [4,5,6]
1046 >>> SLG = SymmetricLinearGame(L, K, e1, e2)
1047 >>> print(SLG.solution())
1048 Game value: 0.031...
1049 Player 1 optimal:
1050 [0.031...]
1051 [0.062...]
1052 [0.093...]
1053 Player 2 optimal:
1054 [0.125...]
1055 [0.156...]
1056 [0.187...]
1057
1058 This is another Gowda/Ravindran example that is supposed to have
1059 a negative game value::
1060
1061 >>> from dunshire import *
1062 >>> from dunshire.options import ABS_TOL
1063 >>> L = [[1, -2], [-2, 1]]
1064 >>> K = NonnegativeOrthant(2)
1065 >>> e1 = [1, 1]
1066 >>> e2 = e1
1067 >>> SLG = SymmetricLinearGame(L, K, e1, e2)
1068 >>> SLG.solution().game_value() < -ABS_TOL
1069 True
1070
1071 The following two games are problematic numerically, but we
1072 should be able to solve them::
1073
1074 >>> from dunshire import *
1075 >>> L = [[-0.95237953890954685221, 1.83474556206462535712],
1076 ... [ 1.30481749924621448500, 1.65278664543326403447]]
1077 >>> K = NonnegativeOrthant(2)
1078 >>> e1 = [0.95477167524644313001, 0.63270781756540095397]
1079 >>> e2 = [0.39633793037154141370, 0.10239281495640320530]
1080 >>> SLG = SymmetricLinearGame(L, K, e1, e2)
1081 >>> print(SLG.solution())
1082 Game value: 18.767...
1083 Player 1 optimal:
1084 [0.000...]
1085 [9.766...]
1086 Player 2 optimal:
1087 [1.047...]
1088 [0.000...]
1089
1090 ::
1091
1092 >>> from dunshire import *
1093 >>> L = [[1.54159395026049472754, 2.21344728574316684799],
1094 ... [1.33147433507846657541, 1.17913616272988108769]]
1095 >>> K = NonnegativeOrthant(2)
1096 >>> e1 = [0.39903040089404784307, 0.12377403622479113410]
1097 >>> e2 = [0.15695181142215544612, 0.85527381344651265405]
1098 >>> SLG = SymmetricLinearGame(L, K, e1, e2)
1099 >>> print(SLG.solution())
1100 Game value: 24.614...
1101 Player 1 optimal:
1102 [6.371...]
1103 [0.000...]
1104 Player 2 optimal:
1105 [2.506...]
1106 [0.000...]
1107
1108 This is another one that was difficult numerically, and caused
1109 trouble even after we fixed the first two::
1110
1111 >>> from dunshire import *
1112 >>> L = [[57.22233908627052301199, 41.70631373437460354126],
1113 ... [83.04512571985074487202, 57.82581810406928468637]]
1114 >>> K = NonnegativeOrthant(2)
1115 >>> e1 = [7.31887017043399268346, 0.89744171905822367474]
1116 >>> e2 = [0.11099824781179848388, 6.12564670639315345113]
1117 >>> SLG = SymmetricLinearGame(L,K,e1,e2)
1118 >>> print(SLG.solution())
1119 Game value: 70.437...
1120 Player 1 optimal:
1121 [9.009...]
1122 [0.000...]
1123 Player 2 optimal:
1124 [0.136...]
1125 [0.000...]
1126
1127 And finally, here's one that returns an "optimal" solution, but
1128 whose primal/dual objective function values are far apart::
1129
1130 >>> from dunshire import *
1131 >>> L = [[ 6.49260076597376212248, -0.60528030227678542019],
1132 ... [ 2.59896077096751731972, -0.97685530240286766457]]
1133 >>> K = IceCream(2)
1134 >>> e1 = [1, 0.43749513972645248661]
1135 >>> e2 = [1, 0.46008379832200291260]
1136 >>> SLG = SymmetricLinearGame(L, K, e1, e2)
1137 >>> print(SLG.solution())
1138 Game value: 11.596...
1139 Player 1 optimal:
1140 [ 1.852...]
1141 [-1.852...]
1142 Player 2 optimal:
1143 [ 1.777...]
1144 [-1.777...]
1145
1146 """
1147 try:
1148 opts = {'show_progress': False}
1149 soln_dict = solvers.conelp(self.c(),
1150 self.G(),
1151 self.h(),
1152 self.C().cvxopt_dims(),
1153 self.A(),
1154 self.b(),
1155 primalstart=self.player1_start(),
1156 dualstart=self.player2_start(),
1157 options=opts)
1158 except ValueError as error:
1159 if str(error) == 'math domain error':
1160 # Oops, CVXOPT tried to take the square root of a
1161 # negative number. Report some details about the game
1162 # rather than just the underlying CVXOPT crash.
1163 printing.options['dformat'] = DEBUG_FLOAT_FORMAT
1164 raise PoorScalingException(self)
1165 else:
1166 raise error
1167
1168 # The optimal strategies are named ``p`` and ``q`` in the
1169 # background documentation, and we need to extract them from
1170 # the CVXOPT ``x`` and ``z`` variables. The objective values
1171 # :math:`nu` and :math:`omega` can also be found in the CVXOPT
1172 # ``x`` and ``y`` variables; however, they're stored
1173 # conveniently as separate entries in the solution dictionary.
1174 p1_value = -soln_dict['primal objective']
1175 p2_value = -soln_dict['dual objective']
1176 p1_optimal = soln_dict['x'][1:]
1177 p2_optimal = soln_dict['z'][self.dimension():]
1178
1179 # The "status" field contains "optimal" if everything went
1180 # according to plan. Other possible values are "primal
1181 # infeasible", "dual infeasible", "unknown", all of which mean
1182 # we didn't get a solution.
1183 #
1184 # The "infeasible" ones are the worst, since they indicate
1185 # that CVXOPT is convinced the problem is infeasible (and that
1186 # cannot happen).
1187 if soln_dict['status'] in ['primal infeasible', 'dual infeasible']:
1188 printing.options['dformat'] = DEBUG_FLOAT_FORMAT
1189 raise GameUnsolvableException(self, soln_dict)
1190
1191 # For the game value, we could use any of:
1192 #
1193 # * p1_value
1194 # * p2_value
1195 # * (p1_value + p2_value)/2
1196 # * the game payoff
1197 #
1198 # We want the game value to be the payoff, however, so it
1199 # makes the most sense to just use that, even if it means we
1200 # can't test the fact that p1_value/p2_value are close to the
1201 # payoff.
1202 payoff = self.payoff(p1_optimal, p2_optimal)
1203 soln = Solution(payoff, p1_optimal, p2_optimal)
1204
1205 # The "optimal" and "unknown" results, we actually treat the
1206 # same. Even if CVXOPT bails out due to numerical difficulty,
1207 # it will have some candidate points in mind. If those
1208 # candidates are good enough, we take them. We do the same
1209 # check for "optimal" results.
1210 #
1211 # First we check that the primal/dual objective values are
1212 # close enough because otherwise CVXOPT might return "unknown"
1213 # and give us two points in the cone that are nowhere near
1214 # optimal. And in fact, we need to ensure that they're close
1215 # for "optimal" results, too, because we need to know how
1216 # lenient to be in our testing.
1217 #
1218 if abs(p1_value - p2_value) > self.tolerance_scale(soln)*ABS_TOL:
1219 printing.options['dformat'] = DEBUG_FLOAT_FORMAT
1220 raise GameUnsolvableException(self, soln_dict)
1221
1222 # And we also check that the points it gave us belong to the
1223 # cone, just in case...
1224 if (p1_optimal not in self._K) or (p2_optimal not in self._K):
1225 printing.options['dformat'] = DEBUG_FLOAT_FORMAT
1226 raise GameUnsolvableException(self, soln_dict)
1227
1228 return soln
1229
1230
1231 def condition(self):
1232 r"""
1233 Return the condition number of this game.
1234
1235 In the CVXOPT construction of this game, two matrices ``G`` and
1236 ``A`` appear. When those matrices are nasty, numerical problems
1237 can show up. We define the condition number of this game to be
1238 the average of the condition numbers of ``G`` and ``A`` in the
1239 CVXOPT construction. If the condition number of this game is
1240 high, then you can expect numerical difficulty (such as
1241 :class:`PoorScalingException`).
1242
1243 Returns
1244 -------
1245
1246 float
1247 A real number greater than or equal to one that measures how
1248 bad this game is numerically.
1249
1250 Examples
1251 --------
1252
1253 >>> from dunshire import *
1254 >>> K = NonnegativeOrthant(1)
1255 >>> L = [[1]]
1256 >>> e1 = [1]
1257 >>> e2 = e1
1258 >>> SLG = SymmetricLinearGame(L, K, e1, e2)
1259 >>> SLG.condition()
1260 1.809...
1261
1262 """
1263 return (condition_number(self.G()) + condition_number(self.A()))/2
1264
1265
1266 def dual(self):
1267 r"""
1268 Return the dual game to this game.
1269
1270 If :math:`G = \left(L,K,e_{1},e_{2}\right)` is a linear game,
1271 then its dual is :math:`G^{*} =
1272 \left(L^{*},K^{*},e_{2},e_{1}\right)`. However, since this cone
1273 is symmetric, :math:`K^{*} = K`.
1274
1275 Examples
1276 --------
1277
1278 >>> from dunshire import *
1279 >>> K = NonnegativeOrthant(3)
1280 >>> L = [[1,-5,-15],[-1,2,-3],[-12,-15,1]]
1281 >>> e1 = [1,1,1]
1282 >>> e2 = [1,2,3]
1283 >>> SLG = SymmetricLinearGame(L, K, e1, e2)
1284 >>> print(SLG.dual())
1285 The linear game (L, K, e1, e2) where
1286 L = [ 1 -1 -12]
1287 [ -5 2 -15]
1288 [-15 -3 1],
1289 K = Nonnegative orthant in the real 3-space,
1290 e1 = [ 1]
1291 [ 2]
1292 [ 3],
1293 e2 = [ 1]
1294 [ 1]
1295 [ 1]
1296
1297 """
1298 # We pass ``self.L()`` right back into the constructor, because
1299 # it will be transposed there. And keep in mind that ``self._K``
1300 # is its own dual.
1301 return SymmetricLinearGame(self.L(),
1302 self.K(),
1303 self.e2(),
1304 self.e1())