Merge pull request #17 from chrhansk/rcond-reporting

Report rcond in solver

pygradflow/display.py

@@ -1,5 +1,6 @@
 from termcolor import colored
 
+from pygradflow.params import Params
 from pygradflow.problem import Problem
 
 
@@ -95,8 +96,7 @@ def __call__(self, state):
         return getattr(iterate, self.name)
 
 
-def problem_display(problem: Problem):
-
+def problem_display(problem: Problem, params: Params):
     is_bounded = problem.var_bounded
 
     cols = []
@@ -112,6 +112,10 @@ def problem_display(problem: Problem):
     cols.append(Column("Primal step", 16, "{:16.8e}", StateAttr("primal_step_norm")))
     cols.append(Column("Dual step", 16, "{:16.8e}", StateAttr("dual_step_norm")))
     cols.append(Column("Lambda", 16, "{:16.8e}", StateAttr("lamb")))
+
+    if params.report_rcond:
+        cols.append(Column("Rcond", 5, "{:5.0e}", StateAttr("rcond")))
+
     cols.append(Column("Type", 8, StepFormatter(), StateAttr("step_accept")))
 
     return Display(cols)

pygradflow/params.py

@@ -90,8 +90,12 @@ class Params:
     time_limit: float = np.inf
     display_interval: float = 0.1
 
+    # lower bound on objective function value,
+    # used to detect unbounded problems
     obj_lower_limit: float = -1e10
 
+    report_rcond: bool = False
+
     @property
     def dtype(self):
         return np.float32 if self.precision == Precision.Single else np.float64

pygradflow/solver.py

@@ -44,6 +44,9 @@ def __init__(self, x, y, d, status):
         self.d = d
         self.status = status
 
+    def __repr__(self):
+        return "SolverResult(status={0})".format(self.status)
+
     @property
     def success(self):
         return SolverStatus.success(self.status)
@@ -123,7 +126,7 @@ def print_result(self,
 
         desc = "{:>30s}".format(status.description)
 
-        status_desc = Format.redgreen(desc, status.success, bold=True)
+        status_desc = Format.redgreen(desc, SolverStatus.success(status), bold=True)
         status_name = Format.bold("{:>30s}".format("Status"))
 
         logger.info("%30s: %30s", status_name, status_desc)
@@ -145,7 +148,7 @@ def solve(self, x_0: np.ndarray, y_0: np.ndarray) -> SolverResult:
         params = self.params
         dtype = params.dtype
 
-        display = problem_display(problem)
+        display = problem_display(problem, params)
 
         x = x_0.astype(dtype)
         y = y_0.astype(dtype)
@@ -234,6 +237,7 @@ def solve(self, x_0: np.ndarray, y_0: np.ndarray) -> SolverResult:
                 state["dual_step_norm"] = lambda: dual_step_norm
                 state["lamb"] = lambda: lamb
                 state["step_accept"] = lambda: accept
+                state["rcond"] = lambda: step_result.rcond
 
                 logger.info(display.row(state))
 

pygradflow/step/cond_estimate.py

@@ -0,0 +1,97 @@
+import math
+import numpy as np
+import scipy as sp
+
+from pygradflow.params import Params
+from pygradflow.log import logger
+from pygradflow.step.linear_solver import LinearSolver
+
+seed = 42
+
+
+class ConditionEstimator:
+    """
+    Estimates condition of matrix using generic lienar solver.
+    Based on
+    "Estimating Extremal Eigenvalues and Condition Numbers of Matrices"
+    by Dixon
+    """
+
+    def __init__(self,
+                 mat: sp.sparse.spmatrix,
+                 solver: LinearSolver,
+                 params: Params,
+                 min_prob: float = .99,
+                 factor: float = 10.):
+        assert 0 < min_prob < 1
+        self.size = mat.shape[0]
+
+        self.params = params
+        assert mat.shape == (self.size, self.size)
+
+        self.min_prob = min_prob
+        self.factor = factor
+        self.mat = mat
+        self.linear_solver = sp.sparse.linalg.splu(mat)
+        self.rng = np.random.default_rng(seed=seed)
+
+    def _required_its(self):
+        factor = (1. - self.min_prob) / 1.6 * math.pow(self.size, -0.5)
+        return -2 * math.ceil(math.log(factor, self.factor))
+
+    def _random_vec(self):
+        rng = self.rng
+        size = self.size
+        vec = rng.normal(size=size)
+
+        while True:
+            if (vec != 0.).any():
+                break
+            vec = rng.normal(size=size)
+
+        norm = np.linalg.norm(vec)
+        return (vec / norm).astype(self.params.dtype)
+
+    def estimate_rcond(self) -> float:
+        mat = self.mat
+        trans_mat = mat.T
+        linear_solver = self.linear_solver
+
+        num_its = self._required_its()
+        assert num_its > 0
+
+        x = self._random_vec()
+        y = self._random_vec()
+
+        xprod = np.copy(x)
+        yprod = np.copy(y)
+
+        logger.debug("Number of iterations: %s", num_its)
+
+        for k in range(num_its):
+
+            xprod = mat @ xprod
+            xprod = trans_mat @ xprod
+
+            yprod = linear_solver.solve(yprod, trans='T')
+            yprod = linear_solver.solve(yprod)
+
+        pow_fac = 1. / (2. * num_its)
+
+        xdot = x.dot(xprod)
+        xdot = math.pow(xdot, pow_fac)
+
+        ydot = y.dot(yprod)
+        ydot = math.pow(ydot, pow_fac)
+
+        if np.isinf(xdot) or np.isinf(ydot):
+            return 0.
+
+        cond = xdot * ydot
+
+        if np.isinf(cond):
+            return 0.
+
+        rcond = 1./cond
+
+        return rcond

pygradflow/step/extended_step_solver.py

@@ -1,5 +1,3 @@
-from typing import Tuple
-
 import numpy as np
 import scipy as sp
 
@@ -85,9 +83,8 @@ def _compute_deriv(self) -> None:
         assert self.deriv.shape == (n + m, n + m)
         assert self.deriv.dtype == self.params.dtype
 
-    def solve_scaled(
-        self, b0: np.ndarray, b1: np.ndarray, b2t: np.ndarray
-    ) -> Tuple[np.ndarray, np.ndarray]:
+    def solve_scaled(self, b0, b1, b2t):
+        params = self.params
         if self.deriv is None:
             self._compute_deriv()
 
@@ -106,4 +103,8 @@ def solve_scaled(
         dx = s[:n]
         dy = s[n:]
 
-        return (dx, dy)
+        rcond = None
+        if params.report_rcond:
+            rcond = self.estimate_rcond(self.deriv, self.solver)
+
+        return (dx, dy, rcond)

pygradflow/step/linear_solver.py

@@ -2,35 +2,46 @@
 
 from pygradflow.params import LinearSolverType
 from numpy import ndarray
-from scipy.sparse._csc import csc_matrix
 
 
 class LinearSolver:
-    def solve(self, b: ndarray) -> ndarray:
+    def solve(self, b: ndarray, trans: bool = False) -> ndarray:
         raise NotImplementedError
 
 
 class MINRESSolver(LinearSolver):
     def __init__(self, mat):
         self.mat = mat
 
-    def solve(self, b):
+    def solve(self, b, trans=False):
+        # matrix should be symmetric anyways
         return sp.sparse.linalg.minres(self.mat, b)[0]
 
 
 class GMRESSolver(LinearSolver):
-    def __init__(self, mat: csc_matrix) -> None:
+    def __init__(self, mat: sp.sparse.spmatrix) -> None:
         self.mat = mat
 
-    def solve(self, b: ndarray) -> ndarray:
-        return sp.sparse.linalg.gmres(self.mat, b, atol=0.0)[0]
+    def solve(self, b: ndarray, trans=False) -> ndarray:
+        mat = self.mat.T if trans else self.mat
+        return sp.sparse.linalg.gmres(mat, b, atol=0.0)[0]
+
+
+class LUSolver(LinearSolver):
+    def __init__(self, mat: sp.sparse.spmatrix) -> None:
+        self.mat = mat
+        self.solver = sp.sparse.linalg.splu(mat)
+
+    def solve(self, b: ndarray, trans=False) -> ndarray:
+        trans_str = 'T' if trans else 'N'
+        return self.solver.solve(b, trans=trans_str)
 
 
 def linear_solver(
     mat: sp.sparse.spmatrix, solver_type: LinearSolverType
 ) -> LinearSolver:
     if solver_type == LinearSolverType.LU:
-        return sp.sparse.linalg.splu(mat)
+        return LUSolver(mat)
     elif solver_type == LinearSolverType.MINRES:
         return MINRESSolver(mat)
     else:

pygradflow/step/scaled_step_solver.py

@@ -1,5 +1,5 @@
 import copy
-from typing import Tuple
+from typing import Tuple, Optional
 
 import numpy as np
 
@@ -52,8 +52,8 @@ def initial_rhs(
 
         return (b0, b1, b2)
 
-    def solve_scaled(self, b0, b1, b2t):
-        raise NotImplementedError
+    def solve_scaled(self, b0, b1, b2t) -> Tuple[np.ndarray, np.ndarray, Optional[float]]:
+        raise NotImplementedError()
 
     def reset_deriv(self) -> None:
         self.deriv = None
@@ -82,12 +82,12 @@ def solve(self, iterate: Iterate) -> StepResult:
 
         b2t = fact * b2
 
-        (sx, sy) = self.solve_scaled(b0, b1, b2t)
+        (sx, sy, rcond) = self.solve_scaled(b0, b1, b2t)
 
         assert sx.shape == (n,)
         assert sy.shape == (m,)
 
         dx = sx
         dy = fact * (sy - rho * b2)
 
-        return StepResult(iterate, dx, dy)
+        return StepResult(iterate, dx, dy, rcond)

pygradflow/step/standard_step_solver.py

@@ -1,5 +1,4 @@
 import copy
-from typing import Tuple
 
 import numpy as np
 
@@ -73,4 +72,10 @@ def solve(self, iterate: Iterate) -> StepResult:
         dx = s[:n]
         dy = s[n:]
 
-        return StepResult(iterate, dx, dy)
+        params = self.params
+
+        rcond = None
+        if params.report_rcond:
+            rcond = self.estimate_rcond(self.deriv, self.solver)
+
+        return StepResult(iterate, dx, dy, rcond)