X-Git-Url: https://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=src%2Fdunshire%2Fgames.py;h=3d4b09ad8c0c12ff1bba0294a1c5c87ae8ff72bf;hb=bd5a4b0c8519de2a938663acdbc080560f829628;hp=e2d29788f2b2612fe0d23e12be07010902c93451;hpb=ff576f2eddae5554e4888d99b092d73e19619a91;p=dunshire.git diff --git a/src/dunshire/games.py b/src/dunshire/games.py index e2d2978..3d4b09a 100644 --- a/src/dunshire/games.py +++ b/src/dunshire/games.py @@ -14,7 +14,8 @@ from unittest import TestCase from cvxopt import matrix, printing, solvers from cones import CartesianProduct, IceCream, NonnegativeOrthant from errors import GameUnsolvableException -from matrices import append_col, append_row, identity, inner_product +from matrices import (append_col, append_row, eigenvalues_re, identity, + inner_product, norm) import options printing.options['dformat'] = options.FLOAT_FORMAT @@ -122,23 +123,91 @@ class Solution: class SymmetricLinearGame: - """ + r""" A representation of a symmetric linear game. - The data for a linear game are, + The data for a symmetric linear game are, * A "payoff" operator ``L``. - * A cone ``K``. - * A point ``e`` in the interior of ``K``. - * A point ``f`` in the interior of the dual of ``K``. - - In a symmetric game, the cone ``K`` is be self-dual. We therefore - name the two interior points ``e1`` and ``e2`` to indicate that - they come from the same cone but are "chosen" by players one and - two respectively. + * A symmetric cone ``K``. + * Two points ``e1`` and ``e2`` in the interior of ``K``. The ambient space is assumed to be the span of ``K``. + With those data understood, the game is played as follows. Players + one and two choose points :math:`x` and :math:`y` respectively, from + their respective strategy sets, + + .. math:: + \begin{aligned} + \Delta_{1} + &= + \left\{ + x \in K \ \middle|\ \left\langle x, e_{2} \right\rangle = 1 + \right\}\\ + \Delta_{2} + &= + \left\{ + y \in K \ \middle|\ \left\langle y, e_{1} \right\rangle = 1 + \right\}. + \end{aligned} + + Afterwards, a "payout" is computed as :math:`\left\langle + L\left(x\right), y \right\rangle` and is paid to player one out of + player two's pocket. The game is therefore zero sum, and we suppose + that player one would like to guarantee himself the largest minimum + payout possible. That is, player one wishes to, + + .. math:: + \begin{aligned} + \text{maximize } + &\underset{y \in \Delta_{2}}{\min}\left( + \left\langle L\left(x\right), y \right\rangle + \right)\\ + \text{subject to } & x \in \Delta_{1}. + \end{aligned} + + Player two has the simultaneous goal to, + + .. math:: + \begin{aligned} + \text{minimize } + &\underset{x \in \Delta_{1}}{\max}\left( + \left\langle L\left(x\right), y \right\rangle + \right)\\ + \text{subject to } & y \in \Delta_{2}. + \end{aligned} + + These goals obviously conflict (the game is zero sum), but an + existence theorem guarantees at least one optimal min-max solution + from which neither player would like to deviate. This class is + able to find such a solution. + + Parameters + ---------- + + L : list of list of float + A matrix represented as a list of ROWS. This representation + agrees with (for example) SageMath and NumPy, but not with CVXOPT + (whose matrix constructor accepts a list of columns). + + K : :class:`SymmetricCone` + The symmetric cone instance over which the game is played. + + e1 : iterable float + The interior point of ``K`` belonging to player one; it + can be of any iterable type having the correct length. + + e2 : iterable float + The interior point of ``K`` belonging to player two; it + can be of any enumerable type having the correct length. + + Raises + ------ + + ValueError + If either ``e1`` or ``e2`` lie outside of the cone ``K``. + Examples -------- @@ -161,7 +230,6 @@ class SymmetricLinearGame: [ 2] [ 3]. - Lists can (and probably should) be used for every argument:: >>> from cones import NonnegativeOrthant @@ -242,23 +310,6 @@ class SymmetricLinearGame: def __init__(self, L, K, e1, e2): """ Create a new SymmetricLinearGame object. - - INPUT: - - - ``L`` -- an square matrix represented as a list of lists - of real numbers. ``L`` itself is interpreted as a list of - ROWS, which agrees with (for example) SageMath and NumPy, - but not with CVXOPT (whose matrix constructor accepts a - list of columns). - - - ``K`` -- a SymmetricCone instance. - - - ``e1`` -- the interior point of ``K`` belonging to player one; - it can be of any enumerable type having the correct length. - - - ``e2`` -- the interior point of ``K`` belonging to player two; - it can be of any enumerable type having the correct length. - """ self._K = K self._e1 = matrix(e1, (K.dimension(), 1)) @@ -269,10 +320,10 @@ class SymmetricLinearGame: # feeding it to CVXOPT. self._L = matrix(L, (K.dimension(), K.dimension())).trans() - if not K.contains_strict(self._e1): + if not self._e1 in K: raise ValueError('the point e1 must lie in the interior of K') - if not K.contains_strict(self._e2): + if not self._e2 in K: raise ValueError('the point e2 must lie in the interior of K') def __str__(self): @@ -292,16 +343,20 @@ class SymmetricLinearGame: def solution(self): """ - Solve this linear game and return a Solution object. + Solve this linear game and return a :class:`Solution`. - OUTPUT: + Returns + ------- - If the cone program associated with this game could be - successfully solved, then a Solution object containing the - game's value and optimal strategies is returned. If the game - could *not* be solved -- which should never happen -- then a - GameUnsolvableException is raised. It can be printed to get the - raw output from CVXOPT. + :class:`Solution` + A :class:`Solution` object describing the game's value and + the optimal strategies of both players. + + Raises + ------ + GameUnsolvableException + If the game could not be solved (if an optimal solution to its + associated cone program was not found). Examples -------- @@ -379,25 +434,44 @@ class SymmetricLinearGame: # what happened. soln_dict = solvers.conelp(c, G, h, C.cvxopt_dims(), A, b) + p1_value = -soln_dict['primal objective'] + p2_value = -soln_dict['dual objective'] + p1_optimal = soln_dict['x'][1:] + p2_optimal = soln_dict['z'][self._K.dimension():] + # The "status" field contains "optimal" if everything went # according to plan. Other possible values are "primal - # infeasible", "dual infeasible", "unknown", all of which - # mean we didn't get a solution. That should never happen, - # because by construction our game has a solution, and thus - # the cone program should too. - if soln_dict['status'] != 'optimal': + # infeasible", "dual infeasible", "unknown", all of which mean + # we didn't get a solution. The "infeasible" ones are the + # worst, since they indicate that CVXOPT is convinced the + # problem is infeasible (and that cannot happen). + if soln_dict['status'] in ['primal infeasible', 'dual infeasible']: raise GameUnsolvableException(soln_dict) - - p1_value = soln_dict['x'][0] - p1_optimal = soln_dict['x'][1:] - p2_optimal = soln_dict['z'][self._K.dimension():] + elif soln_dict['status'] == 'unknown': + # When we get a status of "unknown", we may still be able + # to salvage a solution out of the returned + # dictionary. Often this is the result of numerical + # difficulty and we can simply check that the primal/dual + # objectives match (within a tolerance) and that the + # primal/dual optimal solutions are within the cone (to a + # tolerance as well). + if abs(p1_value - p2_value) > options.ABS_TOL: + raise GameUnsolvableException(soln_dict) + if (p1_optimal not in self._K) or (p2_optimal not in self._K): + raise GameUnsolvableException(soln_dict) return Solution(p1_value, p1_optimal, p2_optimal) + def dual(self): - """ + r""" Return the dual game to this game. + If :math:`G = \left(L,K,e_{1},e_{2}\right)` is a linear game, + then its dual is :math:`G^{*} = + \left(L^{*},K^{*},e_{2},e_{1}\right)`. However, since this cone + is symmetric, :math:`K^{*} = K`. + Examples -------- @@ -430,11 +504,112 @@ class SymmetricLinearGame: self._e1) + +def _random_matrix(dims): + """ + Generate a random square (``dims``-by-``dims``) matrix. This is used + only by the :class:`SymmetricLinearGameTest` class. + """ + return matrix([[uniform(-10, 10) for i in range(dims)] + for j in range(dims)]) + +def _random_nonnegative_matrix(dims): + """ + Generate a random square (``dims``-by-``dims``) matrix with + nonnegative entries. This is used only by the + :class:`SymmetricLinearGameTest` class. + """ + L = _random_matrix(dims) + return matrix([abs(entry) for entry in L], (dims, dims)) + +def _random_diagonal_matrix(dims): + """ + Generate a random square (``dims``-by-``dims``) matrix with nonzero + entries only on the diagonal. This is used only by the + :class:`SymmetricLinearGameTest` class. + """ + return matrix([[uniform(-10, 10)*int(i == j) for i in range(dims)] + for j in range(dims)]) + + +def _random_skew_symmetric_matrix(dims): + """ + Generate a random skew-symmetrix (``dims``-by-``dims``) matrix. + + Examples + -------- + + >>> A = _random_skew_symmetric_matrix(randint(1, 10)) + >>> norm(A + A.trans()) < options.ABS_TOL + True + + """ + strict_ut = [[uniform(-10, 10)*int(i < j) for i in range(dims)] + for j in range(dims)] + + strict_ut = matrix(strict_ut, (dims,dims)) + return (strict_ut - strict_ut.trans()) + + +def _random_lyapunov_like_icecream(dims): + """ + Generate a random Lyapunov-like matrix over the ice-cream cone in + ``dims`` dimensions. + """ + a = matrix([uniform(-10,10)], (1,1)) + b = matrix([uniform(-10,10) for idx in range(dims-1)], (dims-1, 1)) + D = _random_skew_symmetric_matrix(dims-1) + a*identity(dims-1) + row1 = append_col(a, b.trans()) + row2 = append_col(b, D) + return append_row(row1,row2) + + +def _random_orthant_params(): + """ + Generate the ``L``, ``K``, ``e1``, and ``e2`` parameters for a + random game over the nonnegative orthant. This is only used by + the :class:`SymmetricLinearGameTest` class. + """ + ambient_dim = randint(1, 10) + K = NonnegativeOrthant(ambient_dim) + e1 = [uniform(0.5, 10) for idx in range(K.dimension())] + e2 = [uniform(0.5, 10) for idx in range(K.dimension())] + L = _random_matrix(K.dimension()) + return (L, K, matrix(e1), matrix(e2)) + + +def _random_icecream_params(): + """ + Generate the ``L``, ``K``, ``e1``, and ``e2`` parameters for a + random game over the ice cream cone. This is only used by + the :class:`SymmetricLinearGameTest` class. + """ + # Use a minimum dimension of two to avoid divide-by-zero in + # the fudge factor we make up later. + ambient_dim = randint(2, 10) + K = IceCream(ambient_dim) + e1 = [1] # Set the "height" of e1 to one + e2 = [1] # And the same for e2 + + # If we choose the rest of the components of e1,e2 randomly + # between 0 and 1, then the largest the squared norm of the + # non-height part of e1,e2 could be is the 1*(dim(K) - 1). We + # need to make it less than one (the height of the cone) so + # that the whole thing is in the cone. The norm of the + # non-height part is sqrt(dim(K) - 1), and we can divide by + # twice that. + fudge_factor = 1.0 / (2.0*sqrt(K.dimension() - 1.0)) + e1 += [fudge_factor*uniform(0, 1) for idx in range(K.dimension() - 1)] + e2 += [fudge_factor*uniform(0, 1) for idx in range(K.dimension() - 1)] + L = _random_matrix(K.dimension()) + + return (L, K, matrix(e1), matrix(e2)) + + class SymmetricLinearGameTest(TestCase): """ Tests for the SymmetricLinearGame and Solution classes. """ - def assert_within_tol(self, first, second): """ Test that ``first`` and ``second`` are equal within our default @@ -443,19 +618,32 @@ class SymmetricLinearGameTest(TestCase): self.assertTrue(abs(first - second) < options.ABS_TOL) + def assert_norm_within_tol(self, first, second): + """ + Test that ``first`` and ``second`` vectors are equal in the + sense that the norm of their difference is within our default + tolerance. + """ + self.assert_within_tol(norm(first - second), 0) + + def assert_solution_exists(self, L, K, e1, e2): """ Given the parameters needed to construct a SymmetricLinearGame, ensure that that game has a solution. """ - G = SymmetricLinearGame(L, K, e1, e2) + # The matrix() constructor assumes that ``L`` is a list of + # columns, so we transpose it to agree with what + # SymmetricLinearGame() thinks. + G = SymmetricLinearGame(L.trans(), K, e1, e2) soln = G.solution() - L_matrix = matrix(L).trans() - expected = inner_product(L_matrix*soln.player1_optimal(), + + expected = inner_product(L*soln.player1_optimal(), soln.player2_optimal()) self.assert_within_tol(soln.game_value(), expected) - def test_solution_exists_nonnegative_orthant(self): + + def test_solution_exists_orthant(self): """ Every linear game has a solution, so we should be able to solve every symmetric linear game over the NonnegativeOrthant. Pick @@ -463,37 +651,258 @@ class SymmetricLinearGameTest(TestCase): optimal solutions should give us the optimal game value when we apply the payoff operator to them. """ - ambient_dim = randint(1, 10) - K = NonnegativeOrthant(ambient_dim) - e1 = [uniform(0.1, 10) for idx in range(K.dimension())] - e2 = [uniform(0.1, 10) for idx in range(K.dimension())] - L = [[uniform(-10, 10) for i in range(K.dimension())] - for j in range(K.dimension())] + (L, K, e1, e2) = _random_orthant_params() self.assert_solution_exists(L, K, e1, e2) - def test_solution_exists_ice_cream(self): + + def test_solution_exists_icecream(self): """ Like :meth:`test_solution_exists_nonnegative_orthant`, except over the ice cream cone. """ - # Use a minimum dimension of two to avoid divide-by-zero in - # the fudge factor we make up later. - ambient_dim = randint(2, 10) - K = IceCream(ambient_dim) - e1 = [1] # Set the "height" of e1 to one - e2 = [1] # And the same for e2 - - # If we choose the rest of the components of e1,e2 randomly - # between 0 and 1, then the largest the squared norm of the - # non-height part of e1,e2 could be is the 1*(dim(K) - 1). We - # need to make it less than one (the height of the cone) so - # that the whole thing is in the cone. The norm of the - # non-height part is sqrt(dim(K) - 1), and we can divide by - # twice that. - fudge_factor = 1.0 / (2.0*sqrt(K.dimension() - 1.0)) - e1 += [fudge_factor*uniform(0, 1) for idx in range(K.dimension() - 1)] - e2 += [fudge_factor*uniform(0, 1) for idx in range(K.dimension() - 1)] - L = [[uniform(-10, 10) for i in range(K.dimension())] - for j in range(K.dimension())] + (L, K, e1, e2) = _random_icecream_params() self.assert_solution_exists(L, K, e1, e2) + + def test_negative_value_z_operator(self): + """ + Test the example given in Gowda/Ravindran of a Z-matrix with + negative game value on the nonnegative orthant. + """ + K = NonnegativeOrthant(2) + e1 = [1, 1] + e2 = e1 + L = [[1, -2], [-2, 1]] + G = SymmetricLinearGame(L, K, e1, e2) + self.assertTrue(G.solution().game_value() < -options.ABS_TOL) + + + def assert_scaling_works(self, L, K, e1, e2): + """ + Test that scaling ``L`` by a nonnegative number scales the value + of the game by the same number. + """ + game1 = SymmetricLinearGame(L, K, e1, e2) + value1 = game1.solution().game_value() + + alpha = uniform(0.1, 10) + game2 = SymmetricLinearGame(alpha*L, K, e1, e2) + value2 = game2.solution().game_value() + self.assert_within_tol(alpha*value1, value2) + + + def test_scaling_orthant(self): + """ + Test that scaling ``L`` by a nonnegative number scales the value + of the game by the same number over the nonnegative orthant. + """ + (L, K, e1, e2) = _random_orthant_params() + self.assert_scaling_works(L, K, e1, e2) + + + def test_scaling_icecream(self): + """ + The same test as :meth:`test_nonnegative_scaling_orthant`, + except over the ice cream cone. + """ + (L, K, e1, e2) = _random_icecream_params() + self.assert_scaling_works(L, K, e1, e2) + + + def assert_translation_works(self, L, K, e1, e2): + """ + Check that translating ``L`` by alpha*(e1*e2.trans()) increases + the value of the associated game by alpha. + """ + # We need to use ``L`` later, so make sure we transpose it + # before passing it in as a column-indexed matrix. + game1 = SymmetricLinearGame(L.trans(), K, e1, e2) + soln1 = game1.solution() + value1 = soln1.game_value() + x_bar = soln1.player1_optimal() + y_bar = soln1.player2_optimal() + + alpha = uniform(-10, 10) + tensor_prod = e1*e2.trans() + + # This is the "correct" representation of ``M``, but COLUMN + # indexed... + M = L + alpha*tensor_prod + + # so we have to transpose it when we feed it to the constructor. + game2 = SymmetricLinearGame(M.trans(), K, e1, e2) + value2 = game2.solution().game_value() + + self.assert_within_tol(value1 + alpha, value2) + + # Make sure the same optimal pair works. + self.assert_within_tol(value2, inner_product(M*x_bar, y_bar)) + + + def test_translation_orthant(self): + """ + Test that translation works over the nonnegative orthant. + """ + (L, K, e1, e2) = _random_orthant_params() + self.assert_translation_works(L, K, e1, e2) + + + def test_translation_icecream(self): + """ + The same as :meth:`test_translation_orthant`, except over the + ice cream cone. + """ + (L, K, e1, e2) = _random_icecream_params() + self.assert_translation_works(L, K, e1, e2) + + + def assert_opposite_game_works(self, L, K, e1, e2): + """ + Check the value of the "opposite" game that gives rise to a + value that is the negation of the original game. Comes from + some corollary. + """ + # We need to use ``L`` later, so make sure we transpose it + # before passing it in as a column-indexed matrix. + game1 = SymmetricLinearGame(L.trans(), K, e1, e2) + + # This is the "correct" representation of ``M``, but + # COLUMN indexed... + M = -L.trans() + + # so we have to transpose it when we feed it to the constructor. + game2 = SymmetricLinearGame(M.trans(), K, e2, e1) + + soln1 = game1.solution() + x_bar = soln1.player1_optimal() + y_bar = soln1.player2_optimal() + soln2 = game2.solution() + + self.assert_within_tol(-soln1.game_value(), soln2.game_value()) + + # Make sure the switched optimal pair works. + self.assert_within_tol(soln2.game_value(), + inner_product(M*y_bar, x_bar)) + + + def test_opposite_game_orthant(self): + """ + Test the value of the "opposite" game over the nonnegative + orthant. + """ + (L, K, e1, e2) = _random_orthant_params() + self.assert_opposite_game_works(L, K, e1, e2) + + + def test_opposite_game_icecream(self): + """ + Like :meth:`test_opposite_game_orthant`, except over the + ice-cream cone. + """ + (L, K, e1, e2) = _random_icecream_params() + self.assert_opposite_game_works(L, K, e1, e2) + + + def assert_orthogonality(self, L, K, e1, e2): + """ + Two orthogonality relations hold at an optimal solution, and we + check them here. + """ + # We need to use ``L`` later, so make sure we transpose it + # before passing it in as a column-indexed matrix. + game = SymmetricLinearGame(L.trans(), K, e1, e2) + soln = game.solution() + x_bar = soln.player1_optimal() + y_bar = soln.player2_optimal() + value = soln.game_value() + + ip1 = inner_product(y_bar, L*x_bar - value*e1) + self.assert_within_tol(ip1, 0) + + ip2 = inner_product(value*e2 - L.trans()*y_bar, x_bar) + self.assert_within_tol(ip2, 0) + + + def test_orthogonality_orthant(self): + """ + Check the orthgonality relationships that hold for a solution + over the nonnegative orthant. + """ + (L, K, e1, e2) = _random_orthant_params() + self.assert_orthogonality(L, K, e1, e2) + + + def test_orthogonality_icecream(self): + """ + Check the orthgonality relationships that hold for a solution + over the ice-cream cone. + """ + (L, K, e1, e2) = _random_icecream_params() + self.assert_orthogonality(L, K, e1, e2) + + + def test_positive_operator_value(self): + """ + Test that a positive operator on the nonnegative orthant gives + rise to a a game with a nonnegative value. + + This test theoretically applies to the ice-cream cone as well, + but we don't know how to make positive operators on that cone. + """ + (_, K, e1, e2) = _random_orthant_params() + + # Ignore that L, we need a nonnegative one. + L = _random_nonnegative_matrix(K.dimension()) + + game = SymmetricLinearGame(L, K, e1, e2) + self.assertTrue(game.solution().game_value() >= -options.ABS_TOL) + + + def assert_lyapunov_works(self, L, K, e1, e2): + """ + Check that Lyapunov games act the way we expect. + """ + game = SymmetricLinearGame(L, K, e1, e2) + soln = game.solution() + + # We only check for positive/negative stability if the game + # value is not basically zero. If the value is that close to + # zero, we just won't check any assertions. + if soln.game_value() > options.ABS_TOL: + # L should be positive stable + ps = all([eig > -options.ABS_TOL for eig in eigenvalues_re(L)]) + self.assertTrue(ps) + elif soln.game_value() < -options.ABS_TOL: + # L should be negative stable + ns = all([eig < options.ABS_TOL for eig in eigenvalues_re(L)]) + self.assertTrue(ns) + + # The dual game's value should always equal the primal's. + dualsoln = game.dual().solution() + self.assert_within_tol(dualsoln.game_value(), soln.game_value()) + + + def test_lyapunov_orthant(self): + """ + Test that a Lyapunov game on the nonnegative orthant works. + """ + (L, K, e1, e2) = _random_orthant_params() + + # Ignore that L, we need a diagonal (Lyapunov-like) one. + # (And we don't need to transpose those.) + L = _random_diagonal_matrix(K.dimension()) + + self.assert_lyapunov_works(L, K, e1, e2) + + + def test_lyapunov_icecream(self): + """ + Test that a Lyapunov game on the ice-cream cone works. + """ + (L, K, e1, e2) = _random_icecream_params() + + # Ignore that L, we need a diagonal (Lyapunov-like) one. + # (And we don't need to transpose those.) + L = _random_lyapunov_like_icecream(K.dimension()) + + self.assert_lyapunov_works(L, K, e1, e2)