Merge pull request #100 from chrhansk/feature-active-set

Add different active set computations

docs/conf.py

@@ -14,7 +14,7 @@
 project = "pygradflow"
 copyright = "2023, Christoph Hansknecht"
 author = "Christoph Hansknecht"
-release = "0.5.7"
+release = "0.5.8"
 
 # -- General configuration ---------------------------------------------------
 # https://www.sphinx-doc.org/en/master/usage/configuration.html#general-configuration

pygradflow/implicit_func.py

@@ -66,11 +66,15 @@ def value_at(
         raise NotImplementedError()
 
     @abc.abstractmethod
-    def projection_initial(self, iterate: Iterate, rho: float) -> np.ndarray:
+    def projection_initial(self, iterate: Iterate, rho: float, tau=None) -> np.ndarray:
         raise NotImplementedError()
 
+    def compute_active_set(self, iterate: Iterate, rho: float, tau=None) -> np.ndarray:
+        p = self.projection_initial(iterate, rho, tau)
+        return self.active_set_at_point(p)
+
     @abc.abstractmethod
-    def compute_active_set(self, x: np.ndarray) -> np.ndarray:
+    def active_set_at_point(self, p: np.ndarray) -> np.ndarray:
         raise NotImplementedError()
 
     def apply_project_deriv(self, mat: sp.sparse.spmatrix, active_set: np.ndarray):
@@ -106,15 +110,25 @@ def project(self, x: np.ndarray, active_set: np.ndarray):
 
         return super().project_box(x, lb, ub, active_set)
 
-    def compute_active_set(self, x: np.ndarray):
+    def active_set_at_point(self, p):
         problem = self.problem
         lb = problem.var_lb
         ub = problem.var_ub
-        return super().compute_active_set_box(x, lb, ub)
+        return super().compute_active_set_box(p, lb, ub)
 
-    def projection_initial(self, iterate: Iterate, rho: float):
+    def projection_initial(self, iterate: Iterate, rho: float, tau=None):
         x_0 = self.orig_iterate.x
         dt = self.dt
+
+        if tau is not None:
+            lamb = 1.0 / dt
+            x = iterate.x
+            return (
+                (1.0 - tau * lamb) * x
+                + (tau * lamb) * x_0
+                - tau * iterate.aug_lag_deriv_x(rho)
+            )
+
         return x_0 - dt * iterate.aug_lag_deriv_x(rho)
 
     # @override
@@ -125,7 +139,7 @@ def value_at(self, iterate, rho, active_set=None):
         p = self.projection_initial(iterate, rho)
 
         if active_set is None:
-            active_set = self.compute_active_set(p)
+            active_set = self.compute_active_set(iterate, rho)
 
         xval = iterate.x - self.project(p, active_set)
         yval = iterate.y - (y_0 + dt * iterate.aug_lag_deriv_y())
@@ -162,8 +176,7 @@ def deriv_at(
         self, iterate: Iterate, rho: float, active_set: Optional[np.ndarray] = None
     ):
         if active_set is None:
-            p = self.projection_initial(iterate, rho)
-            active_set = self.compute_active_set(p)
+            active_set = self.compute_active_set(iterate, rho)
 
         hess = iterate.aug_lag_deriv_xx(rho)
         jac = iterate.aug_lag_deriv_xy()
@@ -187,24 +200,33 @@ def value_at(self, iterate, rho, active_set=None):
         p = self.projection_initial(iterate, rho)
 
         if active_set is None:
-            active_set = self.compute_active_set(p)
+            active_set = self.compute_active_set(iterate, rho)
 
         xval = lamb * iterate.x - self.project(p, active_set)
         yval = -(lamb * iterate.y - (lamb * y_0 + iterate.aug_lag_deriv_y()))
 
         return np.concatenate([xval, yval])
 
-    def projection_initial(self, iterate: Iterate, rho: float):
+    def projection_initial(self, iterate: Iterate, rho: float, tau=None):
         x_0 = self.orig_iterate.x
         lamb = self.lamb
 
+        if tau is not None:
+            lamb = 1.0 / self.dt
+            x = iterate.x
+            f_x = lamb * (1 - tau * lamb)
+            f_x0 = tau * lamb * lamb
+            f_d = tau * lamb
+
+            return f_x * x + f_x0 * x_0 - f_d * iterate.aug_lag_deriv_x(rho)
+
         return lamb * x_0 - iterate.aug_lag_deriv_x(rho)
 
     def project(self, x: np.ndarray, active_set: np.ndarray):
         return super().project_box(x, self.lb, self.ub, active_set)
 
-    def compute_active_set(self, x: np.ndarray):
-        return super().compute_active_set_box(x, self.lb, self.ub)
+    def active_set_at_point(self, p):
+        return super().compute_active_set_box(p, self.lb, self.ub)
 
     def deriv(
         self, jac: sp.sparse.spmatrix, hess: sp.sparse.spmatrix, active_set: np.ndarray
@@ -241,8 +263,7 @@ def deriv_at(
         self, iterate: Iterate, rho: float, active_set: Optional[np.ndarray] = None
     ):
         if active_set is None:
-            p = self.projection_initial(iterate, rho)
-            active_set = self.compute_active_set(p)
+            active_set = self.compute_active_set(iterate, rho)
 
         hess = iterate.aug_lag_deriv_xx(rho)
         jac = iterate.aug_lag_deriv_xy()

pygradflow/iterate.py

@@ -196,3 +196,13 @@ def cons_nonlin(self, other: "Iterate") -> np.ndarray:
             problem = self.problem
             return np.zeros((problem.num_cons,))
         return (other.cons - next_cons) / dx_dot
+
+    def check_eval(self):
+        # Check for and raise EvalError if any of the
+        # evaluations are fail
+        self.obj
+        self.obj_grad
+
+        if self.problem.num_cons > 0:
+            self.cons
+            self.cons_jac

pygradflow/linear_solver/mumps_solver.py

@@ -15,6 +15,9 @@ def __init__(self, mat: sp.sparse.spmatrix, symmetric=False) -> None:
 
         self.ctx = mumps.DMumpsContext(sym=sym, par=1)
 
+        # Only output error messages
+        self.ctx.set_icntl(4, 1)
+
         self.ctx.set_icntl(13, 1)
 
         self.mat = mat

pygradflow/newton.py

@@ -15,12 +15,13 @@
 
 class NewtonMethod(abc.ABC):
     def __init__(
-        self, problem: Problem, orig_iterate: Iterate, dt: float, rho: float
+        self, problem: Problem, orig_iterate: Iterate, dt: float, rho: float, tau=None
     ) -> None:
         self.problem = problem
         self.orig_iterate = orig_iterate
         self.dt = dt
         self.rho = rho
+        self.tau = tau
 
     @property
     def params(self) -> Params:
@@ -44,13 +45,13 @@ def __init__(
         dt: float,
         rho: float,
         step_solver: StepSolver,
+        tau=None,
     ) -> None:
-        super().__init__(problem, orig_iterate, dt, rho)
+        super().__init__(problem, orig_iterate, dt, rho, tau)
         self.func = step_solver.func
 
         self.step_solver = step_solver
-        p = self.func.projection_initial(orig_iterate, rho)
-        active_set = self.func.compute_active_set(p)
+        active_set = self.func.compute_active_set(orig_iterate, rho, tau)
 
         self.step_solver.update_active_set(active_set)
         self.step_solver.update_derivs(orig_iterate)
@@ -73,14 +74,14 @@ def __init__(
         dt: float,
         rho: float,
         step_solver: StepSolver,
+        tau=None,
     ) -> None:
-        super().__init__(problem, orig_iterate, dt, rho)
+        super().__init__(problem, orig_iterate, dt, rho, tau)
         self.func = step_solver.func
         self.step_solver = step_solver
 
     def step(self, iterate: Iterate) -> StepResult:
-        p = self.func.projection_initial(iterate, self.rho)
-        active_set = self.func.compute_active_set(p)
+        active_set = self.func.compute_active_set(iterate, self.rho, self.tau)
 
         self.step_solver.update_active_set(active_set)
         self.step_solver.update_derivs(iterate)
@@ -95,8 +96,8 @@ class FixedActiveSetNewtonMethod(NewtonMethod):
     at each step.
     """
 
-    def __init__(self, problem, active_set, orig_iterate, dt, rho):
-        super().__init__(problem, orig_iterate, dt, rho)
+    def __init__(self, problem, active_set, orig_iterate, dt, rho, tau=None):
+        super().__init__(problem, orig_iterate, dt, rho, tau)
         self.func = step_solver.func
 
         assert active_set.dtype == bool
@@ -192,16 +193,16 @@ def __init__(
         dt: float,
         rho: float,
         step_solver: StepSolver,
+        tau=None,
     ) -> None:
-        super().__init__(problem, orig_iterate, dt, rho)
+        super().__init__(problem, orig_iterate, dt, rho, tau)
         self.func = step_solver.func
 
         self.step_solver = step_solver
         self.step_solver.update_derivs(orig_iterate)
 
     def step(self, iterate):
-        p = self.func.projection_initial(iterate, self.rho)
-        active_set = self.func.compute_active_set(p)
+        active_set = self.func.compute_active_set(iterate, self.rho, self.tau)
 
         self.step_solver.update_active_set(active_set)
 
@@ -221,15 +222,15 @@ def __init__(
         dt: float,
         rho: float,
         step_solver: StepSolver,
+        tau=None,
     ) -> None:
-        super().__init__(problem, orig_iterate, dt, rho)
+        super().__init__(problem, orig_iterate, dt, rho, tau)
         self.func = step_solver.func
         self.step_solver = step_solver
 
     def _set_iterate(self, iterate):
         self.step_solver.update_derivs(iterate)
-        p = self.func.projection_initial(iterate, self.rho)
-        active_set = self.func.compute_active_set(p)
+        active_set = self.func.compute_active_set(iterate, self.rho, self.tau)
         self.step_solver.update_active_set(active_set)
 
     def step(self, iterate):
@@ -285,26 +286,29 @@ def step(self, iterate):
 
         logger.debug("Line search converged in %d iterations", it + 1)
 
-        next_x = iterate.x + dx
-        active_set = self.func.compute_active_set(next_x)
+        step_result = StepResult(self.orig_iterate, dx, dy, active_set=None, rcond=None)
 
-        return StepResult(self.orig_iterate, dx, dy, active_set, rcond=None)
+        step_result.active_set = self.func.compute_active_set(
+            step_result.iterate, self.rho, self.tau
+        )
+
+        return step_result
 
 
 def newton_method(
-    problem: Problem, params: Params, iterate: Iterate, dt: float, rho: float
+    problem: Problem, params: Params, iterate: Iterate, dt: float, rho: float, tau=None
 ) -> NewtonMethod:
     assert dt > 0.0
     assert rho > 0.0
 
     solver = step_solver(problem, params, iterate, dt, rho)
 
     if params.newton_type == NewtonType.Simplified:
-        return SimplifiedNewtonMethod(problem, iterate, dt, rho, solver)
+        return SimplifiedNewtonMethod(problem, iterate, dt, rho, solver, tau)
     elif params.newton_type == NewtonType.Full:
-        return FullNewtonMethod(problem, iterate, dt, rho, solver)
+        return FullNewtonMethod(problem, iterate, dt, rho, solver, tau)
     elif params.newton_type == NewtonType.ActiveSet:
-        return ActiveSetNewtonMethod(problem, iterate, dt, rho, solver)
+        return ActiveSetNewtonMethod(problem, iterate, dt, rho, solver, tau)
     else:
         assert params.newton_type == NewtonType.Globalized
-        return GlobalizedNewtonMethod(problem, iterate, dt, rho, solver)
+        return GlobalizedNewtonMethod(problem, iterate, dt, rho, solver, tau)

pygradflow/params.py

@@ -11,6 +11,12 @@
 import numpy as np
 
 
+class ActiveSetType(Enum):
+    Standard = auto()
+    SmallestActiveSet = auto()
+    LargestActiveSet = auto()
+
+
 class NewtonType(Enum):
     """
     The Newton method to be used to solve the semi-smooth systems.
@@ -213,6 +219,9 @@ class Params:
 
     local_infeas_tol: float = 1e-8
 
+    active_set_type: ActiveSetType = ActiveSetType.Standard
+    active_set_method: Optional[Callable[..., float]] = None
+
     newton_type: NewtonType = NewtonType.Simplified
     newton_tol: float = 1e-8
 

pygradflow/solver.py

@@ -4,6 +4,7 @@
 
 from pygradflow.callbacks import Callbacks, CallbackType
 from pygradflow.display import Format, StateData, print_problem_stats, solver_display
+from pygradflow.eval import EvalError
 from pygradflow.iterate import Iterate
 from pygradflow.log import logger
 from pygradflow.params import Params
@@ -229,6 +230,11 @@ def solve(
 
         iterate = self.transform.initial_iterate
 
+        try:
+            iterate.check_eval()
+        except EvalError as e:
+            raise Exception("Failed to evaluate initial iterate") from e
+
         print_problem_stats(problem, iterate)
 
         lamb = params.lamb_init