From: Michael Orlitzky Date: Tue, 11 Oct 2016 15:24:30 +0000 (-0400) Subject: Attempt to recover from "unknown" solutions. X-Git-Tag: 0.1.0~160 X-Git-Url: https://gitweb.michael.orlitzky.com/?a=commitdiff_plain;h=6aff913423ca1f1ec2434650a6065d18a2f6415f;p=dunshire.git Attempt to recover from "unknown" solutions. Sometimes CVXOPT was coming back with "unknown" solutions that looked pretty close to real solutions. Now, instead of failing immediately, we check to see if the primal/dual objectives are within ABS_TOL of each other, and whether or not the optimal solutions lie in the cone. If they all do, we take that as our solution and succeed rather than raising an exception. This required the loosening of what it means to be in the cone. Before this commit, a point needed to be "far" inside the cone (past the tolerance). Now it can lie outside of the cone (but still within the tolerance). This weakens the checks that we perform on e1 and e2, but gives us an easy way to check the optimal strategies. --- diff --git a/src/dunshire/cones.py b/src/dunshire/cones.py index d7c9b62..870a882 100644 --- a/src/dunshire/cones.py +++ b/src/dunshire/cones.py @@ -134,12 +134,6 @@ class NonnegativeOrthant(SymmetricCone): distinguish between strict and non-strict containment -- the test uses a tolerance parameter. - This test will err on the side of caution, and return ``True`` - only when the ``point`` lies inside this cone *past* the - tolerance guarantee. That means if you ask whether or not a - boundary point lies in this cone, you will get ``False`` as your - answer. - Parameters ---------- @@ -178,10 +172,10 @@ class NonnegativeOrthant(SymmetricCone): >>> matrix([1,-0.1,3]) in K False - However, not even a boundary point is considered inside of ``K``: + A boundary point is considered inside of ``K``: >>> K = NonnegativeOrthant(3) >>> matrix([1,0,3]) in K - False + True Junk arguments don't work: @@ -203,7 +197,7 @@ class NonnegativeOrthant(SymmetricCone): if not point.size == (self.dimension(), 1): raise TypeError('the given point has the wrong dimensions') - return all([x > options.ABS_TOL for x in point]) + return all([x > -options.ABS_TOL for x in point]) @@ -236,12 +230,6 @@ class IceCream(SymmetricCone): distinguish between strict and non-strict containment -- the test uses a tolerance parameter. - This test will err on the side of caution, and return ``True`` - only when the ``point`` lies inside this cone *past* the - tolerance guarantee. That means if you ask whether or not a - boundary point lies in this cone, you will get ``False`` as your - answer. - Parameters ---------- @@ -274,11 +262,11 @@ class IceCream(SymmetricCone): >>> matrix([1,0.5,0.5]) in K True - But this one lies on its boundary and runs foul of the tolerance: + This one lies on its boundary: >>> K = IceCream(3) >>> matrix([1,0,1]) in K - False + True This point lies entirely outside of the ice cream cone: @@ -310,10 +298,10 @@ class IceCream(SymmetricCone): if self.dimension() == 1: # In one dimension, the ice cream cone is the nonnegative # orthant. - return height > options.ABS_TOL + return height > -options.ABS_TOL else: radius = point[1:] - return norm(radius) < (height - options.ABS_TOL) + return norm(radius) < (height + options.ABS_TOL) @@ -358,12 +346,6 @@ class SymmetricPSD(SymmetricCone): distinguish between strict and non-strict containment -- the test uses a tolerance parameter. - This test will err on the side of caution, and return ``True`` - only when the ``point`` lies inside this cone *past* the - tolerance guarantee. That means if you ask whether or not a - boundary point lies in this cone, you will get ``False`` as your - answer. - Parameters ---------- @@ -409,21 +391,28 @@ class SymmetricPSD(SymmetricCone): >>> A in K True - But any matrix with a zero eigenvalue will lie on the boundary - of the Symmetric PSD cone and run foul of our tolerance: + Boundary points lie in the cone as well: >>> K = SymmetricPSD(2) >>> matrix([[0,0],[0,0]]) in K - False + True >>> K = SymmetricPSD(5) >>> A = matrix([[1,0,0,0,0], - ... [0,1,0,0,0], - ... [0,0,0,0,0], - ... [0,0,0,1,0], - ... [0,0,0,0,1]]) + ... [0,1,0,0,0], + ... [0,0,0,0,0], + ... [0,0,0,1,0], + ... [0,0,0,0,1]]) >>> A in K - False + True + + However, this matrix has a negative eigenvalue: + + >>> K = SymmetricPSD(2) + >>> A = matrix([[1,-2], + ... [-2,1]]) + >>> A in K + False Junk arguments don't work: @@ -446,7 +435,7 @@ class SymmetricPSD(SymmetricCone): raise TypeError('the given point has the wrong dimensions') if not point.typecode == 'd': point = matrix(point, (self.dimension(), self.dimension()), 'd') - return all([e > options.ABS_TOL for e in eigenvalues(point)]) + return all([e > -options.ABS_TOL for e in eigenvalues(point)]) @@ -525,7 +514,7 @@ class CartesianProduct(SymmetricCone): >>> K = CartesianProduct(NonnegativeOrthant(3), IceCream(3)) >>> (matrix([0,0,0]), matrix([1,0,1])) in K - False + True >>> K = CartesianProduct(NonnegativeOrthant(3), IceCream(3)) >>> (matrix([1,1,1]), matrix([1,1,1])) in K diff --git a/src/dunshire/games.py b/src/dunshire/games.py index 4b5383d..e85fe63 100644 --- a/src/dunshire/games.py +++ b/src/dunshire/games.py @@ -433,21 +433,35 @@ 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. @@ -494,6 +508,48 @@ class SymmetricLinearGameTest(TestCase): Tests for the SymmetricLinearGame and Solution classes. """ + def random_orthant_params(self): + """ + Generate the ``L``, ``K``, ``e1``, and ``e2`` parameters for a + random game over the nonnegative orthant. + """ + 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())] + return (L, K, e1, e2) + + + def random_icecream_params(self): + """ + Generate the ``L``, ``K``, ``e1``, and ``e2`` parameters for a + random game 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())] + + return (L, K, e1, e2) + + def assert_within_tol(self, first, second): """ Test that ``first`` and ``second`` are equal within our default @@ -522,12 +578,7 @@ 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) = self.random_orthant_params() self.assert_solution_exists(L, K, e1, e2) def test_solution_exists_ice_cream(self): @@ -535,25 +586,7 @@ class SymmetricLinearGameTest(TestCase): 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) = self.random_icecream_params() self.assert_solution_exists(L, K, e1, e2) @@ -576,12 +609,8 @@ class SymmetricLinearGameTest(TestCase): of the game by the same number. Use the nonnegative orthant as our cone. """ - 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 = matrix([[uniform(-10, 10) for i in range(K.dimension())] - for j in range(K.dimension())]) + (L, K, e1, e2) = self.random_orthant_params() + L = matrix(L) # So that we can scale it by alpha below. G1 = SymmetricLinearGame(L, K, e1, e2) value1 = G1.solution().game_value() alpha = uniform(0.1, 10) @@ -596,25 +625,8 @@ class SymmetricLinearGameTest(TestCase): The same test as :meth:`test_nonnegative_scaling_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 = matrix([[uniform(-10, 10) for i in range(K.dimension())] - for j in range(K.dimension())]) + (L, K, e1, e2) = self.random_icecream_params() + L = matrix(L) # So that we can scale it by alpha below. G1 = SymmetricLinearGame(L, K, e1, e2) value1 = G1.solution().game_value()