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 solvers.options['show_progress'] = options.VERBOSE
  16
  17
  18 class Solution:
  19     """
  20     A representation of the solution of a linear game. It should contain
  21     the value of the game, and both players' strategies.
  22
  23     Examples
  24     --------
  25
  26         >>> print(Solution(10, matrix([1,2]), matrix([3,4])))
  27         Game value: 10.0000000
  28         Player 1 optimal:
  29           [ 1]
  30           [ 2]
  31         Player 2 optimal:
  32           [ 3]
  33           [ 4]
  34
  35     """
  36     def __init__(self, game_value, p1_optimal, p2_optimal):
  37         """
  38         Create a new Solution object from a game value and two optimal
  39         strategies for the players.
  40         """
  41         self._game_value = game_value
  42         self._player1_optimal = p1_optimal
  43         self._player2_optimal = p2_optimal
  44
  45     def __str__(self):
  46         """
  47         Return a string describing the solution of a linear game.
  48
  49         The three data that are described are,
  50
  51           * The value of the game.
  52           * The optimal strategy of player one.
  53           * The optimal strategy of player two.
  54
  55         The two optimal strategy vectors are indented by two spaces.
  56         """
  57         tpl = 'Game value: {:.7f}\n' \
  58               'Player 1 optimal:{:s}\n' \
  59               'Player 2 optimal:{:s}'
  60
  61         p1_str = '\n{!s}'.format(self.player1_optimal())
  62         p1_str = '\n  '.join(p1_str.splitlines())
  63         p2_str = '\n{!s}'.format(self.player2_optimal())
  64         p2_str = '\n  '.join(p2_str.splitlines())
  65
  66         return tpl.format(self.game_value(), p1_str, p2_str)
  67
  68
  69     def game_value(self):
  70         """
  71         Return the game value for this solution.
  72
  73         Examples
  74         --------
  75
  76             >>> s = Solution(10, matrix([1,2]), matrix([3,4]))
  77             >>> s.game_value()
  78             10
  79
  80         """
  81         return self._game_value
  82
  83
  84     def player1_optimal(self):
  85         """
  86         Return player one's optimal strategy in this solution.
  87
  88         Examples
  89         --------
  90
  91             >>> s = Solution(10, matrix([1,2]), matrix([3,4]))
  92             >>> print(s.player1_optimal())
  93             [ 1]
  94             [ 2]
  95             <BLANKLINE>
  96
  97         """
  98         return self._player1_optimal
  99
 100
 101     def player2_optimal(self):
 102         """
 103         Return player two's optimal strategy in this solution.
 104
 105         Examples
 106         --------
 107
 108             >>> s = Solution(10, matrix([1,2]), matrix([3,4]))
 109             >>> print(s.player2_optimal())
 110             [ 3]
 111             [ 4]
 112             <BLANKLINE>
 113
 114         """
 115         return self._player2_optimal
 116
 117
 118 class SymmetricLinearGame:
 119     r"""
 120     A representation of a symmetric linear game.
 121
 122     The data for a symmetric linear game are,
 123
 124       * A "payoff" operator ``L``.
 125       * A symmetric cone ``K``.
 126       * Two points ``e1`` and ``e2`` in the interior of ``K``.
 127
 128     The ambient space is assumed to be the span of ``K``.
 129
 130     With those data understood, the game is played as follows. Players
 131     one and two choose points :math:`x` and :math:`y` respectively, from
 132     their respective strategy sets,
 133
 134     .. math::
 135         \begin{aligned}
 136             \Delta_{1}
 137             &=
 138             \left\{
 139               x \in K \ \middle|\ \left\langle x, e_{2} \right\rangle = 1
 140             \right\}\\
 141             \Delta_{2}
 142             &=
 143             \left\{
 144               y \in K \ \middle|\ \left\langle y, e_{1} \right\rangle = 1
 145             \right\}.
 146         \end{aligned}
 147
 148     Afterwards, a "payout" is computed as :math:`\left\langle
 149     L\left(x\right), y \right\rangle` and is paid to player one out of
 150     player two's pocket. The game is therefore zero sum, and we suppose
 151     that player one would like to guarantee himself the largest minimum
 152     payout possible. That is, player one wishes to,
 153
 154     .. math::
 155         \begin{aligned}
 156             \text{maximize }
 157               &\underset{y \in \Delta_{2}}{\min}\left(
 158                 \left\langle L\left(x\right), y \right\rangle
 159               \right)\\
 160             \text{subject to } & x \in \Delta_{1}.
 161         \end{aligned}
 162
 163     Player two has the simultaneous goal to,
 164
 165     .. math::
 166         \begin{aligned}
 167             \text{minimize }
 168               &\underset{x \in \Delta_{1}}{\max}\left(
 169                 \left\langle L\left(x\right), y \right\rangle
 170               \right)\\
 171             \text{subject to } & y \in \Delta_{2}.
 172         \end{aligned}
 173
 174     These goals obviously conflict (the game is zero sum), but an
 175     existence theorem guarantees at least one optimal min-max solution
 176     from which neither player would like to deviate. This class is
 177     able to find such a solution.
 178
 179     Parameters
 180     ----------
 181
 182     L : list of list of float
 183         A matrix represented as a list of ROWS. This representation
 184         agrees with (for example) SageMath and NumPy, but not with CVXOPT
 185         (whose matrix constructor accepts a list of columns).
 186
 187     K : :class:`SymmetricCone`
 188         The symmetric cone instance over which the game is played.
 189
 190     e1 : iterable float
 191         The interior point of ``K`` belonging to player one; it
 192         can be of any iterable type having the correct length.
 193
 194     e2 : iterable float
 195         The interior point of ``K`` belonging to player two; it
 196         can be of any enumerable type having the correct length.
 197
 198     Raises
 199     ------
 200
 201     ValueError
 202         If either ``e1`` or ``e2`` lie outside of the cone ``K``.
 203
 204     Examples
 205     --------
 206
 207         >>> from dunshire import *
 208         >>> K = NonnegativeOrthant(3)
 209         >>> L = [[1,-5,-15],[-1,2,-3],[-12,-15,1]]
 210         >>> e1 = [1,1,1]
 211         >>> e2 = [1,2,3]
 212         >>> SLG = SymmetricLinearGame(L, K, e1, e2)
 213         >>> print(SLG)
 214         The linear game (L, K, e1, e2) where
 215           L = [  1  -5 -15]
 216               [ -1   2  -3]
 217               [-12 -15   1],
 218           K = Nonnegative orthant in the real 3-space,
 219           e1 = [ 1]
 220                [ 1]
 221                [ 1],
 222           e2 = [ 1]
 223                [ 2]
 224                [ 3],
 225           Condition((L, K, e1, e2)) = 31.834...
 226
 227     Lists can (and probably should) be used for every argument::
 228
 229         >>> from dunshire import *
 230         >>> K = NonnegativeOrthant(2)
 231         >>> L = [[1,0],[0,1]]
 232         >>> e1 = [1,1]
 233         >>> e2 = [1,1]
 234         >>> G = SymmetricLinearGame(L, K, e1, e2)
 235         >>> print(G)
 236         The linear game (L, K, e1, e2) where
 237           L = [ 1  0]
 238               [ 0  1],
 239           K = Nonnegative orthant in the real 2-space,
 240           e1 = [ 1]
 241                [ 1],
 242           e2 = [ 1]
 243                [ 1],
 244           Condition((L, K, e1, e2)) = 1.707...
 245
 246     The points ``e1`` and ``e2`` can also be passed as some other
 247     enumerable type (of the correct length) without much harm, since
 248     there is no row/column ambiguity::
 249
 250         >>> import cvxopt
 251         >>> import numpy
 252         >>> from dunshire import *
 253         >>> K = NonnegativeOrthant(2)
 254         >>> L = [[1,0],[0,1]]
 255         >>> e1 = cvxopt.matrix([1,1])
 256         >>> e2 = numpy.matrix([1,1])
 257         >>> G = SymmetricLinearGame(L, K, e1, e2)
 258         >>> print(G)
 259         The linear game (L, K, e1, e2) where
 260           L = [ 1  0]
 261               [ 0  1],
 262           K = Nonnegative orthant in the real 2-space,
 263           e1 = [ 1]
 264                [ 1],
 265           e2 = [ 1]
 266                [ 1],
 267           Condition((L, K, e1, e2)) = 1.707...
 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           Condition((L, K, e1, e2)) = 6.073...
 290         >>> L = cvxopt.matrix(L)
 291         >>> print(L)
 292         [ 1  3]
 293         [ 2  4]
 294         <BLANKLINE>
 295         >>> G = SymmetricLinearGame(L, K, e1, e2)
 296         >>> print(G)
 297         The linear game (L, K, e1, e2) where
 298           L = [ 1  2]
 299               [ 3  4],
 300           K = Nonnegative orthant in the real 2-space,
 301           e1 = [ 1]
 302                [ 1],
 303           e2 = [ 1]
 304                [ 1],
 305           Condition((L, K, e1, e2)) = 6.073...
 306
 307     """
 308     def __init__(self, L, K, e1, e2):
 309         """
 310         Create a new SymmetricLinearGame object.
 311         """
 312         self._K = K
 313         self._e1 = matrix(e1, (K.dimension(), 1))
 314         self._e2 = matrix(e2, (K.dimension(), 1))
 315
 316         # Our input ``L`` is indexed by rows but CVXOPT matrices are
 317         # indexed by columns, so we need to transpose the input before
 318         # feeding it to CVXOPT.
 319         self._L = matrix(L, (K.dimension(), K.dimension())).trans()
 320
 321         if not self._e1 in K:
 322             raise ValueError('the point e1 must lie in the interior of K')
 323
 324         if not self._e2 in K:
 325             raise ValueError('the point e2 must lie in the interior of K')
 326
 327         # Cached result of the self._zero() method.
 328         self._zero_col = None
 329
 330
 331     def __str__(self):
 332         """
 333         Return a string representation of this game.
 334         """
 335         tpl = 'The linear game (L, K, e1, e2) where\n' \
 336               '  L = {:s},\n' \
 337               '  K = {!s},\n' \
 338               '  e1 = {:s},\n' \
 339               '  e2 = {:s},\n' \
 340               '  Condition((L, K, e1, e2)) = {:f}.'
 341         indented_L = '\n      '.join(str(self._L).splitlines())
 342         indented_e1 = '\n       '.join(str(self._e1).splitlines())
 343         indented_e2 = '\n       '.join(str(self._e2).splitlines())
 344
 345         return tpl.format(indented_L,
 346                           str(self._K),
 347                           indented_e1,
 348                           indented_e2,
 349                           self.condition())
 350
 351
 352     def _zero(self):
 353         """
 354         Return a column of zeros that fits ``K``.
 355
 356         This is used in our CVXOPT construction.
 357         """
 358         if self._zero_col is None:
 359             # Cache it, it's constant.
 360             self._zero_col = matrix(0, (self._K.dimension(), 1), tc='d')
 361         return self._zero_col
 362
 363
 364     def _A(self):
 365         """
 366         Return the matrix ``A`` used in our CVXOPT construction.
 367
 368         This matrix ``A``  appears on the right-hand side of ``Ax = b``
 369         in the statement of the CVXOPT conelp program.
 370         """
 371         return matrix([0, self._e2], (1, self._K.dimension() + 1), 'd')
 372
 373
 374     def _G(self):
 375         r"""
 376         Return the matrix ``G`` used in our CVXOPT construction.
 377
 378         Thus matrix ``G``that appears on the left-hand side of ``Gx + s = h``
 379         in the statement of the CVXOPT conelp program.
 380         """
 381         I = identity(self._K.dimension())
 382         return append_row(append_col(self._zero(), -I),
 383                           append_col(self._e1, -self._L))
 384
 385
 386     def solution(self):
 387         """
 388         Solve this linear game and return a :class:`Solution`.
 389
 390         Returns
 391         -------
 392
 393         :class:`Solution`
 394             A :class:`Solution` object describing the game's value and
 395             the optimal strategies of both players.
 396
 397         Raises
 398         ------
 399         GameUnsolvableException
 400             If the game could not be solved (if an optimal solution to its
 401             associated cone program was not found).
 402
 403         PoorScalingException
 404             If the game could not be solved because CVXOPT crashed while
 405             trying to take the square root of a negative number.
 406
 407         Examples
 408         --------
 409
 410         This example is computed in Gowda and Ravindran in the section
 411         "The value of a Z-transformation"::
 412
 413             >>> from dunshire import *
 414             >>> K = NonnegativeOrthant(3)
 415             >>> L = [[1,-5,-15],[-1,2,-3],[-12,-15,1]]
 416             >>> e1 = [1,1,1]
 417             >>> e2 = [1,1,1]
 418             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
 419             >>> print(SLG.solution())
 420             Game value: -6.1724138
 421             Player 1 optimal:
 422               [ 0.551...]
 423               [-0.000...]
 424               [ 0.448...]
 425             Player 2 optimal:
 426               [0.448...]
 427               [0.000...]
 428               [0.551...]
 429
 430         The value of the following game can be computed using the fact
 431         that the identity is invertible::
 432
 433             >>> from dunshire import *
 434             >>> K = NonnegativeOrthant(3)
 435             >>> L = [[1,0,0],[0,1,0],[0,0,1]]
 436             >>> e1 = [1,2,3]
 437             >>> e2 = [4,5,6]
 438             >>> SLG = SymmetricLinearGame(L, K, e1, e2)
 439             >>> print(SLG.solution())
 440             Game value: 0.0312500
 441             Player 1 optimal:
 442               [0.031...]
 443               [0.062...]
 444               [0.093...]
 445             Player 2 optimal:
 446               [0.125...]
 447               [0.156...]
 448               [0.187...]
 449
 450         """
 451         # The cone "C" that appears in the statement of the CVXOPT
 452         # conelp program.
 453         C = CartesianProduct(self._K, self._K)
 454
 455         # The column vector "b" that appears on the right-hand side of
 456         # Ax = b in the statement of the CVXOPT conelp program.
 457         b = matrix([1], tc='d')
 458
 459         # The column vector "h" that appears on the right-hand side of
 460         # Gx + s = h in the statement of the CVXOPT conelp program.
 461         h = matrix([self._zero(), self._zero()])
 462
 463         # The column vector "c" that appears in the objective function
 464         # value <c,x> in the statement of the CVXOPT conelp program.
 465         c = matrix([-1, self._zero()])
 466
 467         # Actually solve the thing and obtain a dictionary describing
 468         # what happened.
 469         try:
 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) > 2*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)