]> gitweb.michael.orlitzky.com - dunshire.git/commitdiff
Add the condition number of the game to its string representation.
authorMichael Orlitzky <michael@orlitzky.com>
Sun, 30 Oct 2016 21:24:30 +0000 (17:24 -0400)
committerMichael Orlitzky <michael@orlitzky.com>
Sun, 30 Oct 2016 21:24:30 +0000 (17:24 -0400)
dunshire/errors.py
dunshire/games.py

index 1ecf3dc7ffad232bf3fd540e119778f2eb42be95..61fe43d243b65d51595baa09c75e797c002fe306 100644 (file)
@@ -90,7 +90,8 @@ class GameUnsolvableException(Exception):
          e1 = [1.0000000]
               [0.1000000],
          e2 = [3.0000000]
-              [0.1000000].
+              [0.1000000],
+         Condition((L, K, e1, e2)) = 8.311277.
        CVXOPT returned:
          dual infeasibility: None
          dual objective: 1.0
@@ -190,7 +191,8 @@ class PoorScalingException(Exception):
          e1 = [1.0000000]
               [0.1000000],
          e2 = [3.0000000]
-              [0.1000000].
+              [0.1000000],
+         Condition((L, K, e1, e2)) = 8.311277.
        <BLANKLINE>
     """
     def __init__(self, game):
index 9610802e4ffd4108d2c244b2c36a9e44084427cd..46092c380eca141ff993313bd30ec55989a32ed8 100644 (file)
@@ -8,7 +8,7 @@ knows how to solve a linear game.
 from cvxopt import matrix, printing, solvers
 from .cones import CartesianProduct
 from .errors import GameUnsolvableException, PoorScalingException
-from .matrices import append_col, append_row, identity
+from .matrices import append_col, append_row, condition_number, identity
 from . import options
 
 printing.options['dformat'] = options.FLOAT_FORMAT
@@ -221,7 +221,8 @@ class SymmetricLinearGame:
                [ 1],
           e2 = [ 1]
                [ 2]
-               [ 3].
+               [ 3],
+          Condition((L, K, e1, e2)) = 63.669790.
 
     Lists can (and probably should) be used for every argument::
 
@@ -239,7 +240,8 @@ class SymmetricLinearGame:
           e1 = [ 1]
                [ 1],
           e2 = [ 1]
-               [ 1].
+               [ 1],
+          Condition((L, K, e1, e2)) = 3.414214.
 
     The points ``e1`` and ``e2`` can also be passed as some other
     enumerable type (of the correct length) without much harm, since
@@ -261,7 +263,8 @@ class SymmetricLinearGame:
           e1 = [ 1]
                [ 1],
           e2 = [ 1]
-               [ 1].
+               [ 1],
+          Condition((L, K, e1, e2)) = 3.414214.
 
     However, ``L`` will always be intepreted as a list of rows, even
     if it is passed as a :class:`cvxopt.base.matrix` which is
@@ -282,7 +285,8 @@ class SymmetricLinearGame:
           e1 = [ 1]
                [ 1],
           e2 = [ 1]
-               [ 1].
+               [ 1],
+          Condition((L, K, e1, e2)) = 12.147542.
         >>> L = cvxopt.matrix(L)
         >>> print(L)
         [ 1  3]
@@ -297,7 +301,8 @@ class SymmetricLinearGame:
           e1 = [ 1]
                [ 1],
           e2 = [ 1]
-               [ 1].
+               [ 1],
+          Condition((L, K, e1, e2)) = 12.147542.
 
     """
     def __init__(self, L, K, e1, e2):
@@ -319,6 +324,10 @@ class SymmetricLinearGame:
         if not self._e2 in K:
             raise ValueError('the point e2 must lie in the interior of K')
 
+        # Cached result of the self._zero() method.
+        self._zero_col = None
+
+
     def __str__(self):
         """
         Return a string representation of this game.
@@ -327,11 +336,51 @@ class SymmetricLinearGame:
               '  L = {:s},\n' \
               '  K = {!s},\n' \
               '  e1 = {:s},\n' \
-              '  e2 = {:s}.'
+              '  e2 = {:s},\n' \
+              '  Condition((L, K, e1, e2)) = {:f}.'
         indented_L = '\n      '.join(str(self._L).splitlines())
         indented_e1 = '\n       '.join(str(self._e1).splitlines())
         indented_e2 = '\n       '.join(str(self._e2).splitlines())
-        return tpl.format(indented_L, str(self._K), indented_e1, indented_e2)
+
+        return tpl.format(indented_L,
+                          str(self._K),
+                          indented_e1,
+                          indented_e2,
+                          self._condition())
+
+
+    def _zero(self):
+        """
+        Return a column of zeros that fits ``K``.
+
+        This is used in our CVXOPT construction.
+        """
+        if self._zero_col is None:
+            # Cache it, it's constant.
+            self._zero_col = matrix(0, (self._K.dimension(), 1), tc='d')
+        return self._zero_col
+
+
+    def _A(self):
+        """
+        Return the matrix ``A`` used in our CVXOPT construction.
+
+        This matrix ``A``  appears on the right-hand side of ``Ax = b``
+        in the statement of the CVXOPT conelp program.
+        """
+        return matrix([0, self._e2], (1, self._K.dimension() + 1), 'd')
+
+
+    def _G(self):
+        r"""
+        Return the matrix ``G`` used in our CVXOPT construction.
+
+        Thus matrix ``G``that appears on the left-hand side of ``Gx + s = h``
+        in the statement of the CVXOPT conelp program.
+        """
+        I = identity(self._K.dimension())
+        return append_row(append_col(self._zero(), -I),
+                          append_col(self._e1, -self._L))
 
 
     def solution(self):
@@ -407,30 +456,19 @@ class SymmetricLinearGame:
         # Ax = b in the statement of the CVXOPT conelp program.
         b = matrix([1], tc='d')
 
-        # A column of zeros that fits K.
-        zero = matrix(0, (self._K.dimension(), 1), tc='d')
-
         # The column vector "h" that appears on the right-hand side of
         # Gx + s = h in the statement of the CVXOPT conelp program.
-        h = matrix([zero, zero])
+        h = matrix([self._zero(), self._zero()])
 
         # The column vector "c" that appears in the objective function
         # value <c,x> in the statement of the CVXOPT conelp program.
-        c = matrix([-1, zero])
-
-        # The matrix "G" that appears on the left-hand side of Gx + s = h
-        # in the statement of the CVXOPT conelp program.
-        G = append_row(append_col(zero, -identity(self._K.dimension())),
-                       append_col(self._e1, -self._L))
-
-        # The matrix "A" that appears on the right-hand side of Ax = b
-        # in the statement of the CVXOPT conelp program.
-        A = matrix([0, self._e2], (1, self._K.dimension() + 1), 'd')
+        c = matrix([-1, self._zero()])
 
         # Actually solve the thing and obtain a dictionary describing
         # what happened.
         try:
-            soln_dict = solvers.conelp(c, G, h, C.cvxopt_dims(), A, b)
+            soln_dict = solvers.conelp(c, self._G(), h,
+                                       C.cvxopt_dims(), self._A(), b)
         except ValueError as e:
             if str(e) == 'math domain error':
                 # Oops, CVXOPT tried to take the square root of a
@@ -475,6 +513,36 @@ class SymmetricLinearGame:
         return Solution(p1_value, p1_optimal, p2_optimal)
 
 
+    def _condition(self):
+        r"""
+        Return the condition number of this game.
+
+        In the CVXOPT construction of this game, two matrices ``G`` and
+        ``A`` appear. When those matrices are nasty, numerical problems
+        can show up. We define the condition number of this game to be
+        the sum of the condition numbers of ``G`` and ``A`` in the
+        CVXOPT construction. If the condition number of this game is
+        high, then you can expect numerical difficulty (such as
+        :class:`PoorScalingException`).
+
+        Examples
+        --------
+
+        >>> from dunshire import *
+        >>> K = NonnegativeOrthant(1)
+        >>> L = [[1]]
+        >>> e1 = [1]
+        >>> e2 = e1
+        >>> SLG = SymmetricLinearGame(L, K, e1, e2)
+        >>> actual = SLG._condition()
+        >>> expected = 3.6180339887498953
+        >>> abs(actual - expected) < options.ABS_TOL
+        True
+
+        """
+        return condition_number(self._G()) + condition_number(self._A())
+
+
     def dual(self):
         r"""
         Return the dual game to this game.
@@ -504,7 +572,8 @@ class SymmetricLinearGame:
                    [ 3],
               e2 = [ 1]
                    [ 1]
-                   [ 1].
+                   [ 1],
+              Condition((L, K, e1, e2)) = 88.953530.
 
         """
         # We pass ``self._L`` right back into the constructor, because