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         # Cached result of the self._zero() method.
 326         self._zero_col = 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},\n' \
 338               '  Condition((L, K, e1, e2)) = {:f}.'
 339         indented_L = '\n      '.join(str(self._L).splitlines())
 340         indented_e1 = '\n       '.join(str(self._e1).splitlines())
 341         indented_e2 = '\n       '.join(str(self._e2).splitlines())
 342
 343         return tpl.format(indented_L,
 344                           str(self._K),
 345                           indented_e1,
 346                           indented_e2,
 347                           self.condition())
 348
 349
 350     def _zero(self):
 351         """
 352         Return a column of zeros that fits ``K``.
 353
 354         This is used in our CVXOPT construction.
 355         """
 356         if self._zero_col is None:
 357             # Cache it, it's constant.
 358             self._zero_col = matrix(0, (self._K.dimension(), 1), tc='d')
 359         return self._zero_col
 360
 361
 362     def _A(self):
 363         """
 364         Return the matrix ``A`` used in our CVXOPT construction.
 365
 366         This matrix ``A``  appears on the right-hand side of ``Ax = b``
 367         in the statement of the CVXOPT conelp program.
 368         """
 369         return matrix([0, self._e2], (1, self._K.dimension() + 1), 'd')
 370
 371
 372     def _G(self):
 373         r"""
 374         Return the matrix ``G`` used in our CVXOPT construction.
 375
 376         Thus matrix ``G``that appears on the left-hand side of ``Gx + s = h``
 377         in the statement of the CVXOPT conelp program.
 378         """
 379         I = identity(self._K.dimension())
 380         return append_row(append_col(self._zero(), -I),
 381                           append_col(self._e1, -self._L))
 382
 383
 384     def solution(self):
 385         """
 386         Solve this linear game and return a :class:`Solution`.
 387
 388         Returns
 389         -------
 390
 391         :class:`Solution`
 392             A :class:`Solution` object describing the game's value and
 393             the optimal strategies of both players.
 394
 395         Raises
 396         ------
 397         GameUnsolvableException
 398             If the game could not be solved (if an optimal solution to its
 399             associated cone program was not found).
 400
 401         PoorScalingException
 402             If the game could not be solved because CVXOPT crashed while
 403             trying to take the square root of a negative number.
 404
 405         Examples
 406         --------
 407
 408         This example is computed in Gowda and Ravindran in the section
 409         "The value of a Z-transformation"::
 410
 411             >>> from dunshire import *
 412             >>> K = NonnegativeOrthant(3)
 413             >>> L = [[1,-5,-15],[-1,2,-3],[-12,-15,1]]
 414             >>> e1 = [1,1,1]
 415             >>> e2 = [1,1,1]
 416             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
 417             >>> print(SLG.solution())
 418             Game value: -6.1724138
 419             Player 1 optimal:
 420               [ 0.551...]
 421               [-0.000...]
 422               [ 0.448...]
 423             Player 2 optimal:
 424               [0.448...]
 425               [0.000...]
 426               [0.551...]
 427
 428         The value of the following game can be computed using the fact
 429         that the identity is invertible::
 430
 431             >>> from dunshire import *
 432             >>> K = NonnegativeOrthant(3)
 433             >>> L = [[1,0,0],[0,1,0],[0,0,1]]
 434             >>> e1 = [1,2,3]
 435             >>> e2 = [4,5,6]
 436             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
 437             >>> print(SLG.solution())
 438             Game value: 0.0312500
 439             Player 1 optimal:
 440               [0.031...]
 441               [0.062...]
 442               [0.093...]
 443             Player 2 optimal:
 444               [0.125...]
 445               [0.156...]
 446               [0.187...]
 447
 448         """
 449         # The cone "C" that appears in the statement of the CVXOPT
 450         # conelp program.
 451         C = CartesianProduct(self._K, self._K)
 452
 453         # The column vector "b" that appears on the right-hand side of
 454         # Ax = b in the statement of the CVXOPT conelp program.
 455         b = matrix([1], tc='d')
 456
 457         # The column vector "h" that appears on the right-hand side of
 458         # Gx + s = h in the statement of the CVXOPT conelp program.
 459         h = matrix([self._zero(), self._zero()])
 460
 461         # The column vector "c" that appears in the objective function
 462         # value <c,x> in the statement of the CVXOPT conelp program.
 463         c = matrix([-1, self._zero()])
 464
 465         # Actually solve the thing and obtain a dictionary describing
 466         # what happened.
 467         try:
 468             solvers.options['show_progress'] = options.VERBOSE
 469             solvers.options['abs_tol'] = options.ABS_TOL
 470             soln_dict = solvers.conelp(c, self._G(), h,
 471                                        C.cvxopt_dims(), self._A(), b)
 472         except ValueError as e:
 473             if str(e) == 'math domain error':
 474                 # Oops, CVXOPT tried to take the square root of a
 475                 # negative number. Report some details about the game
 476                 # rather than just the underlying CVXOPT crash.
 477                 raise PoorScalingException(self)
 478             else:
 479                 raise e
 480
 481         # The optimal strategies are named ``p`` and ``q`` in the
 482         # background documentation, and we need to extract them from
 483         # the CVXOPT ``x`` and ``z`` variables. The objective values
 484         # :math:`nu` and :math:`omega` can also be found in the CVXOPT
 485         # ``x`` and ``y`` variables; however, they're stored
 486         # conveniently as separate entries in the solution dictionary.
 487         p1_value = -soln_dict['primal objective']
 488         p2_value = -soln_dict['dual objective']
 489         p1_optimal = soln_dict['x'][1:]
 490         p2_optimal = soln_dict['z'][self._K.dimension():]
 491
 492         # The "status" field contains "optimal" if everything went
 493         # according to plan. Other possible values are "primal
 494         # infeasible", "dual infeasible", "unknown", all of which mean
 495         # we didn't get a solution. The "infeasible" ones are the
 496         # worst, since they indicate that CVXOPT is convinced the
 497         # problem is infeasible (and that cannot happen).
 498         if soln_dict['status'] in ['primal infeasible', 'dual infeasible']:
 499             raise GameUnsolvableException(self, soln_dict)
 500         elif soln_dict['status'] == 'unknown':
 501             # When we get a status of "unknown", we may still be able
 502             # to salvage a solution out of the returned
 503             # dictionary. Often this is the result of numerical
 504             # difficulty and we can simply check that the primal/dual
 505             # objectives match (within a tolerance) and that the
 506             # primal/dual optimal solutions are within the cone (to a
 507             # tolerance as well).
 508             #
 509             # The fudge factor of two is basically unjustified, but
 510             # makes intuitive sense when you imagine that the primal
 511             # value could be under the true optimal by ``ABS_TOL``
 512             # and the dual value could be over by the same amount.
 513             #
 514             if abs(p1_value - p2_value) > options.ABS_TOL:
 515                 raise GameUnsolvableException(self, soln_dict)
 516             if (p1_optimal not in self._K) or (p2_optimal not in self._K):
 517                 raise GameUnsolvableException(self, soln_dict)
 518
 519         return Solution(p1_value, p1_optimal, p2_optimal)
 520
 521
 522     def condition(self):
 523         r"""
 524         Return the condition number of this game.
 525
 526         In the CVXOPT construction of this game, two matrices ``G`` and
 527         ``A`` appear. When those matrices are nasty, numerical problems
 528         can show up. We define the condition number of this game to be
 529         the average of the condition numbers of ``G`` and ``A`` in the
 530         CVXOPT construction. If the condition number of this game is
 531         high, then you can expect numerical difficulty (such as
 532         :class:`PoorScalingException`).
 533
 534         Returns
 535         -------
 536
 537         float
 538             A real number greater than or equal to one that measures how
 539             bad this game is numerically.
 540
 541         Examples
 542         --------
 543
 544         >>> from dunshire import *
 545         >>> K = NonnegativeOrthant(1)
 546         >>> L = [[1]]
 547         >>> e1 = [1]
 548         >>> e2 = e1
 549         >>> SLG = SymmetricLinearGame(L, K, e1, e2)
 550         >>> actual = SLG.condition()
 551         >>> expected = 1.8090169943749477
 552         >>> abs(actual - expected) < options.ABS_TOL
 553         True
 554
 555         """
 556         return (condition_number(self._G()) + condition_number(self._A()))/2
 557
 558
 559     def dual(self):
 560         r"""
 561         Return the dual game to this game.
 562
 563         If :math:`G = \left(L,K,e_{1},e_{2}\right)` is a linear game,
 564         then its dual is :math:`G^{*} =
 565         \left(L^{*},K^{*},e_{2},e_{1}\right)`. However, since this cone
 566         is symmetric, :math:`K^{*} = K`.
 567
 568         Examples
 569         --------
 570
 571             >>> from dunshire import *
 572             >>> K = NonnegativeOrthant(3)
 573             >>> L = [[1,-5,-15],[-1,2,-3],[-12,-15,1]]
 574             >>> e1 = [1,1,1]
 575             >>> e2 = [1,2,3]
 576             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
 577             >>> print(SLG.dual())
 578             The linear game (L, K, e1, e2) where
 579               L = [  1  -1 -12]
 580                   [ -5   2 -15]
 581                   [-15  -3   1],
 582               K = Nonnegative orthant in the real 3-space,
 583               e1 = [ 1]
 584                    [ 2]
 585                    [ 3],
 586               e2 = [ 1]
 587                    [ 1]
 588                    [ 1],
 589               Condition((L, K, e1, e2)) = 44.476...
 590
 591         """
 592         # We pass ``self._L`` right back into the constructor, because
 593         # it will be transposed there. And keep in mind that ``self._K``
 594         # is its own dual.
 595         return SymmetricLinearGame(self._L,
 596                                    self._K,
 597                                    self._e2,
 598                                    self._e1)