src/test/symmetric_linear_game_test.py

   1 # These few are used only for tests.
   2 from math import sqrt
   3 from random import randint, uniform
   4 from unittest import TestCase
   5
   6 from cvxopt import matrix
   7 from dunshire.cones import NonnegativeOrthant, IceCream
   8 from dunshire.games import SymmetricLinearGame
   9 from dunshire.matrices import (append_col, append_row, eigenvalues_re,
  10                                identity, inner_product)
  11 from dunshire import options
  12
  13
  14 def random_matrix(dims):
  15     """
  16     Generate a random square matrix.
  17
  18     Parameters
  19     ----------
  20
  21     dims : int
  22         The number of rows/columns you want in the returned matrix.
  23
  24     Returns
  25     -------
  26
  27     matrix
  28         A new matrix whose entries are random floats chosen uniformly from
  29         the interval [-10, 10].
  30
  31     Examples
  32     --------
  33
  34         >>> A = random_matrix(3)
  35         >>> A.size
  36         (3, 3)
  37
  38     """
  39     return matrix([[uniform(-10, 10) for i in range(dims)]
  40                    for j in range(dims)])
  41
  42
  43 def random_nonnegative_matrix(dims):
  44     """
  45     Generate a random square matrix with nonnegative entries.
  46
  47     Parameters
  48     ----------
  49
  50     dims : int
  51         The number of rows/columns you want in the returned matrix.
  52
  53     Returns
  54     -------
  55
  56     matrix
  57         A new matrix whose entries are random floats chosen uniformly from
  58         the interval [0, 10].
  59
  60     Examples
  61     --------
  62
  63         >>> A = random_nonnegative_matrix(3)
  64         >>> A.size
  65         (3, 3)
  66         >>> all([entry >= 0 for entry in A])
  67         True
  68
  69     """
  70     L = random_matrix(dims)
  71     return matrix([abs(entry) for entry in L], (dims, dims))
  72
  73
  74 def random_diagonal_matrix(dims):
  75     """
  76     Generate a random square matrix with zero off-diagonal entries.
  77
  78     These matrices are Lyapunov-like on the nonnegative orthant, as is
  79     fairly easy to see.
  80
  81     Parameters
  82     ----------
  83
  84     dims : int
  85         The number of rows/columns you want in the returned matrix.
  86
  87     Returns
  88     -------
  89
  90     matrix
  91         A new matrix whose diagonal entries are random floats chosen
  92         uniformly from the interval [-10, 10] and whose off-diagonal
  93         entries are zero.
  94
  95     Examples
  96     --------
  97
  98         >>> A = random_diagonal_matrix(3)
  99         >>> A.size
 100         (3, 3)
 101         >>> A[0,1] == A[0,2] == A[1,0] == A[2,0] == A[1,2] == A[2,1] == 0
 102         True
 103
 104     """
 105     return matrix([[uniform(-10, 10)*int(i == j) for i in range(dims)]
 106                    for j in range(dims)])
 107
 108
 109 def random_skew_symmetric_matrix(dims):
 110     """
 111     Generate a random skew-symmetrix matrix.
 112
 113     Parameters
 114     ----------
 115
 116     dims : int
 117         The number of rows/columns you want in the returned matrix.
 118
 119     Returns
 120     -------
 121
 122     matrix
 123         A new skew-matrix whose strictly above-diagonal entries are
 124         random floats chosen uniformly from the interval [-10, 10].
 125
 126     Examples
 127     --------
 128
 129         >>> A = random_skew_symmetric_matrix(3)
 130         >>> A.size
 131         (3, 3)
 132
 133         >>> from dunshire.matrices import norm
 134         >>> A = random_skew_symmetric_matrix(randint(1, 10))
 135         >>> norm(A + A.trans()) < options.ABS_TOL
 136         True
 137
 138     """
 139     strict_ut = [[uniform(-10, 10)*int(i < j) for i in range(dims)]
 140                  for j in range(dims)]
 141
 142     strict_ut = matrix(strict_ut, (dims, dims))
 143     return strict_ut - strict_ut.trans()
 144
 145
 146 def random_lyapunov_like_icecream(dims):
 147     r"""
 148     Generate a random matrix Lyapunov-like on the ice-cream cone.
 149
 150     The form of these matrices is cited in Gowda and Tao
 151     [GowdaTao]_. The scalar ``a`` and the vector ``b`` (using their
 152     notation) are easy to generate. The submatrix ``D`` is a little
 153     trickier, but it can be found noticing that :math:`C + C^{T} = 0`
 154     for a skew-symmetric matrix :math:`C` implying that :math:`C + C^{T}
 155     + \left(2a\right)I = \left(2a\right)I`. Thus we can stick an
 156     :math:`aI` with each of :math:`C,C^{T}` and let those be our
 157     :math:`D,D^{T}`.
 158
 159     Parameters
 160     ----------
 161
 162     dims : int
 163         The dimension of the ice-cream cone (not of the matrix you want!)
 164         on which the returned matrix should be Lyapunov-like.
 165
 166     Returns
 167     -------
 168
 169     matrix
 170         A new matrix, Lyapunov-like on the ice-cream cone in ``dims``
 171         dimensions, whose free entries are random floats chosen uniformly
 172         from the interval [-10, 10].
 173
 174     References
 175     ----------
 176
 177     .. [GowdaTao] M. S. Gowda and J. Tao. On the bilinearity rank of a
 178        proper cone and Lyapunov-like transformations. Mathematical
 179        Programming, 147:155–170, 2014.
 180
 181     Examples
 182     --------
 183
 184         >>> L = random_lyapunov_like_icecream(3)
 185         >>> L.size
 186         (3, 3)
 187         >>> x = matrix([1,1,0])
 188         >>> s = matrix([1,-1,0])
 189         >>> abs(inner_product(L*x, s)) < options.ABS_TOL
 190         True
 191
 192     """
 193     a = matrix([uniform(-10, 10)], (1, 1))
 194     b = matrix([uniform(-10, 10) for idx in range(dims-1)], (dims-1, 1))
 195     D = random_skew_symmetric_matrix(dims-1) + a*identity(dims-1)
 196     row1 = append_col(a, b.trans())
 197     row2 = append_col(b, D)
 198     return append_row(row1, row2)
 199
 200
 201 def random_orthant_params():
 202     """
 203     Generate the ``L``, ``K``, ``e1``, and ``e2`` parameters for a
 204     random game over the nonnegative orthant.
 205     """
 206     ambient_dim = randint(1, 10)
 207     K = NonnegativeOrthant(ambient_dim)
 208     e1 = [uniform(0.5, 10) for idx in range(K.dimension())]
 209     e2 = [uniform(0.5, 10) for idx in range(K.dimension())]
 210     L = random_matrix(K.dimension())
 211     return (L, K, matrix(e1), matrix(e2))
 212
 213
 214 def random_icecream_params():
 215     """
 216     Generate the ``L``, ``K``, ``e1``, and ``e2`` parameters for a
 217     random game over the ice-cream cone.
 218     """
 219     # Use a minimum dimension of two to avoid divide-by-zero in
 220     # the fudge factor we make up later.
 221     ambient_dim = randint(2, 10)
 222     K = IceCream(ambient_dim)
 223     e1 = [1] # Set the "height" of e1 to one
 224     e2 = [1] # And the same for e2
 225
 226     # If we choose the rest of the components of e1,e2 randomly
 227     # between 0 and 1, then the largest the squared norm of the
 228     # non-height part of e1,e2 could be is the 1*(dim(K) - 1). We
 229     # need to make it less than one (the height of the cone) so
 230     # that the whole thing is in the cone. The norm of the
 231     # non-height part is sqrt(dim(K) - 1), and we can divide by
 232     # twice that.
 233     fudge_factor = 1.0 / (2.0*sqrt(K.dimension() - 1.0))
 234     e1 += [fudge_factor*uniform(0, 1) for idx in range(K.dimension() - 1)]
 235     e2 += [fudge_factor*uniform(0, 1) for idx in range(K.dimension() - 1)]
 236     L = random_matrix(K.dimension())
 237
 238     return (L, K, matrix(e1), matrix(e2))
 239
 240
 241 # Tell pylint to shut up about the large number of methods.
 242 class SymmetricLinearGameTest(TestCase): # pylint: disable=R0904
 243     """
 244     Tests for the SymmetricLinearGame and Solution classes.
 245     """
 246     def assert_within_tol(self, first, second):
 247         """
 248         Test that ``first`` and ``second`` are equal within our default
 249         tolerance.
 250         """
 251         self.assertTrue(abs(first - second) < options.ABS_TOL)
 252
 253
 254     def assert_norm_within_tol(self, first, second):
 255         """
 256         Test that ``first`` and ``second`` vectors are equal in the
 257         sense that the norm of their difference is within our default
 258         tolerance.
 259         """
 260         self.assert_within_tol(norm(first - second), 0)
 261
 262
 263     def assert_solution_exists(self, L, K, e1, e2):
 264         """
 265         Given the parameters needed to construct a SymmetricLinearGame,
 266         ensure that that game has a solution.
 267         """
 268         # The matrix() constructor assumes that ``L`` is a list of
 269         # columns, so we transpose it to agree with what
 270         # SymmetricLinearGame() thinks.
 271         G = SymmetricLinearGame(L.trans(), K, e1, e2)
 272         soln = G.solution()
 273
 274         expected = inner_product(L*soln.player1_optimal(),
 275                                  soln.player2_optimal())
 276         self.assert_within_tol(soln.game_value(), expected)
 277
 278
 279     def test_solution_exists_orthant(self):
 280         """
 281         Every linear game has a solution, so we should be able to solve
 282         every symmetric linear game over the NonnegativeOrthant. Pick
 283         some parameters randomly and give it a shot. The resulting
 284         optimal solutions should give us the optimal game value when we
 285         apply the payoff operator to them.
 286         """
 287         (L, K, e1, e2) = random_orthant_params()
 288         self.assert_solution_exists(L, K, e1, e2)
 289
 290
 291     def test_solution_exists_icecream(self):
 292         """
 293         Like :meth:`test_solution_exists_nonnegative_orthant`, except
 294         over the ice cream cone.
 295         """
 296         (L, K, e1, e2) = random_icecream_params()
 297         self.assert_solution_exists(L, K, e1, e2)
 298
 299
 300     def test_negative_value_z_operator(self):
 301         """
 302         Test the example given in Gowda/Ravindran of a Z-matrix with
 303         negative game value on the nonnegative orthant.
 304         """
 305         K = NonnegativeOrthant(2)
 306         e1 = [1, 1]
 307         e2 = e1
 308         L = [[1, -2], [-2, 1]]
 309         G = SymmetricLinearGame(L, K, e1, e2)
 310         self.assertTrue(G.solution().game_value() < -options.ABS_TOL)
 311
 312
 313     def assert_scaling_works(self, L, K, e1, e2):
 314         """
 315         Test that scaling ``L`` by a nonnegative number scales the value
 316         of the game by the same number.
 317         """
 318         game1 = SymmetricLinearGame(L, K, e1, e2)
 319         value1 = game1.solution().game_value()
 320
 321         alpha = uniform(0.1, 10)
 322         game2 = SymmetricLinearGame(alpha*L, K, e1, e2)
 323         value2 = game2.solution().game_value()
 324         self.assert_within_tol(alpha*value1, value2)
 325
 326
 327     def test_scaling_orthant(self):
 328         """
 329         Test that scaling ``L`` by a nonnegative number scales the value
 330         of the game by the same number over the nonnegative orthant.
 331         """
 332         (L, K, e1, e2) = random_orthant_params()
 333         self.assert_scaling_works(L, K, e1, e2)
 334
 335
 336     def test_scaling_icecream(self):
 337         """
 338         The same test as :meth:`test_nonnegative_scaling_orthant`,
 339         except over the ice cream cone.
 340         """
 341         (L, K, e1, e2) = random_icecream_params()
 342         self.assert_scaling_works(L, K, e1, e2)
 343
 344
 345     def assert_translation_works(self, L, K, e1, e2):
 346         """
 347         Check that translating ``L`` by alpha*(e1*e2.trans()) increases
 348         the value of the associated game by alpha.
 349         """
 350         # We need to use ``L`` later, so make sure we transpose it
 351         # before passing it in as a column-indexed matrix.
 352         game1 = SymmetricLinearGame(L.trans(), K, e1, e2)
 353         soln1 = game1.solution()
 354         value1 = soln1.game_value()
 355         x_bar = soln1.player1_optimal()
 356         y_bar = soln1.player2_optimal()
 357
 358         alpha = uniform(-10, 10)
 359         tensor_prod = e1*e2.trans()
 360
 361         # This is the "correct" representation of ``M``, but COLUMN
 362         # indexed...
 363         M = L + alpha*tensor_prod
 364
 365         # so we have to transpose it when we feed it to the constructor.
 366         game2 = SymmetricLinearGame(M.trans(), K, e1, e2)
 367         value2 = game2.solution().game_value()
 368
 369         self.assert_within_tol(value1 + alpha, value2)
 370
 371         # Make sure the same optimal pair works.
 372         self.assert_within_tol(value2, inner_product(M*x_bar, y_bar))
 373
 374
 375     def test_translation_orthant(self):
 376         """
 377         Test that translation works over the nonnegative orthant.
 378         """
 379         (L, K, e1, e2) = random_orthant_params()
 380         self.assert_translation_works(L, K, e1, e2)
 381
 382
 383     def test_translation_icecream(self):
 384         """
 385         The same as :meth:`test_translation_orthant`, except over the
 386         ice cream cone.
 387         """
 388         (L, K, e1, e2) = random_icecream_params()
 389         self.assert_translation_works(L, K, e1, e2)
 390
 391
 392     def assert_opposite_game_works(self, L, K, e1, e2):
 393         """
 394         Check the value of the "opposite" game that gives rise to a
 395         value that is the negation of the original game. Comes from
 396         some corollary.
 397         """
 398         # We need to use ``L`` later, so make sure we transpose it
 399         # before passing it in as a column-indexed matrix.
 400         game1 = SymmetricLinearGame(L.trans(), K, e1, e2)
 401
 402         # This is the "correct" representation of ``M``, but
 403         # COLUMN indexed...
 404         M = -L.trans()
 405
 406         # so we have to transpose it when we feed it to the constructor.
 407         game2 = SymmetricLinearGame(M.trans(), K, e2, e1)
 408
 409         soln1 = game1.solution()
 410         x_bar = soln1.player1_optimal()
 411         y_bar = soln1.player2_optimal()
 412         soln2 = game2.solution()
 413
 414         self.assert_within_tol(-soln1.game_value(), soln2.game_value())
 415
 416         # Make sure the switched optimal pair works.
 417         self.assert_within_tol(soln2.game_value(),
 418                                inner_product(M*y_bar, x_bar))
 419
 420
 421     def test_opposite_game_orthant(self):
 422         """
 423         Test the value of the "opposite" game over the nonnegative
 424         orthant.
 425         """
 426         (L, K, e1, e2) = random_orthant_params()
 427         self.assert_opposite_game_works(L, K, e1, e2)
 428
 429
 430     def test_opposite_game_icecream(self):
 431         """
 432         Like :meth:`test_opposite_game_orthant`, except over the
 433         ice-cream cone.
 434         """
 435         (L, K, e1, e2) = random_icecream_params()
 436         self.assert_opposite_game_works(L, K, e1, e2)
 437
 438
 439     def assert_orthogonality(self, L, K, e1, e2):
 440         """
 441         Two orthogonality relations hold at an optimal solution, and we
 442         check them here.
 443         """
 444         # We need to use ``L`` later, so make sure we transpose it
 445         # before passing it in as a column-indexed matrix.
 446         game = SymmetricLinearGame(L.trans(), K, e1, e2)
 447         soln = game.solution()
 448         x_bar = soln.player1_optimal()
 449         y_bar = soln.player2_optimal()
 450         value = soln.game_value()
 451
 452         ip1 = inner_product(y_bar, L*x_bar - value*e1)
 453         self.assert_within_tol(ip1, 0)
 454
 455         ip2 = inner_product(value*e2 - L.trans()*y_bar, x_bar)
 456         self.assert_within_tol(ip2, 0)
 457
 458
 459     def test_orthogonality_orthant(self):
 460         """
 461         Check the orthgonality relationships that hold for a solution
 462         over the nonnegative orthant.
 463         """
 464         (L, K, e1, e2) = random_orthant_params()
 465         self.assert_orthogonality(L, K, e1, e2)
 466
 467
 468     def test_orthogonality_icecream(self):
 469         """
 470         Check the orthgonality relationships that hold for a solution
 471         over the ice-cream cone.
 472         """
 473         (L, K, e1, e2) = random_icecream_params()
 474         self.assert_orthogonality(L, K, e1, e2)
 475
 476
 477     def test_positive_operator_value(self):
 478         """
 479         Test that a positive operator on the nonnegative orthant gives
 480         rise to a a game with a nonnegative value.
 481
 482         This test theoretically applies to the ice-cream cone as well,
 483         but we don't know how to make positive operators on that cone.
 484         """
 485         (K, e1, e2) = random_orthant_params()[1:]
 486         L = random_nonnegative_matrix(K.dimension())
 487
 488         game = SymmetricLinearGame(L, K, e1, e2)
 489         self.assertTrue(game.solution().game_value() >= -options.ABS_TOL)
 490
 491
 492     def assert_lyapunov_works(self, L, K, e1, e2):
 493         """
 494         Check that Lyapunov games act the way we expect.
 495         """
 496         game = SymmetricLinearGame(L, K, e1, e2)
 497         soln = game.solution()
 498
 499         # We only check for positive/negative stability if the game
 500         # value is not basically zero. If the value is that close to
 501         # zero, we just won't check any assertions.
 502         eigs = eigenvalues_re(L)
 503         if soln.game_value() > options.ABS_TOL:
 504             # L should be positive stable
 505             positive_stable = all([eig > -options.ABS_TOL for eig in eigs])
 506             self.assertTrue(positive_stable)
 507         elif soln.game_value() < -options.ABS_TOL:
 508             # L should be negative stable
 509             negative_stable = all([eig < options.ABS_TOL for eig in eigs])
 510             self.assertTrue(negative_stable)
 511
 512         # The dual game's value should always equal the primal's.
 513         dualsoln = game.dual().solution()
 514         self.assert_within_tol(dualsoln.game_value(), soln.game_value())
 515
 516
 517     def test_lyapunov_orthant(self):
 518         """
 519         Test that a Lyapunov game on the nonnegative orthant works.
 520         """
 521         (K, e1, e2) = random_orthant_params()[1:]
 522         L = random_diagonal_matrix(K.dimension())
 523
 524         self.assert_lyapunov_works(L, K, e1, e2)
 525
 526
 527     def test_lyapunov_icecream(self):
 528         """
 529         Test that a Lyapunov game on the ice-cream cone works.
 530         """
 531         (K, e1, e2) = random_icecream_params()[1:]
 532         L = random_lyapunov_like_icecream(K.dimension())
 533
 534         self.assert_lyapunov_works(L, K, e1, e2)