Fixed default cheap initialisation behaviour for inverse problems

Author: lindonroberts

dfols/solver.py

@@ -94,7 +94,7 @@ def __str__(self):
 
 
 def solve_main(objfun, x0, args, xl, xu, npt, rhobeg, rhoend, maxfun, nruns_so_far, nf_so_far, nx_so_far, nsamples, params,
-               diagnostic_info, scaling_changes, r0_avg_old=None, r0_nsamples_old=None):
+               diagnostic_info, scaling_changes, r0_avg_old=None, r0_nsamples_old=None, default_growing_method_set_by_user=None):
     # Evaluate at x0 (keep nf, nx correct and check for f < 1e-12)
     # The hard bit is determining what m = len(r0) should be, and allocating memory appropriately
     if r0_avg_old is None:
@@ -142,6 +142,17 @@ def solve_main(objfun, x0, args, xl, xu, npt, rhobeg, rhoend, maxfun, nruns_so_f
         num_samples_run = r0_nsamples_old
         nf = nf_so_far
         nx = nx_so_far
+    
+    # On the first run, set default growing method (unless the user has already done this)
+    if default_growing_method_set_by_user is not None and (not default_growing_method_set_by_user):
+        # If m>=n, the default growing method (use_full_rank_interp) is best
+        # However, this can fail for m<n, so need to use an alternative method (perturb_trust_region_step)
+        if m < len(x0):
+            logging.debug("Inverse problem (m<n), switching default growing method")
+            params('growing.full_rank.use_full_rank_interp', new_value=False)
+            params('growing.perturb_trust_region_step', new_value=True)
+            if not params.params_changed['growing.delta_scale_new_dirns']:
+                params('growing.delta_scale_new_dirns', new_value=0.1)
 
     # Initialise controller
     control = Controller(objfun, args, x0, r0_avg, num_samples_run, xl, xu, npt, rhobeg, rhoend, nf, nx, maxfun, params, scaling_changes)
@@ -836,6 +847,11 @@ def solve(objfun, x0, args=(), bounds=None, npt=None, rhobeg=None, rhoend=1e-8,
     if user_params is not None:
         for (key, val) in user_params.items():
             params(key, new_value=val)
+    
+    # Default growing method depends on if m>=n or m<n - need to set this once we know m
+    # But should only do this if the user hasn't forced a choice on us
+    default_growing_method_set_by_user = user_params is not None and \
+        ('growing.full_rank.use_full_rank_interp' in user_params or 'growing.perturb_trust_region_step' in user_params)
 
     scaling_changes = None
     if scaling_within_bounds:
@@ -948,7 +964,7 @@ def solve(objfun, x0, args=(), bounds=None, npt=None, rhobeg=None, rhoend=1e-8,
     nx = 0
     xmin, rmin, fmin, jacmin, nsamples_min, nf, nx, nruns, exit_info, diagnostic_info = \
         solve_main(objfun, x0, args, xl, xu, npt, rhobeg, rhoend, maxfun, nruns, nf, nx, nsamples, params,
-                    diagnostic_info, scaling_changes)
+                    diagnostic_info, scaling_changes, default_growing_method_set_by_user=default_growing_method_set_by_user)
 
     # Hard restarts loop
     last_successful_run = nruns

dfols/tests/test_solver.py

@@ -133,3 +133,41 @@ def runTest(self):
         self.assertTrue(array_compare(soln.jacobian, rosenbrock_jacobian(soln.x), thresh=3e-0), "Wrong Jacobian")
         self.assertTrue(abs(soln.f) < 1e-10, "Wrong fmin")
 
+
+class TestInverseProblem(unittest.TestCase):
+    # Minimise an inverse problem
+    def runTest(self):
+        # Simple problem with global minimum at the origin, for instance
+        objfun = lambda x: np.array([np.sin(x[0])**2, np.sin(x[1])**2, np.sum(np.sin(x[2:]**2))])
+        jac = lambda x: np.array([[2*np.sin(x[0])*np.cos(x[0]), 0, 0, 0, 0], 
+                                  [0, 2*np.sin(x[1])*np.cos(x[1]), 0, 0, 0], 
+                                  [0, 0, 2*np.sin(x[2])*np.cos(x[2]), 2*np.sin(x[3])*np.cos(x[3]), 2*np.sin(x[4])*np.cos(x[4])]])  # for n=5 only
+        x0 = np.ones((5,))
+        np.random.seed(0)
+        soln = dfols.solve(objfun, x0)
+        self.assertTrue(array_compare(soln.x, np.zeros((5,)), thresh=1e-3), "Wrong xmin")
+        self.assertTrue(array_compare(soln.resid, objfun(soln.x), thresh=1e-10), "Wrong resid")
+        print(soln.jacobian)
+        print(jac(soln.x))
+        self.assertTrue(array_compare(soln.jacobian, jac(soln.x), thresh=1e-2), "Wrong Jacobian")
+        self.assertTrue(abs(soln.f) < 1e-10, "Wrong fmin")
+
+
+class TestInverseProblemGrowing(unittest.TestCase):
+    # Minimise an inverse problem, with growing
+    def runTest(self):
+        # Simple problem with global minimum at the origin, for instance
+        objfun = lambda x: np.array([np.sin(x[0])**2, np.sin(x[1])**2, np.sum(np.sin(x[2:]**2))])
+        jac = lambda x: np.array([[2*np.sin(x[0])*np.cos(x[0]), 0, 0, 0, 0], 
+                                  [0, 2*np.sin(x[1])*np.cos(x[1]), 0, 0, 0], 
+                                  [0, 0, 2*np.sin(x[2])*np.cos(x[2]), 2*np.sin(x[3])*np.cos(x[3]), 2*np.sin(x[4])*np.cos(x[4])]])  # for n=5 only
+        x0 = np.ones((5,))
+        np.random.seed(0)
+        soln = dfols.solve(objfun, x0, user_params={'growing.ndirs_initial':2})
+        self.assertTrue(array_compare(soln.x, np.zeros((5,)), thresh=1e-3), "Wrong xmin")
+        self.assertTrue(array_compare(soln.resid, objfun(soln.x), thresh=1e-10), "Wrong resid")
+        print(soln.jacobian)
+        print(jac(soln.x))
+        self.assertTrue(array_compare(soln.jacobian, jac(soln.x), thresh=1e-2), "Wrong Jacobian")
+        self.assertTrue(abs(soln.f) < 1e-10, "Wrong fmin")
+

dfols/version.py

@@ -22,4 +22,4 @@
 
 """
 
-__version__ = '1.0.2'
+__version__ = '1.1'

docs/advanced.rst

@@ -85,22 +85,22 @@ Multiple Restarts
 
 Dynamically Growing Initial Set
 -------------------------------
-* :code:`growing.ndirs_initial` - Number of initial points to add (excluding :math:`x_k`). Default is :code:`npt-1`. 
-* :code:`growing.num_new_dirns_each_iter` - Number of new search directions to add with each iteration where we do not have a full set of search directions. Default is 1. 
-* :code:`growing.delta_scale_new_dirns` - When adding new search directions, the length of the step as a multiple of :math:`\Delta_k`. Default is 1, but if setting :code:`growing.perturb_trust_region_step=True` should be made smaller (e.g. 0.1). 
-* :code:`growing.do_geom_steps` - While still growing the initial set, whether to do geometry-improving steps in the trust region algorithm, as per the usual algorithm. Default is :code:`False`. 
-* :code:`growing.safety.do_safety_step` - While still growing the initial set, whether to perform safety steps, or the regular trust region steps. Default is :code:`True`. 
-* :code:`growing.safety.reduce_delta` - While still growing the initial set, whether to reduce :math:`\Delta_k` in safety steps. Default is :code:`False`. 
-* :code:`growing.safety.full_geom_step` - While still growing the initial set, whether to do a full geometry-improving step within safety steps (the same as the post-growing phase of the algorithm). Since this involves reducing :math:`\Delta_k`, cannot be :code:`True` if :code:`growing.safety.reduce_delta` is :code:`True`. Default is :code:`False`. 
-* :code:`growing.full_rank.use_full_rank_interp` - Whether to perturb the interpolated :math:`J_k` to make it full rank, allowing the trust region step to include components in the full search space. If :code:`True`, setting :code:`growing.num_new_dirns_each_iter` to 0 is recommended. Default is :code:`False`. 
+* :code:`growing.ndirs_initial` - Number of initial points to add (excluding :math:`x_k`). This should only be changed to a value less than :math:`n`, and only if the default setup cost of :math:`n+1` evaluations of :code:`objfun` is impractical. If this is set to be less than the default, the input value :code:`npt` should be set to :math:`n`. If the default is used, all the below parameters have no effect on DFO-LS. Default is :code:`npt-1`. 
+* :code:`growing.full_rank.use_full_rank_interp` - If :code:`growing.ndirs_initial` is less than :code:`npt`, whether to perturb the interpolated :math:`J_k` to make it full rank, allowing the trust region step to include components in the full search space. Default is :code:`True` if :math:`m\geq n` and :code:`False` otherwise (opposite to :code:`growing.perturb_trust_region_step`). 
+* :code:`growing.perturb_trust_region_step` - Whether to perturb the trust region step by an orthogonal direction not yet searched. This is an alternative to :code:`growing.full_rank.use_full_rank_interp`. Default is :code:`False` if :math:`m\geq n` and :code:`True` otherwise (opposite to :code:`growing.full_rank.use_full_rank_interp`).
+* :code:`growing.delta_scale_new_dirns` - When adding new search directions, the length of the step as a multiple of :math:`\Delta_k`. Default is 1, or 0.1 if :code:`growing.perturb_trust_region_step=True`. 
 * :code:`growing.full_rank.scale_factor` - Magnitude of extra components added to :math:`J_k`. Default is :math:`10^{-2}`. 
 * :code:`growing.full_rank.svd_scale_factor` - Floor singular values of :math:`J_k` at this factor of the last nonzero value. Default is 1.
 * :code:`growing.full_rank.min_sing_val` - Absolute floor on singular values of :math:`J_k`. Default is :math:`10^{-6}`.
 * :code:`growing.full_rank.svd_max_jac_cond` - Cap on condition number of :math:`J_k` after applying floors to singular values (effectively another floor on the smallest singular value, since the largest singular value is fixed). Default is :math:`10^8`.
+* :code:`growing.do_geom_steps` - While still growing the initial set, whether to do geometry-improving steps in the trust region algorithm, as per the usual algorithm. Default is :code:`False`. 
+* :code:`growing.safety.do_safety_step` - While still growing the initial set, whether to perform safety steps, or the regular trust region steps. Default is :code:`True`. 
+* :code:`growing.safety.reduce_delta` - While still growing the initial set, whether to reduce :math:`\Delta_k` in safety steps. Default is :code:`False`. 
+* :code:`growing.safety.full_geom_step` - While still growing the initial set, whether to do a full geometry-improving step within safety steps (the same as the post-growing phase of the algorithm). Since this involves reducing :math:`\Delta_k`, cannot be :code:`True` if :code:`growing.safety.reduce_delta` is :code:`True`. Default is :code:`False`. 
 * :code:`growing.reset_delta` - Whether or not to reset trust region radius :math:`\Delta_k` to its initial value at the end of the growing phase. Default is :code:`False`.
 * :code:`growing.reset_rho` - Whether or not to reset trust region radius lower bound :math:`\rho_k` to its initial value at the end of the growing phase. Default is :code:`False`.
 * :code:`growing.gamma_dec` - Trust region decrease parameter during the growing phase. Default is :code:`tr_radius.gamma_dec`.
-* :code:`growing.perturb_trust_region_step` - Whether to perturb the trust region step by an orthogonal direction not yet searched. This is an alternative to :code:`growing.full_rank.use_full_rank_interp`. Default is :code:`False`.
+* :code:`growing.num_new_dirns_each_iter` - Number of new search directions to add with each iteration where we do not have a full set of search directions. Default is 0, as this approach is not recommended. 
 
 
 References

docs/build/html/.buildinfo

@@ -1,4 +1,4 @@
 # Sphinx build info version 1
 # This file hashes the configuration used when building these files. When it is not found, a full rebuild will be done.
-config: aa64240bdba1bf7c73696f2a8006ee55
+config: f6da778e78f2e0fc2cb0ba72618f9760
 tags: 645f666f9bcd5a90fca523b33c5a78b7

docs/build/html/_sources/advanced.rst.txt

@@ -85,22 +85,22 @@ Multiple Restarts
 
 Dynamically Growing Initial Set
 -------------------------------
-* :code:`growing.ndirs_initial` - Number of initial points to add (excluding :math:`x_k`). Default is :code:`npt-1`. 
-* :code:`growing.num_new_dirns_each_iter` - Number of new search directions to add with each iteration where we do not have a full set of search directions. Default is 1. 
-* :code:`growing.delta_scale_new_dirns` - When adding new search directions, the length of the step as a multiple of :math:`\Delta_k`. Default is 1, but if setting :code:`growing.perturb_trust_region_step=True` should be made smaller (e.g. 0.1). 
-* :code:`growing.do_geom_steps` - While still growing the initial set, whether to do geometry-improving steps in the trust region algorithm, as per the usual algorithm. Default is :code:`False`. 
-* :code:`growing.safety.do_safety_step` - While still growing the initial set, whether to perform safety steps, or the regular trust region steps. Default is :code:`True`. 
-* :code:`growing.safety.reduce_delta` - While still growing the initial set, whether to reduce :math:`\Delta_k` in safety steps. Default is :code:`False`. 
-* :code:`growing.safety.full_geom_step` - While still growing the initial set, whether to do a full geometry-improving step within safety steps (the same as the post-growing phase of the algorithm). Since this involves reducing :math:`\Delta_k`, cannot be :code:`True` if :code:`growing.safety.reduce_delta` is :code:`True`. Default is :code:`False`. 
-* :code:`growing.full_rank.use_full_rank_interp` - Whether to perturb the interpolated :math:`J_k` to make it full rank, allowing the trust region step to include components in the full search space. If :code:`True`, setting :code:`growing.num_new_dirns_each_iter` to 0 is recommended. Default is :code:`False`. 
+* :code:`growing.ndirs_initial` - Number of initial points to add (excluding :math:`x_k`). This should only be changed to a value less than :math:`n`, and only if the default setup cost of :math:`n+1` evaluations of :code:`objfun` is impractical. If this is set to be less than the default, the input value :code:`npt` should be set to :math:`n`. If the default is used, all the below parameters have no effect on DFO-LS. Default is :code:`npt-1`. 
+* :code:`growing.full_rank.use_full_rank_interp` - If :code:`growing.ndirs_initial` is less than :code:`npt`, whether to perturb the interpolated :math:`J_k` to make it full rank, allowing the trust region step to include components in the full search space. Default is :code:`True` if :math:`m\geq n` and :code:`False` otherwise (opposite to :code:`growing.perturb_trust_region_step`). 
+* :code:`growing.perturb_trust_region_step` - Whether to perturb the trust region step by an orthogonal direction not yet searched. This is an alternative to :code:`growing.full_rank.use_full_rank_interp`. Default is :code:`False` if :math:`m\geq n` and :code:`True` otherwise (opposite to :code:`growing.full_rank.use_full_rank_interp`).
+* :code:`growing.delta_scale_new_dirns` - When adding new search directions, the length of the step as a multiple of :math:`\Delta_k`. Default is 1, or 0.1 if :code:`growing.perturb_trust_region_step=True`. 
 * :code:`growing.full_rank.scale_factor` - Magnitude of extra components added to :math:`J_k`. Default is :math:`10^{-2}`. 
 * :code:`growing.full_rank.svd_scale_factor` - Floor singular values of :math:`J_k` at this factor of the last nonzero value. Default is 1.
 * :code:`growing.full_rank.min_sing_val` - Absolute floor on singular values of :math:`J_k`. Default is :math:`10^{-6}`.
 * :code:`growing.full_rank.svd_max_jac_cond` - Cap on condition number of :math:`J_k` after applying floors to singular values (effectively another floor on the smallest singular value, since the largest singular value is fixed). Default is :math:`10^8`.
+* :code:`growing.do_geom_steps` - While still growing the initial set, whether to do geometry-improving steps in the trust region algorithm, as per the usual algorithm. Default is :code:`False`. 
+* :code:`growing.safety.do_safety_step` - While still growing the initial set, whether to perform safety steps, or the regular trust region steps. Default is :code:`True`. 
+* :code:`growing.safety.reduce_delta` - While still growing the initial set, whether to reduce :math:`\Delta_k` in safety steps. Default is :code:`False`. 
+* :code:`growing.safety.full_geom_step` - While still growing the initial set, whether to do a full geometry-improving step within safety steps (the same as the post-growing phase of the algorithm). Since this involves reducing :math:`\Delta_k`, cannot be :code:`True` if :code:`growing.safety.reduce_delta` is :code:`True`. Default is :code:`False`. 
 * :code:`growing.reset_delta` - Whether or not to reset trust region radius :math:`\Delta_k` to its initial value at the end of the growing phase. Default is :code:`False`.
 * :code:`growing.reset_rho` - Whether or not to reset trust region radius lower bound :math:`\rho_k` to its initial value at the end of the growing phase. Default is :code:`False`.
 * :code:`growing.gamma_dec` - Trust region decrease parameter during the growing phase. Default is :code:`tr_radius.gamma_dec`.
-* :code:`growing.perturb_trust_region_step` - Whether to perturb the trust region step by an orthogonal direction not yet searched. This is an alternative to :code:`growing.full_rank.use_full_rank_interp`. Default is :code:`False`.
+* :code:`growing.num_new_dirns_each_iter` - Number of new search directions to add with each iteration where we do not have a full set of search directions. Default is 0, as this approach is not recommended. 
 
 
 References

docs/build/html/_sources/history.rst.txt

@@ -13,5 +13,10 @@ Version 1.0.1 (20 Feb 2018)
 Version 1.0.2 (20 Jun 2018)
 ---------------------------
 * Extra optional input :code:`args` which passes through arguments for :code:`objfun`.
-* Bug fixes: default parameters for reduced initialization cost regime, returning correct value from safety steps, retrieving dependencies during installation.
+* Bug fixes: default parameters for reduced initialization cost regime, returning correct value if exiting from within a safety step, retrieving dependencies during installation.
+
+Version 1.1 (16 Jan 2019)
+-------------------------
+* Use different default reduced initialization cost method for inverse problems to ensure whole space is searched correctly.
+* Bug fixes: default trust region radius when scaling feasible region, exit correctly when no Jacobian returned, handling overflow at initial value