Add step result with step differences to avoid cancellation errors

pygradflow/newton.py

@@ -8,7 +8,7 @@
 from pygradflow.params import Params, NewtonType
 from pygradflow.problem import Problem
 from pygradflow.step import step_solver
-from pygradflow.step.step_solver import StepSolver
+from pygradflow.step.step_solver import StepSolver, StepResult
 
 
 logger = lgg.getChild("newton")
@@ -29,7 +29,7 @@ def params(self) -> Params:
         return self.orig_iterate.params
 
     @abc.abstractmethod
-    def step(self, iterate):
+    def step(self, iterate: Iterate) -> StepResult:
         raise NotImplementedError()
 
 
@@ -56,14 +56,8 @@ def __init__(
         self.step_solver.update_active_set(active_set)
         self.step_solver.update_derivs(orig_iterate)
 
-    def step(self, iterate: Iterate) -> Iterate:
-        (xn, yn) = self.step_solver.solve(iterate)
-
-        return Iterate(self.problem,
-                       self.params,
-                       xn,
-                       yn,
-                       self.orig_iterate.eval)
+    def step(self, iterate):
+        return self.step_solver.solve(iterate)
 
 
 class FullNewtonMethod(NewtonMethod):
@@ -84,20 +78,14 @@ def __init__(
         super().__init__(problem, orig_iterate, dt, rho)
         self.step_solver = step_solver
 
-    def step(self, iterate: Iterate) -> Iterate:
+    def step(self, iterate: Iterate) -> StepResult:
         p = self.func.projection_initial(iterate, self.rho)
         active_set = self.func.compute_active_set(p)
 
         self.step_solver.update_active_set(active_set)
         self.step_solver.update_derivs(iterate)
 
-        (xn, yn) = self.step_solver.solve(iterate)
-
-        return Iterate(self.problem,
-                       self.params,
-                       xn,
-                       yn,
-                       self.orig_iterate.eval)
+        return self.step_solver.solve(iterate)
 
 
 class FixedActiveSetNewtonMethod(NewtonMethod):
@@ -121,23 +109,22 @@ def __init__(self, problem, active_set, orig_iterate, dt, rho):
             np.sum(active_set),
         )
 
-    def split_sol(self, s):
+    def split_sol(self, s: np.ndarray):
         n = self.problem.num_vars
         m = self.problem.num_cons
 
         assert s.shape == (n + m,)
         return s[:n], s[n:]
 
-    def create_iterate(self, iterate, s):
+    def create_step(self, iterate: Iterate, s: np.ndarray) -> StepResult:
         xn, yn = self.split_sol(s)
         x = iterate.x
         y = iterate.y
 
-        return Iterate(self.problem,
-                       self.params,
-                       x - xn,
-                       y - yn,
-                       self.orig_iterate.eval)
+        dx = x - xn
+        dy = y - yn
+
+        return StepResult(iterate, dx, dy)
 
     @staticmethod
     def active_set_from_iterate(problem, iterate):
@@ -177,15 +164,16 @@ def step(self, iterate):
         solver = sp.sparse.linalg.splu(mat)
 
         s = solver.solve(rhs)
-        next_iterate = self.create_iterate(iterate, s)
+        next_step = self.create_step(iterate, s)
+        next_iterate = next_step.iterate
 
         logger.info(
             "Initial rhs norm: {0}, final: {1}".format(
                 np.linalg.norm(rhs), np.linalg.norm(self.func.value_at(next_iterate))
             )
         )
 
-        return next_iterate
+        return next_step
 
 
 class ActiveSetNewtonMethod(NewtonMethod):
@@ -209,19 +197,13 @@ def __init__(
         self.step_solver = step_solver
         self.step_solver.update_derivs(orig_iterate)
 
-    def step(self, iterate: Iterate) -> Iterate:
+    def step(self, iterate):
         p = self.func.projection_initial(iterate, self.rho)
         active_set = self.func.compute_active_set(p)
 
         self.step_solver.update_active_set(active_set)
 
-        (xn, yn) = self.step_solver.solve(iterate)
-
-        return Iterate(self.problem,
-                       self.params,
-                       xn,
-                       yn,
-                       self.orig_iterate.eval)
+        return self.step_solver.solve(iterate)
 
 
 def newton_method(

pygradflow/solver.py

@@ -82,14 +82,14 @@ def compute_step(
 
         method = newton_method(problem, params, iterate, dt, self.rho)
 
-        def next_iterates():
+        def next_steps():
             curr_iterate = iterate
             while True:
-                next_iterate = method.step(curr_iterate)
-                yield next_iterate
-                curr_iterate = next_iterate
+                next_step = method.step(curr_iterate)
+                yield next_step
+                curr_iterate = next_step.iterate
 
-        return controller.step(iterate, self.rho, dt, next_iterates())
+        return controller.step(iterate, self.rho, dt, next_steps())
 
     def _deriv_check(self, x: np.ndarray, y: np.ndarray) -> None:
         from pygradflow.deriv_check import deriv_check

pygradflow/step/scaled_step_solver.py

@@ -4,7 +4,7 @@
 import numpy as np
 
 from pygradflow.implicit_func import ScaledImplicitFunc
-from pygradflow.step.step_solver import StepSolver
+from pygradflow.step.step_solver import StepSolver, StepResult
 from pygradflow.iterate import Iterate
 from pygradflow.params import Params
 from pygradflow.problem import Problem
@@ -68,7 +68,7 @@ def update_active_set(self, active_set: np.ndarray) -> None:
         self.active_set = copy.copy(active_set)
         self.reset_deriv()
 
-    def solve(self, iterate: Iterate) -> Tuple[np.ndarray, np.ndarray]:
+    def solve(self, iterate: Iterate) -> StepResult:
         (b0, b1, b2) = self.initial_rhs(iterate)
 
         n = self.n
@@ -90,7 +90,4 @@ def solve(self, iterate: Iterate) -> Tuple[np.ndarray, np.ndarray]:
         dx = sx
         dy = fact * (sy - rho * b2)
 
-        x = iterate.x
-        y = iterate.y
-
-        return (x - dx, y - dy)
+        return StepResult(iterate, dx, dy)

pygradflow/step/standard_step_solver.py

@@ -4,7 +4,7 @@
 import numpy as np
 
 from pygradflow.implicit_func import ImplicitFunc
-from pygradflow.step.step_solver import StepSolver
+from pygradflow.step.step_solver import StepSolver, StepResult
 from pygradflow.iterate import Iterate
 from pygradflow.params import Params
 from pygradflow.problem import Problem
@@ -54,7 +54,7 @@ def update_active_set(self, active_set: np.ndarray) -> None:
         self.active_set = copy.copy(active_set)
         self._reset_deriv()
 
-    def solve(self, iterate: Iterate) -> Tuple[np.ndarray, np.ndarray]:
+    def solve(self, iterate: Iterate) -> StepResult:
         if self.deriv is None:
             self._compute_deriv()
 
@@ -73,7 +73,4 @@ def solve(self, iterate: Iterate) -> Tuple[np.ndarray, np.ndarray]:
         dx = s[:n]
         dy = s[n:]
 
-        x = iterate.x
-        y = iterate.y
-
-        return (x - dx, y - dy)
+        return StepResult(iterate, dx, dy)

pygradflow/step/step_control.py

@@ -5,14 +5,15 @@
 
 from pygradflow.log import logger
 
+from pygradflow.step.step_solver import StepResult
 from pygradflow.controller import ControllerSettings, LogController
 from pygradflow.implicit_func import ImplicitFunc
 from pygradflow.iterate import Iterate
 from pygradflow.params import Params
 from pygradflow.problem import Problem
 
 
-class StepResult:
+class StepControlResult:
     def __init__(self, iterate: Iterate, lamb: float, accepted: bool) -> None:
         self.iterate = iterate
         self.lamb = lamb
@@ -26,15 +27,19 @@ def __init__(self, problem: Problem, params: Params) -> None:
         self.lamb = params.lamb_init
 
     @abc.abstractmethod
-    def step(self, iterate, rho, dt, next_iterates):
+    def step(self,
+             iterate: Iterate,
+             rho: float,
+             dt: float,
+             next_steps: Iterator[StepResult]) -> StepControlResult:
         raise NotImplementedError()
 
 
 class ExactController(StepController):
     def __init__(self, problem, params):
         super().__init__(problem, params)
 
-    def step(self, iterate, rho, dt, next_iterates):
+    def step(self, iterate, rho, dt, next_steps):
         assert dt > 0.0
         lamb = 1.0 / dt
 
@@ -46,20 +51,20 @@ def func_val(iterate):
         cur_func_val = func_val(iterate)
 
         for i in range(10):
-            next_iterate = next(next_iterates)
+            next_iterate = next(next_steps).iterate
 
             next_func_val = func_val(next_iterate)
             logger.info(f"Func val: {next_func_val}")
 
             if next_func_val <= self.params.newton_tol:
                 logger.debug("Newton method converged in %d iterations", i + 1)
-                return StepResult(next_iterate, 0.5 * lamb, True)
+                return StepControlResult(next_iterate, 0.5 * lamb, True)
             elif next_func_val > cur_func_val:
                 break
 
         logger.debug("Newton method did not converge")
 
-        return StepResult(next_iterate, 2.0 * lamb, False)
+        return StepControlResult(next_iterate, 2.0 * lamb, False)
 
 
 class DistanceRatioController(StepController):
@@ -68,9 +73,7 @@ def __init__(self, problem: Problem, params: Params) -> None:
         settings = ControllerSettings.from_params(params)
         self.controller = LogController(settings, params.theta_ref)
 
-    def step(
-        self, iterate: Iterate, rho: float, dt: float, next_iterates: Iterator[Iterate]
-    ) -> StepResult:
+    def step(self, iterate, rho, dt, next_steps):
         assert dt > 0.0
         lamb = 1.0 / dt
 
@@ -79,24 +82,30 @@ def step(
 
         func = ImplicitFunc(problem, iterate, dt)
 
-        mid_iterate = next(next_iterates)
+        # TODO: Fix
+        mid_step = next(next_steps)
+        mid_iterate = mid_step.iterate
 
-        if np.linalg.norm(func.value_at(mid_iterate, rho)) <= params.newton_tol:
+        mid_func_norm = np.linalg.norm(func.value_at(mid_iterate, rho))
+
+        if mid_func_norm <= params.newton_tol:
             lamb_n = max(lamb * params.lamb_red, params.lamb_min)
-            logger.info("Newton converged during first iteration, lamb_n = %f", lamb_n)
-            return StepResult(mid_iterate, lamb_n, True)
+            logger.info("Newton converged during first iteration, lamb_n = %f",
+                        lamb_n)
+            return StepControlResult(mid_iterate, lamb_n, True)
 
-        first_diff = mid_iterate.dist(iterate)
+        first_diff = mid_step.diff
 
         if first_diff == 0.:
-            return StepResult(mid_iterate, lamb, True)
+            return StepControlResult(mid_iterate, lamb, True)
 
-        final_iterate = next(next_iterates)
+        final_step = next(next_steps)
+        final_iterate = final_step.iterate
 
-        second_diff = final_iterate.dist(mid_iterate)
+        second_diff = final_step.diff
 
         if second_diff == 0.:
-            return StepResult(final_iterate, lamb, True)
+            return StepControlResult(final_iterate, lamb, True)
 
         theta = second_diff / first_diff
 
@@ -120,4 +129,4 @@ def step(
 
         self.lamb = lamb_n
 
-        return StepResult(final_iterate, lamb_n, accepted)
+        return StepControlResult(final_iterate, lamb_n, accepted)

pygradflow/step/step_solver.py

@@ -1,4 +1,5 @@
 import abc
+import functools
 
 import numpy as np
 import scipy as sp
@@ -7,6 +8,33 @@
 from pygradflow.params import Params
 from pygradflow.problem import Problem
 from pygradflow.step.linear_solver import LinearSolver
+from pygradflow.util import norm_mult
+
+
+class StepResult:
+    def __init__(self, orig_iterate, dx, dy):
+        self.orig_iterate = orig_iterate
+        self.dx = dx
+        self.dy = dy
+
+    @functools.cached_property
+    def iterate(self):
+        iterate = self.orig_iterate
+
+        return Iterate(iterate.problem,
+                       iterate.params,
+                       iterate.x - self.dx,
+                       iterate.y - self.dy,
+                       iterate.eval)
+
+    @functools.cached_property
+    def diff(self):
+        diff = norm_mult(self.dx, self.dy)
+
+        assert np.allclose(self.iterate.dist(self.orig_iterate),
+                           diff)
+
+        return diff
 
 
 class StepSolver(abc.ABC):
@@ -24,10 +52,10 @@ def linear_solver(self, mat: sp.sparse.spmatrix) -> LinearSolver:
 
     @abc.abstractmethod
     def update_active_set(self, active_set: np.ndarray):
-        raise NotImplementedError
+        raise NotImplementedError()
 
     def update_derivs(self, iterate: Iterate):
-        raise NotImplementedError
+        raise NotImplementedError()
 
-    def solve(self, iterate: Iterate):
-        raise NotImplementedError
+    def solve(self, iterate: Iterate) -> StepResult:
+        raise NotImplementedError()

pyproject.toml

@@ -1,6 +1,6 @@
 [tool.poetry]
 name = "pygradflow"
-version = "0.2.1"
+version = "0.2.2"
 description = "PyGradFlow is a simple implementation of the sequential homotopy method to be used to solve general nonlinear programs."
 authors = ["Christoph Hansknecht <[email protected]>"]
 readme = "README.md"