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
   8 from cvxopt import matrix, printing, solvers
   9 from .cones import CartesianProduct
  10 from .errors import GameUnsolvableException, PoorScalingException
  11 from .matrices import append_col, append_row, condition_number, identity
  12 from . import options
  13
  14 printing.options['dformat'] = options.FLOAT_FORMAT
  15
  16 class Solution:
  17     """
  18     A representation of the solution of a linear game. It should contain
  19     the value of the game, and both players' strategies.
  20
  21     Examples
  22     --------
  23
  24         >>> print(Solution(10, matrix([1,2]), matrix([3,4])))
  25         Game value: 10.0000000
  26         Player 1 optimal:
  27           [ 1]
  28           [ 2]
  29         Player 2 optimal:
  30           [ 3]
  31           [ 4]
  32
  33     """
  34     def __init__(self, game_value, p1_optimal, p2_optimal):
  35         """
  36         Create a new Solution object from a game value and two optimal
  37         strategies for the players.
  38         """
  39         self._game_value = game_value
  40         self._player1_optimal = p1_optimal
  41         self._player2_optimal = p2_optimal
  42
  43     def __str__(self):
  44         """
  45         Return a string describing the solution of a linear game.
  46
  47         The three data that are described are,
  48
  49           * The value of the game.
  50           * The optimal strategy of player one.
  51           * The optimal strategy of player two.
  52
  53         The two optimal strategy vectors are indented by two spaces.
  54         """
  55         tpl = 'Game value: {:.7f}\n' \
  56               'Player 1 optimal:{:s}\n' \
  57               'Player 2 optimal:{:s}'
  58
  59         p1_str = '\n{!s}'.format(self.player1_optimal())
  60         p1_str = '\n  '.join(p1_str.splitlines())
  61         p2_str = '\n{!s}'.format(self.player2_optimal())
  62         p2_str = '\n  '.join(p2_str.splitlines())
  63
  64         return tpl.format(self.game_value(), p1_str, p2_str)
  65
  66
  67     def game_value(self):
  68         """
  69         Return the game value for this solution.
  70
  71         Examples
  72         --------
  73
  74             >>> s = Solution(10, matrix([1,2]), matrix([3,4]))
  75             >>> s.game_value()
  76             10
  77
  78         """
  79         return self._game_value
  80
  81
  82     def player1_optimal(self):
  83         """
  84         Return player one's optimal strategy in this solution.
  85
  86         Examples
  87         --------
  88
  89             >>> s = Solution(10, matrix([1,2]), matrix([3,4]))
  90             >>> print(s.player1_optimal())
  91             [ 1]
  92             [ 2]
  93             <BLANKLINE>
  94
  95         """
  96         return self._player1_optimal
  97
  98
  99     def player2_optimal(self):
 100         """
 101         Return player two's optimal strategy in this solution.
 102
 103         Examples
 104         --------
 105
 106             >>> s = Solution(10, matrix([1,2]), matrix([3,4]))
 107             >>> print(s.player2_optimal())
 108             [ 3]
 109             [ 4]
 110             <BLANKLINE>
 111
 112         """
 113         return self._player2_optimal
 114
 115
 116 class SymmetricLinearGame:
 117     r"""
 118     A representation of a symmetric linear game.
 119
 120     The data for a symmetric linear game are,
 121
 122       * A "payoff" operator ``L``.
 123       * A symmetric cone ``K``.
 124       * Two points ``e1`` and ``e2`` in the interior of ``K``.
 125
 126     The ambient space is assumed to be the span of ``K``.
 127
 128     With those data understood, the game is played as follows. Players
 129     one and two choose points :math:`x` and :math:`y` respectively, from
 130     their respective strategy sets,
 131
 132     .. math::
 133         \begin{aligned}
 134             \Delta_{1}
 135             &=
 136             \left\{
 137               x \in K \ \middle|\ \left\langle x, e_{2} \right\rangle = 1
 138             \right\}\\
 139             \Delta_{2}
 140             &=
 141             \left\{
 142               y \in K \ \middle|\ \left\langle y, e_{1} \right\rangle = 1
 143             \right\}.
 144         \end{aligned}
 145
 146     Afterwards, a "payout" is computed as :math:`\left\langle
 147     L\left(x\right), y \right\rangle` and is paid to player one out of
 148     player two's pocket. The game is therefore zero sum, and we suppose
 149     that player one would like to guarantee himself the largest minimum
 150     payout possible. That is, player one wishes to,
 151
 152     .. math::
 153         \begin{aligned}
 154             \text{maximize }
 155               &\underset{y \in \Delta_{2}}{\min}\left(
 156                 \left\langle L\left(x\right), y \right\rangle
 157               \right)\\
 158             \text{subject to } & x \in \Delta_{1}.
 159         \end{aligned}
 160
 161     Player two has the simultaneous goal to,
 162
 163     .. math::
 164         \begin{aligned}
 165             \text{minimize }
 166               &\underset{x \in \Delta_{1}}{\max}\left(
 167                 \left\langle L\left(x\right), y \right\rangle
 168               \right)\\
 169             \text{subject to } & y \in \Delta_{2}.
 170         \end{aligned}
 171
 172     These goals obviously conflict (the game is zero sum), but an
 173     existence theorem guarantees at least one optimal min-max solution
 174     from which neither player would like to deviate. This class is
 175     able to find such a solution.
 176
 177     Parameters
 178     ----------
 179
 180     L : list of list of float
 181         A matrix represented as a list of ROWS. This representation
 182         agrees with (for example) SageMath and NumPy, but not with CVXOPT
 183         (whose matrix constructor accepts a list of columns).
 184
 185     K : :class:`SymmetricCone`
 186         The symmetric cone instance over which the game is played.
 187
 188     e1 : iterable float
 189         The interior point of ``K`` belonging to player one; it
 190         can be of any iterable type having the correct length.
 191
 192     e2 : iterable float
 193         The interior point of ``K`` belonging to player two; it
 194         can be of any enumerable type having the correct length.
 195
 196     Raises
 197     ------
 198
 199     ValueError
 200         If either ``e1`` or ``e2`` lie outside of the cone ``K``.
 201
 202     Examples
 203     --------
 204
 205         >>> from dunshire import *
 206         >>> K = NonnegativeOrthant(3)
 207         >>> L = [[1,-5,-15],[-1,2,-3],[-12,-15,1]]
 208         >>> e1 = [1,1,1]
 209         >>> e2 = [1,2,3]
 210         >>> SLG = SymmetricLinearGame(L, K, e1, e2)
 211         >>> print(SLG)
 212         The linear game (L, K, e1, e2) where
 213           L = [  1  -5 -15]
 214               [ -1   2  -3]
 215               [-12 -15   1],
 216           K = Nonnegative orthant in the real 3-space,
 217           e1 = [ 1]
 218                [ 1]
 219                [ 1],
 220           e2 = [ 1]
 221                [ 2]
 222                [ 3],
 223           Condition((L, K, e1, e2)) = 31.834...
 224
 225     Lists can (and probably should) be used for every argument::
 226
 227         >>> from dunshire import *
 228         >>> K = NonnegativeOrthant(2)
 229         >>> L = [[1,0],[0,1]]
 230         >>> e1 = [1,1]
 231         >>> e2 = [1,1]
 232         >>> G = SymmetricLinearGame(L, K, e1, e2)
 233         >>> print(G)
 234         The linear game (L, K, e1, e2) where
 235           L = [ 1  0]
 236               [ 0  1],
 237           K = Nonnegative orthant in the real 2-space,
 238           e1 = [ 1]
 239                [ 1],
 240           e2 = [ 1]
 241                [ 1],
 242           Condition((L, K, e1, e2)) = 1.707...
 243
 244     The points ``e1`` and ``e2`` can also be passed as some other
 245     enumerable type (of the correct length) without much harm, since
 246     there is no row/column ambiguity::
 247
 248         >>> import cvxopt
 249         >>> import numpy
 250         >>> from dunshire import *
 251         >>> K = NonnegativeOrthant(2)
 252         >>> L = [[1,0],[0,1]]
 253         >>> e1 = cvxopt.matrix([1,1])
 254         >>> e2 = numpy.matrix([1,1])
 255         >>> G = SymmetricLinearGame(L, K, e1, e2)
 256         >>> print(G)
 257         The linear game (L, K, e1, e2) where
 258           L = [ 1  0]
 259               [ 0  1],
 260           K = Nonnegative orthant in the real 2-space,
 261           e1 = [ 1]
 262                [ 1],
 263           e2 = [ 1]
 264                [ 1],
 265           Condition((L, K, e1, e2)) = 1.707...
 266
 267     However, ``L`` will always be intepreted as a list of rows, even
 268     if it is passed as a :class:`cvxopt.base.matrix` which is
 269     otherwise indexed by columns::
 270
 271         >>> import cvxopt
 272         >>> from dunshire import *
 273         >>> K = NonnegativeOrthant(2)
 274         >>> L = [[1,2],[3,4]]
 275         >>> e1 = [1,1]
 276         >>> e2 = e1
 277         >>> G = SymmetricLinearGame(L, K, e1, e2)
 278         >>> print(G)
 279         The linear game (L, K, e1, e2) where
 280           L = [ 1  2]
 281               [ 3  4],
 282           K = Nonnegative orthant in the real 2-space,
 283           e1 = [ 1]
 284                [ 1],
 285           e2 = [ 1]
 286                [ 1],
 287           Condition((L, K, e1, e2)) = 6.073...
 288         >>> L = cvxopt.matrix(L)
 289         >>> print(L)
 290         [ 1  3]
 291         [ 2  4]
 292         <BLANKLINE>
 293         >>> G = SymmetricLinearGame(L, K, e1, e2)
 294         >>> print(G)
 295         The linear game (L, K, e1, e2) where
 296           L = [ 1  2]
 297               [ 3  4],
 298           K = Nonnegative orthant in the real 2-space,
 299           e1 = [ 1]
 300                [ 1],
 301           e2 = [ 1]
 302                [ 1],
 303           Condition((L, K, e1, e2)) = 6.073...
 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
 326
 327     def __str__(self):
 328         """
 329         Return a string representation of this game.
 330         """
 331         tpl = 'The linear game (L, K, e1, e2) where\n' \
 332               '  L = {:s},\n' \
 333               '  K = {!s},\n' \
 334               '  e1 = {:s},\n' \
 335               '  e2 = {:s},\n' \
 336               '  Condition((L, K, e1, e2)) = {:f}.'
 337         indented_L = '\n      '.join(str(self._L).splitlines())
 338         indented_e1 = '\n       '.join(str(self._e1).splitlines())
 339         indented_e2 = '\n       '.join(str(self._e2).splitlines())
 340
 341         return tpl.format(indented_L,
 342                           str(self._K),
 343                           indented_e1,
 344                           indented_e2,
 345                           self.condition())
 346
 347
 348     def _zero(self):
 349         """
 350         Return a column of zeros that fits ``K``.
 351
 352         This is used in our CVXOPT construction.
 353
 354         .. warning::
 355
 356             It is not safe to cache any of the matrices passed to
 357             CVXOPT, because it can clobber them.
 358
 359         Returns
 360         -------
 361
 362         matrix
 363             A ``K.dimension()``-by-``1`` column vector of zeros.
 364
 365         Examples
 366         --------
 367
 368             >>> from dunshire import *
 369             >>> K = NonnegativeOrthant(3)
 370             >>> L = identity(3)
 371             >>> e1 = [1,1,1]
 372             >>> e2 = e1
 373             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
 374             >>> print(SLG._zero())
 375             [0.0000000]
 376             [0.0000000]
 377             [0.0000000]
 378             <BLANKLINE>
 379
 380         """
 381         return matrix(0, (self._K.dimension(), 1), tc='d')
 382
 383
 384     def _A(self):
 385         """
 386         Return the matrix ``A`` used in our CVXOPT construction.
 387
 388         This matrix ``A``  appears on the right-hand side of ``Ax = b``
 389         in the statement of the CVXOPT conelp program.
 390
 391         .. warning::
 392
 393             It is not safe to cache any of the matrices passed to
 394             CVXOPT, because it can clobber them.
 395
 396         Returns
 397         -------
 398
 399         matrix
 400             A ``1``-by-``(1 + K.dimension())`` row vector. Its first
 401             entry is zero, and the rest are the entries of ``e2``.
 402
 403         Examples
 404         --------
 405
 406             >>> from dunshire import *
 407             >>> K = NonnegativeOrthant(3)
 408             >>> L = [[1,1,1],[1,1,1],[1,1,1]]
 409             >>> e1 = [1,1,1]
 410             >>> e2 = [1,2,3]
 411             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
 412             >>> print(SLG._A())
 413             [0.0000000 1.0000000 2.0000000 3.0000000]
 414             <BLANKLINE>
 415
 416         """
 417         return matrix([0, self._e2], (1, self._K.dimension() + 1), 'd')
 418
 419
 420
 421     def _G(self):
 422         r"""
 423         Return the matrix ``G`` used in our CVXOPT construction.
 424
 425         Thus matrix ``G``that appears on the left-hand side of ``Gx + s = h``
 426         in the statement of the CVXOPT conelp program.
 427
 428         .. warning::
 429
 430             It is not safe to cache any of the matrices passed to
 431             CVXOPT, because it can clobber them.
 432
 433         Returns
 434         -------
 435
 436         matrix
 437             A ``2*K.dimension()``-by-``1 + K.dimension()`` matrix.
 438
 439         Examples
 440         --------
 441
 442             >>> from dunshire import *
 443             >>> K = NonnegativeOrthant(3)
 444             >>> L = [[4,5,6],[7,8,9],[10,11,12]]
 445             >>> e1 = [1,2,3]
 446             >>> e2 = [1,1,1]
 447             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
 448             >>> print(SLG._G())
 449             [  0.0000000  -1.0000000   0.0000000   0.0000000]
 450             [  0.0000000   0.0000000  -1.0000000   0.0000000]
 451             [  0.0000000   0.0000000   0.0000000  -1.0000000]
 452             [  1.0000000  -4.0000000  -5.0000000  -6.0000000]
 453             [  2.0000000  -7.0000000  -8.0000000  -9.0000000]
 454             [  3.0000000 -10.0000000 -11.0000000 -12.0000000]
 455             <BLANKLINE>
 456
 457         """
 458         I = identity(self._K.dimension())
 459         return append_row(append_col(self._zero(), -I),
 460                           append_col(self._e1, -self._L))
 461
 462
 463     def _try_solution(self, c, h, C, b, tolerance):
 464         # Actually solve the thing and obtain a dictionary describing
 465         # what happened.
 466         try:
 467             solvers.options['show_progress'] = options.VERBOSE
 468             solvers.options['abs_tol'] = tolerance
 469             soln_dict = solvers.conelp(c,self._G(),h,C,self._A(),b)
 470         except ValueError as e:
 471             if str(e) == 'math domain error':
 472                 # Oops, CVXOPT tried to take the square root of a
 473                 # negative number. Report some details about the game
 474                 # rather than just the underlying CVXOPT crash.
 475                 raise PoorScalingException(self)
 476             else:
 477                 raise e
 478
 479         # The optimal strategies are named ``p`` and ``q`` in the
 480         # background documentation, and we need to extract them from
 481         # the CVXOPT ``x`` and ``z`` variables. The objective values
 482         # :math:`nu` and :math:`omega` can also be found in the CVXOPT
 483         # ``x`` and ``y`` variables; however, they're stored
 484         # conveniently as separate entries in the solution dictionary.
 485         p1_value = -soln_dict['primal objective']
 486         p2_value = -soln_dict['dual objective']
 487         p1_optimal = soln_dict['x'][1:]
 488         p2_optimal = soln_dict['z'][self._K.dimension():]
 489
 490         # The "status" field contains "optimal" if everything went
 491         # according to plan. Other possible values are "primal
 492         # infeasible", "dual infeasible", "unknown", all of which mean
 493         # we didn't get a solution. The "infeasible" ones are the
 494         # worst, since they indicate that CVXOPT is convinced the
 495         # problem is infeasible (and that cannot happen).
 496         if soln_dict['status'] in ['primal infeasible', 'dual infeasible']:
 497             raise GameUnsolvableException(self, soln_dict)
 498         elif soln_dict['status'] == 'unknown':
 499             # When we get a status of "unknown", we may still be able
 500             # to salvage a solution out of the returned
 501             # dictionary. Often this is the result of numerical
 502             # difficulty and we can simply check that the primal/dual
 503             # objectives match (within a tolerance) and that the
 504             # primal/dual optimal solutions are within the cone (to a
 505             # tolerance as well).
 506             #
 507             # The fudge factor of two is basically unjustified, but
 508             # makes intuitive sense when you imagine that the primal
 509             # value could be under the true optimal by ``ABS_TOL``
 510             # and the dual value could be over by the same amount.
 511             #
 512             if abs(p1_value - p2_value) > tolerance:
 513                 raise GameUnsolvableException(self, soln_dict)
 514             if (p1_optimal not in self._K) or (p2_optimal not in self._K):
 515                 raise GameUnsolvableException(self, soln_dict)
 516
 517         return Solution(p1_value, p1_optimal, p2_optimal)
 518
 519
 520     def solution(self):
 521         """
 522         Solve this linear game and return a :class:`Solution`.
 523
 524         Returns
 525         -------
 526
 527         :class:`Solution`
 528             A :class:`Solution` object describing the game's value and
 529             the optimal strategies of both players.
 530
 531         Raises
 532         ------
 533         GameUnsolvableException
 534             If the game could not be solved (if an optimal solution to its
 535             associated cone program was not found).
 536
 537         PoorScalingException
 538             If the game could not be solved because CVXOPT crashed while
 539             trying to take the square root of a negative number.
 540
 541         Examples
 542         --------
 543
 544         This example is computed in Gowda and Ravindran in the section
 545         "The value of a Z-transformation"::
 546
 547             >>> from dunshire import *
 548             >>> K = NonnegativeOrthant(3)
 549             >>> L = [[1,-5,-15],[-1,2,-3],[-12,-15,1]]
 550             >>> e1 = [1,1,1]
 551             >>> e2 = [1,1,1]
 552             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
 553             >>> print(SLG.solution())
 554             Game value: -6.1724138
 555             Player 1 optimal:
 556               [ 0.551...]
 557               [-0.000...]
 558               [ 0.448...]
 559             Player 2 optimal:
 560               [0.448...]
 561               [0.000...]
 562               [0.551...]
 563
 564         The value of the following game can be computed using the fact
 565         that the identity is invertible::
 566
 567             >>> from dunshire import *
 568             >>> K = NonnegativeOrthant(3)
 569             >>> L = [[1,0,0],[0,1,0],[0,0,1]]
 570             >>> e1 = [1,2,3]
 571             >>> e2 = [4,5,6]
 572             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
 573             >>> print(SLG.solution())
 574             Game value: 0.0312500
 575             Player 1 optimal:
 576               [0.031...]
 577               [0.062...]
 578               [0.093...]
 579             Player 2 optimal:
 580               [0.125...]
 581               [0.156...]
 582               [0.187...]
 583
 584         """
 585         # The cone "C" that appears in the statement of the CVXOPT
 586         # conelp program.
 587         C = CartesianProduct(self._K, self._K)
 588
 589         # The column vector "b" that appears on the right-hand side of
 590         # Ax = b in the statement of the CVXOPT conelp program.
 591         b = matrix([1], tc='d')
 592
 593         # The column vector "h" that appears on the right-hand side of
 594         # Gx + s = h in the statement of the CVXOPT conelp program.
 595         h = matrix([self._zero(), self._zero()])
 596
 597         # The column vector "c" that appears in the objective function
 598         # value <c,x> in the statement of the CVXOPT conelp program.
 599         c = matrix([-1, self._zero()])
 600
 601         try:
 602             # First try with a stricter tolerance. Who knows, it might work.
 603             return self._try_solution(c, h, C.cvxopt_dims(), b,
 604                                       tolerance = options.ABS_TOL / 10)
 605
 606         except (PoorScalingException, GameUnsolvableException):
 607             # Ok, that didn't work. Let's try it with the default.
 608             return self._try_solution(c, h, C.cvxopt_dims(), b,
 609                                       tolerance = options.ABS_TOL)
 610
 611
 612     def condition(self):
 613         r"""
 614         Return the condition number of this game.
 615
 616         In the CVXOPT construction of this game, two matrices ``G`` and
 617         ``A`` appear. When those matrices are nasty, numerical problems
 618         can show up. We define the condition number of this game to be
 619         the average of the condition numbers of ``G`` and ``A`` in the
 620         CVXOPT construction. If the condition number of this game is
 621         high, then you can expect numerical difficulty (such as
 622         :class:`PoorScalingException`).
 623
 624         Returns
 625         -------
 626
 627         float
 628             A real number greater than or equal to one that measures how
 629             bad this game is numerically.
 630
 631         Examples
 632         --------
 633
 634         >>> from dunshire import *
 635         >>> K = NonnegativeOrthant(1)
 636         >>> L = [[1]]
 637         >>> e1 = [1]
 638         >>> e2 = e1
 639         >>> SLG = SymmetricLinearGame(L, K, e1, e2)
 640         >>> actual = SLG.condition()
 641         >>> expected = 1.8090169943749477
 642         >>> abs(actual - expected) < options.ABS_TOL
 643         True
 644
 645         """
 646         return (condition_number(self._G()) + condition_number(self._A()))/2
 647
 648
 649     def dual(self):
 650         r"""
 651         Return the dual game to this game.
 652
 653         If :math:`G = \left(L,K,e_{1},e_{2}\right)` is a linear game,
 654         then its dual is :math:`G^{*} =
 655         \left(L^{*},K^{*},e_{2},e_{1}\right)`. However, since this cone
 656         is symmetric, :math:`K^{*} = K`.
 657
 658         Examples
 659         --------
 660
 661             >>> from dunshire import *
 662             >>> K = NonnegativeOrthant(3)
 663             >>> L = [[1,-5,-15],[-1,2,-3],[-12,-15,1]]
 664             >>> e1 = [1,1,1]
 665             >>> e2 = [1,2,3]
 666             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
 667             >>> print(SLG.dual())
 668             The linear game (L, K, e1, e2) where
 669               L = [  1  -1 -12]
 670                   [ -5   2 -15]
 671                   [-15  -3   1],
 672               K = Nonnegative orthant in the real 3-space,
 673               e1 = [ 1]
 674                    [ 2]
 675                    [ 3],
 676               e2 = [ 1]
 677                    [ 1]
 678                    [ 1],
 679               Condition((L, K, e1, e2)) = 44.476...
 680
 681         """
 682         # We pass ``self._L`` right back into the constructor, because
 683         # it will be transposed there. And keep in mind that ``self._K``
 684         # is its own dual.
 685         return SymmetricLinearGame(self._L,
 686                                    self._K,
 687                                    self._e2,
 688                                    self._e1)