From 13b585b28aaacb1d603c3ae41614bacf613afa14 Mon Sep 17 00:00:00 2001 From: Michael Orlitzky Date: Wed, 9 Nov 2016 09:23:23 -0500 Subject: [PATCH] Replace _try_solution() with something more reliable and update tests. This huge change eliminates _try_solution() entirely. It turns out that if we just solve with the default tolerance and then *test* what we get back (regardless of the "status" field), we do better than the two-phase _try_solution() approach. Well, we fail less, anyway. The code for solution is now fairly simple. It solves the problem, and if it isn't infeasible, checks the solution for sanity: the primal/dual values are close (within 2*ABS_TOL) and the optimal points are in the cone. If those two things are true, we return the solution even if CVXOPT said "unknown". This fixes two test failures, which are now included as doctests to ensure that they can be solved. Moreover, the value of a game is now set to be the payoff at the optimal points. Before, we simply took the primal optimal value from CVXOPT. That was causing some test failures though, and either one is just as good as the other and mathematically unjustified as yet. There were existing tests to check the payoff at the optimal points, but they became redundant and were removed. Finally, all of the tests have been updated to use more conservative modifiers, not based on the condition number of the game. Some failures are still being ironed out, but they are rare. Note: the _c(), _h(), etc. methods on the game are now overkill since they are only used once, but they don't hurt I guess. --- dunshire/games.py | 273 +++++++++++++---------------- test/failures/001.py | 7 - test/failures/002.py | 7 - test/symmetric_linear_game_test.py | 106 +++-------- 4 files changed, 141 insertions(+), 252 deletions(-) delete mode 100644 test/failures/001.py delete mode 100644 test/failures/002.py diff --git a/dunshire/games.py b/dunshire/games.py index 2d6d6da..130176b 100644 --- a/dunshire/games.py +++ b/dunshire/games.py @@ -809,29 +809,10 @@ class SymmetricLinearGame: return matrix([1], tc='d') - def _try_solution(self, tolerance): - """ - Solve this linear game within ``tolerance``, if possible. - - This private function is the one that does all of the actual - work for :meth:`solution`. This method accepts a ``tolerance``, - and what :meth:`solution` does is call this method twice with - two different tolerances. First it tries a strict tolerance, and - then it tries a looser one. - - .. warning:: - - If you try to be smart and precompute the matrices used by - this function (the ones passed to ``conelp``), then you're - going to shoot yourself in the foot. CVXOPT can and will - clobber some (but not all) of its input matrices. This isn't - performance sensitive, so play it safe. - Parameters - ---------- - - tolerance : float - The absolute tolerance to pass to the CVXOPT solver. + def solution(self): + """ + Solve this linear game and return a :class:`Solution`. Returns ------- @@ -853,53 +834,102 @@ class SymmetricLinearGame: Examples -------- - This game can be solved easily, so the first attempt in - :meth:`solution` should succeed:: + This example is computed in Gowda and Ravindran in the section + "The value of a Z-transformation":: >>> from dunshire import * - >>> from dunshire.matrices import norm - >>> from dunshire.options import ABS_TOL >>> K = NonnegativeOrthant(3) >>> L = [[1,-5,-15],[-1,2,-3],[-12,-15,1]] >>> e1 = [1,1,1] >>> e2 = [1,1,1] >>> SLG = SymmetricLinearGame(L, K, e1, e2) - >>> s1 = SLG.solution() - >>> s2 = SLG._try_solution(options.ABS_TOL) - >>> abs(s1.game_value() - s2.game_value()) < ABS_TOL - True - >>> norm(s1.player1_optimal() - s2.player1_optimal()) < ABS_TOL - True - >>> norm(s1.player2_optimal() - s2.player2_optimal()) < ABS_TOL - True + >>> print(SLG.solution()) + Game value: -6.1724138 + Player 1 optimal: + [ 0.551...] + [-0.000...] + [ 0.448...] + Player 2 optimal: + [0.448...] + [0.000...] + [0.551...] - This game cannot be solved with the default tolerance, but it - can be solved with a weaker one:: + The value of the following game can be computed using the fact + that the identity is invertible:: >>> from dunshire import * - >>> from dunshire.options import ABS_TOL - >>> L = [[ 0.58538005706658102767, 1.53764301129883040886], - ... [-1.34901059721452210027, 1.50121179114155500756]] - >>> K = NonnegativeOrthant(2) - >>> e1 = [1.04537193228494995623, 1.39699624965841895374] - >>> e2 = [0.35326554172108337593, 0.11795703527854853321] - >>> SLG = SymmetricLinearGame(L,K,e1,e2) - >>> print(SLG._try_solution(ABS_TOL / 10)) - Traceback (most recent call last): - ... - dunshire.errors.GameUnsolvableException: Solution failed... - >>> print(SLG._try_solution(ABS_TOL)) - Game value: 9.1100945 + >>> K = NonnegativeOrthant(3) + >>> L = [[1,0,0],[0,1,0],[0,0,1]] + >>> e1 = [1,2,3] + >>> e2 = [4,5,6] + >>> SLG = SymmetricLinearGame(L, K, e1, e2) + >>> print(SLG.solution()) + Game value: 0.0312500 Player 1 optimal: - [-0.0000000] - [ 8.4776631] + [0.031...] + [0.062...] + [0.093...] Player 2 optimal: - [0.0000000] - [0.7158216] + [0.125...] + [0.156...] + [0.187...] + + This is another Gowda/Ravindran example that is supposed to have + a negative game value:: + + >>> from dunshire import * + >>> from dunshire.options import ABS_TOL + >>> L = [[1, -2], [-2, 1]] + >>> K = NonnegativeOrthant(2) + >>> e1 = [1, 1] + >>> e2 = e1 + >>> SLG = SymmetricLinearGame(L, K, e1, e2) + >>> SLG.solution().game_value() < -ABS_TOL + True + + Tests + ----- + + The following two games are problematic numerically, but we + should be able to solve them:: + + >>> from dunshire import * + >>> L = [[-0.95237953890954685221, 1.83474556206462535712], + ... [ 1.30481749924621448500, 1.65278664543326403447]] + >>> K = NonnegativeOrthant(2) + >>> e1 = [0.95477167524644313001, 0.63270781756540095397] + >>> e2 = [0.39633793037154141370, 0.10239281495640320530] + >>> SLG = SymmetricLinearGame(L, K, e1, e2) + >>> print(SLG.solution()) + Game value: 18.767... + Player 1 optimal: + [-0.000...] + [ 9.766...] + Player 2 optimal: + [1.047...] + [0.000...] + + :: + + >>> from dunshire import * + >>> L = [[1.54159395026049472754, 2.21344728574316684799], + ... [1.33147433507846657541, 1.17913616272988108769]] + >>> K = NonnegativeOrthant(2) + >>> e1 = [0.39903040089404784307, 0.12377403622479113410] + >>> e2 = [0.15695181142215544612, 0.85527381344651265405] + >>> SLG = SymmetricLinearGame(L, K, e1, e2) + >>> print(SLG.solution()) + Game value: 24.614... + Player 1 optimal: + [ 6.371...] + [-0.000...] + Player 2 optimal: + [2.506...] + [0.000...] """ try: - opts = {'show_progress': options.VERBOSE, 'abstol': tolerance} + opts = {'show_progress': options.VERBOSE} soln_dict = solvers.conelp(self._c(), self._G(), self._h(), @@ -930,112 +960,45 @@ class SymmetricLinearGame: # 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. The "infeasible" ones are the - # worst, since they indicate that CVXOPT is convinced the - # problem is infeasible (and that cannot happen). + # 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(self, soln_dict) - 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). - # - # The fudge factor of two is basically unjustified, but - # makes intuitive sense when you imagine that the primal - # value could be under the true optimal by ``ABS_TOL`` - # and the dual value could be over by the same amount. - # - if abs(p1_value - p2_value) > tolerance: - raise GameUnsolvableException(self, soln_dict) - if (p1_optimal not in self._K) or (p2_optimal not in self._K): - raise GameUnsolvableException(self, soln_dict) - - return Solution(p1_value, p1_optimal, p2_optimal) - - def solution(self): - """ - Solve this linear game and return a :class:`Solution`. - - Returns - ------- - - :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). - - PoorScalingException - If the game could not be solved because CVXOPT crashed while - trying to take the square root of a negative number. - - Examples - -------- - - This example is computed in Gowda and Ravindran in the section - "The value of a Z-transformation":: - - >>> from dunshire import * - >>> K = NonnegativeOrthant(3) - >>> L = [[1,-5,-15],[-1,2,-3],[-12,-15,1]] - >>> e1 = [1,1,1] - >>> e2 = [1,1,1] - >>> SLG = SymmetricLinearGame(L, K, e1, e2) - >>> print(SLG.solution()) - Game value: -6.1724138 - Player 1 optimal: - [ 0.551...] - [-0.000...] - [ 0.448...] - Player 2 optimal: - [0.448...] - [0.000...] - [0.551...] - - The value of the following game can be computed using the fact - that the identity is invertible:: + # The "optimal" and "unknown" results, we actually treat the + # same. Even if CVXOPT bails out due to numerical difficulty, + # it will have some candidate points in mind. If those + # candidates are good enough, we take them. We do the same + # check (perhaps pointlessly so) for "optimal" results. + # + # First we check that the primal/dual objective values are + # close enough (one could be low by ABS_TOL, the other high by + # it) because otherwise CVXOPT might return "unknown" and give + # us two points in the cone that are nowhere near optimal. + if abs(p1_value - p2_value) > 2*options.ABS_TOL: + raise GameUnsolvableException(self, soln_dict) - >>> from dunshire import * - >>> K = NonnegativeOrthant(3) - >>> L = [[1,0,0],[0,1,0],[0,0,1]] - >>> e1 = [1,2,3] - >>> e2 = [4,5,6] - >>> SLG = SymmetricLinearGame(L, K, e1, e2) - >>> print(SLG.solution()) - Game value: 0.0312500 - Player 1 optimal: - [0.031...] - [0.062...] - [0.093...] - Player 2 optimal: - [0.125...] - [0.156...] - [0.187...] + # And we also check that the points it gave us belong to the + # cone, just in case... + if (p1_optimal not in self._K) or (p2_optimal not in self._K): + raise GameUnsolvableException(self, soln_dict) - """ - try: - # First try with a stricter tolerance. Who knows, it might - # work. If it does, we prefer that solution. - return self._try_solution(options.ABS_TOL / 10) - - except (PoorScalingException, GameUnsolvableException): - # Ok, that didn't work. Let's try it with the default tolerance.. - try: - return self._try_solution(options.ABS_TOL / 10) - except (PoorScalingException, GameUnsolvableException) as error: - # Well, that didn't work either. Let's verbosify the matrix - # output format before we allow the exception to be raised. - printing.options['dformat'] = options.DEBUG_FLOAT_FORMAT - raise error + # For the game value, we could use any of: + # + # * p1_value + # * p2_value + # * (p1_value + p2_value)/2 + # * the game payoff + # + # We want the game value to be the payoff, however, so it + # makes the most sense to just use that, even if it means we + # can't test the fact that p1_value/p2_value are close to the + # payoff. + payoff = self.payoff(p1_optimal,p2_optimal) + return Solution(payoff, p1_optimal, p2_optimal) def condition(self): diff --git a/test/failures/001.py b/test/failures/001.py deleted file mode 100644 index 5e45ee8..0000000 --- a/test/failures/001.py +++ /dev/null @@ -1,7 +0,0 @@ -from dunshire import * -L = [[-0.95237953890954685221, 1.83474556206462535712], - [ 1.30481749924621448500, 1.65278664543326403447]] -K = NonnegativeOrthant(2) -e1 = [0.95477167524644313001, 0.63270781756540095397] -e2 = [0.39633793037154141370, 0.10239281495640320530] -SLG = SymmetricLinearGame(L, K, e1, e2) diff --git a/test/failures/002.py b/test/failures/002.py deleted file mode 100644 index e1945b6..0000000 --- a/test/failures/002.py +++ /dev/null @@ -1,7 +0,0 @@ -from dunshire import * -L = [[1.54159395026049472754, 2.21344728574316684799], - [1.33147433507846657541, 1.17913616272988108769]] -K = NonnegativeOrthant(2) -e1 = [0.39903040089404784307, 0.12377403622479113410] -e2 = [0.15695181142215544612, 0.85527381344651265405] -SLG = SymmetricLinearGame(L, K, e1, e2) diff --git a/test/symmetric_linear_game_test.py b/test/symmetric_linear_game_test.py index bba2f7c..f5b48f4 100644 --- a/test/symmetric_linear_game_test.py +++ b/test/symmetric_linear_game_test.py @@ -6,30 +6,13 @@ from unittest import TestCase from dunshire.cones import NonnegativeOrthant from dunshire.games import SymmetricLinearGame -from dunshire.matrices import eigenvalues_re, inner_product +from dunshire.matrices import eigenvalues_re, inner_product, norm from dunshire import options from .randomgen import (RANDOM_MAX, random_icecream_game, random_ll_icecream_game, random_ll_orthant_game, random_nn_scaling, random_orthant_game, random_positive_orthant_game, random_translation) -EPSILON = (1 + RANDOM_MAX)*options.ABS_TOL -""" -This is the tolerance constant including fudge factors that we use to -determine whether or not two numbers are equal in tests. - -Often we will want to compare two solutions, say for games that are -equivalent. If the first game value is low by ``ABS_TOL`` and the second -is high by ``ABS_TOL``, then the total could be off by ``2*ABS_TOL``. We -also subject solutions to translations and scalings, which adds to or -scales their error. If the first game is low by ``ABS_TOL`` and the -second is high by ``ABS_TOL`` before scaling, then after scaling, the -second could be high by ``RANDOM_MAX*ABS_TOL``. That is the rationale -for the factor of ``1 + RANDOM_MAX`` in ``EPSILON``. Since ``1 + -RANDOM_MAX`` is greater than ``2*ABS_TOL``, we don't need to handle the -first issue mentioned (both solutions off by the same amount in opposite -directions). -""" # Tell pylint to shut up about the large number of methods. class SymmetricLinearGameTest(TestCase): # pylint: disable=R0904 @@ -52,22 +35,12 @@ class SymmetricLinearGameTest(TestCase): # pylint: disable=R0904 modifier : float A scaling factor (default: 1) applied to the default - ``EPSILON`` for this comparison. If you have a poorly- + tolerance for this comparison. If you have a poorly- conditioned matrix, for example, you may want to set this greater than one. """ - self.assertTrue(abs(first - second) < EPSILON*modifier) - - - def assert_solution_exists(self, G): - """ - Given a SymmetricLinearGame, ensure that it has a solution. - """ - soln = G.solution() - - expected = G.payoff(soln.player1_optimal(), soln.player2_optimal()) - self.assert_within_tol(soln.game_value(), expected, G.condition()) + self.assertTrue(abs(first - second) < options.ABS_TOL*modifier) @@ -86,40 +59,6 @@ class SymmetricLinearGameTest(TestCase): # pylint: disable=R0904 self.assertTrue(G.condition() >= 1.0) - 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 - some parameters randomly and give it a shot. The resulting - optimal solutions should give us the optimal game value when we - apply the payoff operator to them. - """ - G = random_orthant_game() - self.assert_solution_exists(G) - - - def test_solution_exists_icecream(self): - """ - Like :meth:`test_solution_exists_nonnegative_orthant`, except - over the ice cream cone. - """ - G = random_icecream_game() - self.assert_solution_exists(G) - - - 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, G): """ Test that scaling ``L`` by a nonnegative number scales the value @@ -128,7 +67,8 @@ class SymmetricLinearGameTest(TestCase): # pylint: disable=R0904 (alpha, H) = random_nn_scaling(G) value1 = G.solution().game_value() value2 = H.solution().game_value() - self.assert_within_tol(alpha*value1, value2, H.condition()) + modifier = 4*max(abs(alpha), 1) + self.assert_within_tol(alpha*value1, value2, modifier) def test_scaling_orthant(self): @@ -166,12 +106,11 @@ class SymmetricLinearGameTest(TestCase): # pylint: disable=R0904 (alpha, H) = random_translation(G) value2 = H.solution().game_value() - self.assert_within_tol(value1 + alpha, value2, H.condition()) + modifier = 4*max(abs(alpha), 1) + self.assert_within_tol(value1 + alpha, value2, modifier) # Make sure the same optimal pair works. - self.assert_within_tol(value2, - H.payoff(x_bar, y_bar), - H.condition()) + self.assert_within_tol(value2, H.payoff(x_bar, y_bar), modifier) def test_translation_orthant(self): @@ -210,14 +149,15 @@ class SymmetricLinearGameTest(TestCase): # pylint: disable=R0904 y_bar = soln1.player2_optimal() soln2 = H.solution() - self.assert_within_tol(-soln1.game_value(), - soln2.game_value(), - H.condition()) + # The modifier of 4 is because each could be off by 2*ABS_TOL, + # which is how far apart the primal/dual objectives have been + # observed being. + self.assert_within_tol(-soln1.game_value(), soln2.game_value(), 4) + + # Make sure the switched optimal pair works. Since x_bar and + # y_bar come from G, we use the same modifier. + self.assert_within_tol(soln2.game_value(), H.payoff(y_bar, x_bar), 4) - # Make sure the switched optimal pair works. - self.assert_within_tol(soln2.game_value(), - H.payoff(y_bar, x_bar), - H.condition()) def test_opposite_game_orthant(self): @@ -249,10 +189,10 @@ class SymmetricLinearGameTest(TestCase): # pylint: disable=R0904 value = soln.game_value() ip1 = inner_product(y_bar, G.L()*x_bar - value*G.e1()) - self.assert_within_tol(ip1, 0, G.condition()) + self.assert_within_tol(ip1, 0) ip2 = inner_product(value*G.e2() - G.L().trans()*y_bar, x_bar) - self.assert_within_tol(ip2, 0, G.condition()) + self.assert_within_tol(ip2, 0) def test_orthogonality_orthant(self): @@ -299,20 +239,20 @@ class SymmetricLinearGameTest(TestCase): # pylint: disable=R0904 # fudge factors. eigs = eigenvalues_re(G.L()) - if soln.game_value() > EPSILON: + if soln.game_value() > options.ABS_TOL: # L should be positive stable positive_stable = all([eig > -options.ABS_TOL for eig in eigs]) self.assertTrue(positive_stable) - elif soln.game_value() < -EPSILON: + elif soln.game_value() < -options.ABS_TOL: # L should be negative stable negative_stable = all([eig < options.ABS_TOL for eig in eigs]) self.assertTrue(negative_stable) # The dual game's value should always equal the primal's. + # The modifier of 4 is because even though the games are dual, + # CVXOPT doesn't know that, and each could be off by 2*ABS_TOL. dualsoln = G.dual().solution() - self.assert_within_tol(dualsoln.game_value(), - soln.game_value(), - G.condition()) + self.assert_within_tol(dualsoln.game_value(), soln.game_value(), 4) def test_lyapunov_orthant(self): -- 2.44.2