Merge pull request #30 from fides-dev/develop

Fides 0.4.0

Author: FFroehlich

LICENSE.md

@@ -1,11 +1,11 @@
 # License conditions
 
-## pytr
+## fides
 
-pytr is released under the 3-Clause BSD License (BSD-3-Clause) with the
+fides is released under the 3-Clause BSD License (BSD-3-Clause) with the
 following terms:
 
-Copyright (c) 2020, Fabian Fröhlich
+Copyright (c) 2020-2021, Fabian Fröhlich
 All rights reserved.
 
 Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:

README.md

@@ -11,18 +11,20 @@
 
 ## About Fides
 
-Fides implements an Interior Trust Region Reflective for boundary costrained
-optimization problems based on the papers [ColemanLi1994] and [ColemanLi1996
-]. Accordingly, Fides is named after the Roman goddess of trust and
+Fides implements an Interior Trust Region Reflective for boundary constrained
+optimization problems based on the papers
+[ColemanLi1994](https://doi.org/10.1007/BF01582221) and
+[ColemanLi1996](http://dx.doi.org/10.1137/0806023). Accordingly, Fides is named
+after the Roman goddess of trust and
 reliability. In contrast to other optimizers, Fides solves the full trust
--region subproblem exactly, which can yields higher quality proposal steps, but
+-region subproblem exactly, which can yield higher quality proposal steps, but
 is computationally more expensive. This makes Fides particularly attractive
 for optimization problems with objective functions that are computationally
 expensive to evaluate and the computational cost of solving the trust
  -region subproblem is negligible.
 
 Fides can be installed via `pip install fides`. Further documentation is
- avaliable at [Read the Docs](fides-optimizer.readthedocs.io).
+ available at [Read the Docs](https://fides-optimizer.readthedocs.io/).
 
 
 ## Features

fides/constants.py

@@ -39,6 +39,7 @@ class SubSpaceDim(str, enum.Enum):
     """
     TWO = '2D'  #: Two dimensional Newton/Gradient subspace
     FULL = 'full'  #: Full :math:`\mathbb{R}^n`
+    STEIHAUG = 'scg'  #: CG subspace via Steihaug's method
 
 
 class StepBackStrategy(str, enum.Enum):

fides/hessian_approximation.py

@@ -3,7 +3,7 @@
 -------------------------
 This module provides various generic Hessian approximation strategies that
 can be employed when the calculating the exact Hessian or an approximation
-is computationally too demandind.
+is computationally too demanding.
 """
 
 
@@ -21,7 +21,8 @@ def __init__(self, hess_init: Optional[np.ndarray] = None):
         Create a Hessian update strategy instance
 
         :param hess_init:
-            Inital guess for the Hessian, if empty Identity matrix will be used
+            Initial guess for the Hessian, if empty Identity matrix will be
+            used
         """
         self.hess_init = None
         if hess_init is not None:
@@ -33,7 +34,8 @@ def set_init(self, hess_init: np.ndarray):
         Create a Hessian update strategy instance
 
         :param hess_init:
-            Inital guess for the Hessian, if empty Identity matrix will be used
+            Initial guess for the Hessian, if empty Identity matrix will be
+            used
         """
         if not isinstance(hess_init, np.ndarray):
             raise ValueError('Cannot initialize with hess_init of type'
@@ -62,7 +64,7 @@ def init_mat(self, dim: int):
         else:
             self._hess = self.hess_init.copy()
             if self._hess.shape[0] != dim:
-                raise ValueError('Inital approximation had inconsistent '
+                raise ValueError('Initial approximation had inconsistent '
                                  f'dimension, was {self._hess.shape[0]}, '
                                  f'but should be {dim}.')
 

fides/minimize.py

@@ -128,15 +128,16 @@ def _reset(self):
 
     def minimize(self, x0: np.ndarray):
         """
-        Minimize the objective function the interior trust-region reflective
-        algorithm described by [ColemanLi1994] and [ColemanLi1996]
+        Minimize the objective function using the interior trust-region
+        reflective algorithm described by [ColemanLi1994] and [ColemanLi1996]
         Convergence with respect to function value is achieved when
         math:`|f_{k+1} - f_k|` < options[`fatol`] - :math:`f_k` options[
-        `frtol`]. Similarly,  convergence with respect to optimization
+        `frtol`]. Similarly, convergence with respect to optimization
         variables is achieved when :math:`||x_{k+1} - x_k||` < options[
-        `xatol`] - :math:`x_k` options[`xrtol`].  Convergence with respect
-        to the gradient is achieved when :math:`||g_k||` <
-        options[`gatol`] or `||g_k||` < options[`grtol`] * `f_k`.  Other than
+        `xtol`] :math:`x_k` (note that this is checked in transformed
+        coordinates that account for distance to boundaries).  Convergence
+        with respect to the gradient is achieved when :math:`||g_k||` <
+        options[`gatol`] or `||g_k||` < options[`grtol`] * `f_k`. Other than
         that, optimization can be terminated when iterations exceed
         options[ `maxiter`] or the elapsed time is expected to exceed
         options[`maxtime`] on the next iteration.
@@ -194,12 +195,6 @@ def minimize(self, x0: np.ndarray):
                              f'x has {len(self.x)} entries but gradient has '
                              f'{len(self.grad)}!')
 
-        if not len(self.grad) == len(self.x):
-            raise ValueError('Provided objective function must return a '
-                             'gradient vector of the same shape as x, '
-                             f'x has {len(self.x)} entries but gradient has '
-                             f'{len(self.grad)}!')
-
         # hessian approximation would error on these earlier
         if not self.hess.ndim == 2:
             raise ValueError('Provided objective function must return a '
@@ -212,7 +207,7 @@ def minimize(self, x0: np.ndarray):
 
         if not self.hess.shape[0] == len(self.x):
             raise ValueError('Provided objective function must return a '
-                             'square Hessian matrix with same dimension'
+                             'square Hessian matrix with same dimension as x. '
                              f'x has {len(self.x)} entries but Hessian has '
                              f'{self.hess.shape[0]}!')
 
@@ -365,16 +360,17 @@ def update_tr_radius(self,
 
             # values as proposed in algorithm 4.1 in Nocedal & Wright
             if self.tr_ratio >= self.get_option(Options.ETA) \
-                    and not interior_solution:
+                    and not interior_solution and step.qpval <= 0:
                 # increase radius
                 self.delta = self.get_option(Options.GAMMA2) * self.delta
-            elif self.tr_ratio <= self.get_option(Options.MU):
+            elif self.tr_ratio <= self.get_option(Options.MU) or \
+                    step.qpval > 0:
                 # decrease radius
                 self.delta = np.nanmin([
                     self.delta * self.get_option(Options.GAMMA1),
                     nsx / 4
                 ])
-            return self.tr_ratio > 0.0
+            return self.tr_ratio > 0.0 and step.qpval <= 0
 
     def check_convergence(self, step: Step, fval: float,
                           grad: np.ndarray) -> None:
@@ -475,8 +471,7 @@ def check_continue(self) -> bool:
             self.exitflag = ExitFlag.DELTA_TOO_SMALL
             self.logger.warning(
                 f'Stopping as trust region radius {self.delta:.2E} is '
-                f'smaller '
-                'than machine precision.'
+                f'smaller than machine precision.'
             )
             return False
 
@@ -498,7 +493,7 @@ def make_non_degenerate(self, eps=1e2 * np.spacing(1)) -> None:
 
     def get_affine_scaling(self) -> Tuple[np.ndarray, np.ndarray]:
         """
-        Computes the vector v and dv, the diagonal of it's Jacobian. For the
+        Computes the vector v and dv, the diagonal of its Jacobian. For the
         definition of v, see Definition 2 in [Coleman-Li1994]
 
         :return:

fides/stepback.py

@@ -27,7 +27,7 @@ def stepback_reflect(tr_step: Step,
     Compute new proposal steps according to a reflection strategy.
 
     :param tr_step:
-        Reference trust region step that will be reflect
+        Reference trust region step that will be reflected
     :param x:
         Current values of the optimization variables
     :param sg:

fides/steps.py

@@ -8,33 +8,21 @@
 
 import numpy as np
 import scipy.linalg as linalg
+import scipy.sparse.linalg as slinalg
 
 from numpy.linalg import norm
 from scipy.sparse import csc_matrix
 from scipy.optimize import Bounds, NonlinearConstraint, minimize
 
 from logging import Logger
 from .subproblem import (
-    solve_1d_trust_region_subproblem, solve_nd_trust_region_subproblem
+    solve_1d_trust_region_subproblem, solve_nd_trust_region_subproblem,
+    get_1d_trust_region_boundary_solution, quadratic_form
 )
 
 from typing import Union
 
 
-def quadratic_form(Q: np.ndarray, p: np.ndarray, x: np.ndarray) -> float:
-    """
-    Computes the quadratic form :math:`x^TQx + x^Tp`
-
-    :param Q: Matrix
-    :param p: Vector
-    :param x: Input
-
-    :return:
-        Value of form
-    """
-    return 0.5 * x.T.dot(Q).dot(x) + p.T.dot(x)
-
-
 def normalize(v: np.ndarray) -> None:
     """
     Inplace normalization of a vector
@@ -57,7 +45,7 @@ class Step:
         transformed step ss: `ss = subspace * sc`
     :ivar ss: Affine transformed step: `s = scaling * ss`
     :ivar og_s: `s` without step back
-    :ivar og_sc: `st` without step back
+    :ivar og_sc: `sc` without step back
     :ivar og_ss: `ss` without step back
     :ivar sg: Rescaled gradient `scaling * g`
     :ivar hess: Hessian of the objective function at `x`
@@ -238,7 +226,7 @@ class TRStepFull(Step):
     the trust region subproblem.
     """
 
-    type = 'trnd'
+    type = 'nd'
 
     def __init__(self, x, sg, hess, scaling, g_dscaling, delta, theta,
                  ub, lb, logger):
@@ -253,7 +241,7 @@ class TRStep2D(Step):
     the trust region subproblem according to a 2D subproblem
     """
 
-    type = 'tr2d'
+    type = '2d'
 
     def __init__(self, x, sg, hess, scaling, g_dscaling, delta, theta,
                  ub, lb, logger):
@@ -270,8 +258,11 @@ def __init__(self, x, sg, hess, scaling, g_dscaling, delta, theta,
                 # in this case we are in Case 2 of Fig 12 in
                 # [Coleman-Li1994]
                 logger.debug('Newton direction did not have negative '
-                             'curvature adding scaling * np.sign(sg) to '
-                             '2D subspace.')
+                             'using scaling * np.sign(sg) and ev to smallest '
+                             'eigenvalue instead.')
+                e, v = slinalg.eigs(self.shess, k=1, which='SR')
+                s_newt = v[:, np.argmin(e)]
+                normalize(s_newt)
                 s_grad = scaling * np.sign(sg) + (sg == 0)
             else:
                 s_grad = sg.copy()
@@ -290,6 +281,78 @@ def __init__(self, x, sg, hess, scaling, g_dscaling, delta, theta,
         self.subspace = np.expand_dims(s_newt, 1)
 
 
+class CGStep(Step):
+    """
+        This class provides the machinery to compute an approximate solution of
+        the trust region subproblem using the Steihaug Method
+        """
+
+    type = 'cg'
+
+    def calculate(self):
+        nsg = norm(self.sg)
+        self.conj_grad(min(0.5, np.sqrt(nsg)) * nsg)
+        self.s = self.scaling.dot(self.ss + self.ss0)
+        self.step_back()
+        self.qpval = quadratic_form(self.shess, self.sg, self.ss + self.ss0)
+
+    def conj_grad(self, eps):
+        """
+        Compute step proposal using conjugate gradient method
+
+        :param eps:
+            tolerance for residual norm
+        """
+        raise NotImplementedError()
+
+
+class TRStepSteihaug(CGStep):
+    """
+    This class provides the machinery to compute an approximate solution of
+    the trust region subproblem using the Steihaug Method
+    """
+
+    type = 'cgs'
+
+    def conj_grad(self, eps):
+        z = np.zeros_like(self.sg)
+        r = self.sg.copy()
+        d = -self.sg.copy()
+        if norm(r) < eps:
+            self.ss = z
+            return
+
+        while True:
+            bd = self.shess.dot(d)
+            c = d.dot(bd)
+            r2 = np.dot(r, r)
+            alpha = r2 / c
+            zp = z + alpha * d
+            if c <= 0 or norm(zp) >= self.delta:
+                self.subspace = np.expand_dims(d, 1)
+                self.ss0 = z
+                self.sc = get_1d_trust_region_boundary_solution(
+                    self.shess, self.sg, self.subspace[:, 0], self.ss0,
+                    self.delta
+                ) * np.ones((1,))
+                self.ss = self.subspace.dot(self.sc)
+                return
+            rp = r + alpha*bd
+            rp2 = np.dot(rp, rp)
+            if np.sqrt(rp2) < eps:
+                normalize(d)
+                self.subspace = np.expand_dims(d, 1)
+                self.sc = zp.dot(d) * np.ones((1,))
+                self.ss = self.subspace.dot(self.sc)
+                self.ss0 = zp - self.ss
+                return
+            beta = rp2 / r2
+
+            d = -rp + beta * d
+            z = zp
+            r = rp
+
+
 class TRStepReflected(Step):
     """
     This class provides the machinery to compute a reflected step based on
@@ -374,6 +437,22 @@ def __init__(self, x, sg, hess, scaling, g_dscaling, delta, theta,
         self.subspace = np.expand_dims(s_grad, 1)
 
 
+class ScaledGradientStep(Step):
+    """
+    This class provides the machinery to compute a gradient step.
+    """
+
+    type = 'dg'
+
+    def __init__(self, x, sg, hess, scaling, g_dscaling, delta, theta,
+                 ub, lb, logger):
+        super().__init__(x, sg, hess, scaling, g_dscaling, delta, theta,
+                         ub, lb, logger)
+        s_grad = scaling*sg.copy()
+        normalize(s_grad)
+        self.subspace = np.expand_dims(s_grad, 1)
+
+
 class RefinedStep(Step):
     """
     This class provides the machinery to refine a step based on interior
@@ -393,8 +472,8 @@ def __init__(self, x, sg, hess, scaling, g_dscaling, delta, theta,
             NonlinearConstraint(
                 fun=lambda xs: (norm(xs) - delta) * np.ones((1,)),
                 jac=lambda xs: np.expand_dims(xs, 1).T / norm(xs),
-                lb=np.zeros((1,)),
-                ub=np.ones((1,)) * np.inf,
+                lb=-np.ones((1,)) * np.inf,
+                ub=np.zeros((1,)),
             )
         ]
         self.guess = step.ss + step.ss0