Bugfix for 1-sided bounds, deterministic initialization, faster subproblem solves, faster linear algebra in Hessian and interpolation problems. This closes #3 and closes #5

Author: lindonroberts

dfols/controller.py

@@ -326,22 +326,25 @@ def get_new_direction_for_growing(self, step_length):
 
     def trust_region_step(self):
         # Build model for full least squares objectives
-        gopt, hq = self.model.build_full_model()
-        d, gnew, crvmin = trsbox(self.model.xopt(), gopt, hq, self.model.sl, self.model.su, self.delta)
-        return d, gopt, hq, gnew, crvmin
+        gopt, H = self.model.build_full_model()
+        d, gnew, crvmin = trsbox(self.model.xopt(), gopt, H, self.model.sl, self.model.su, self.delta)
+        return d, gopt, H, gnew, crvmin
 
     def geometry_step(self, knew, adelt, number_of_samples, params):
         logging.debug("Running geometry-fixing step")
         try:
             c, g = self.model.lagrange_gradient(knew)
             # c = 1.0 if knew == self.model.kopt else 0.0  # based at xopt, just like d
             # Solve problem: bounds are sl <= xnew <= su, and ||xnew-xopt|| <= adelt
-            xnew = trsbox_geometry(self.model.xopt(), c, g, self.model.sl, self.model.su, adelt)
+            logging.debug("xopt = %s" % str(self.model.xopt()))
+            logging.debug("sl = %s" % str(self.model.sl))
+            logging.debug("su = %s" % str(self.model.su))
+            xnew = trsbox_geometry(self.model.xopt(), c, g, np.minimum(self.model.sl, 0.0), np.maximum(self.model.su, 0.0), adelt)
         except LA.LinAlgError:
             exit_info = ExitInformation(EXIT_LINALG_ERROR, "Singular matrix encountered in geometry step")
             return exit_info  # didn't fix geometry - return & quit
 
-        gopt, hq = self.model.build_full_model()  # save here, to calculate predicted value from geometry step
+        gopt, H = self.model.build_full_model()  # save here, to calculate predicted value from geometry step
         fopt = self.model.fopt()  # again, evaluate now, before model.change_point()
         d = xnew - self.model.xopt()
         x = self.model.as_absolute_coordinates(xnew)
@@ -362,8 +365,8 @@ def geometry_step(self, knew, adelt, number_of_samples, params):
         # Estimate actual reduction to add to diffs vector
         f = sumsq(np.mean(rvec_list[:num_samples_run, :], axis=0))  # estimate actual objective value
 
-        # pred_reduction = - calculate_model_value(gopt, hq, d)
-        pred_reduction = - model_value(gopt, hq, d)
+        # pred_reduction = - calculate_model_value(gopt, H, d)
+        pred_reduction = - model_value(gopt, H, d)
         actual_reduction = fopt - f
         self.diffs = [abs(pred_reduction - actual_reduction), self.diffs[0], self.diffs[1]]
         return None  # exit_info = None
@@ -424,16 +427,19 @@ def choose_point_to_replace(self, d, skip_kopt=True):
         knew = None  # may knew never be set here?
         exit_info = None
 
+        try:
+            cs, gs = self.model.lagrange_gradient(k=None)  # find all Lagrange polynomials for k in range(self.model.npt())
+        except LA.LinAlgError:
+            exit_info = ExitInformation(EXIT_LINALG_ERROR, "Singular matrix when choosing point to replace")
+            return knew, exit_info
+
         for k in range(self.model.npt()):
             if skip_kopt and k == self.model.kopt:
                 continue  # skip this k
 
-            # Build Lagrange polynomial
-            try:
-                c, g = self.model.lagrange_gradient(k)
-            except LA.LinAlgError:
-                exit_info = ExitInformation(EXIT_LINALG_ERROR, "Singular matrix when choosing point to replace")
-                break  # end & quit
+            # Extract Lagrange polynomial (based at xopt)
+            c = cs[k]
+            g = gs[:, k]
 
             den = c + np.dot(g, d)
 
@@ -445,9 +451,9 @@ def choose_point_to_replace(self, d, skip_kopt=True):
 
         return knew, exit_info
 
-    def done_with_current_rho(self, xnew, gnew, crvmin, hq, current_iter):
+    def done_with_current_rho(self, xnew, gnew, crvmin, H, current_iter):
         # (xnew, gnew, crvmin) come from trust region step
-        # hq is Hessian of model for the full objective
+        # H is Hessian of model for the full objective
 
         # Wait at least 3 iterations between reductions of rho
         if current_iter <= self.last_successful_iter + 2:
@@ -466,7 +472,7 @@ def done_with_current_rho(self, xnew, gnew, crvmin, hq, current_iter):
             if xnew[j] == self.model.su[j]:
                 bdtest = -gnew[j]
             if bdtest < bdtol:
-                curv = hq.get_element(j, j)  # curv = Hessian(j, j)
+                curv = H[j,j]
                 bdtest += 0.5 * curv * self.rho
                 if bdtest < bdtol:
                     return False
@@ -489,10 +495,10 @@ def reduce_rho(self, current_iter, params):
         self.last_successful_iter = current_iter  # reset successful iteration check
         return
 
-    def calculate_ratio(self, current_iter, rvec_list, d, gopt, hq):
+    def calculate_ratio(self, current_iter, rvec_list, d, gopt, H):
         exit_info = None
         f = sumsq(np.mean(rvec_list, axis=0))  # estimate actual objective value
-        pred_reduction = - model_value(gopt, hq, d)
+        pred_reduction = - model_value(gopt, H, d)
         actual_reduction = self.model.fopt() - f
         self.diffs = [abs(actual_reduction - pred_reduction), self.diffs[0], self.diffs[1]]
         if min(sqrt(sumsq(d)), self.delta) > self.rho:  # if ||d|| >= rho, successful!

dfols/model.py

@@ -35,7 +35,6 @@
 import numpy as np
 import scipy.linalg as LA
 
-from .hessian import Hessian
 from .trust_region import trsbox_geometry
 from .util import sumsq
 
@@ -119,7 +118,7 @@ def fopt(self):
     def xpt(self, k, abs_coordinates=False):
         assert 0 <= k < self.npt(), "Invalid index %g" % k
         if not abs_coordinates:
-            return self.points[k, :].copy()
+            return np.minimum(np.maximum(self.sl, self.points[k, :].copy()), self.su)
         else:
             # Apply bounds and convert back to absolute coordinates
             return self.xbase + np.minimum(np.maximum(self.sl, self.points[k, :]), self.su)
@@ -285,20 +284,23 @@ def factorise_geom_system(self):
         return
 
     def solve_geom_system(self, rhs):
+        # To do preconditioning below, we will need to scale each column of A elementwise by the entries of some vector
+        col_scale = lambda A, scale: (A.T*scale).T  # Uses the trick that A*x scales the 0th column of A by x[0], etc.
+
         if self.factorisation_current:
             if self.qr_of_transpose:
                 # Growing case: solve underdetermined system W*x=rhs with W.T = Q*R
                 # Golub & Van Loan (3rd edn), Algorithm 5.7.2
-                Rb = LA.solve_triangular(self.R, rhs * self.left_scaling, trans='T')  # R.T \ rhs
-                return np.dot(self.Q, Rb) * self.right_scaling  # minimal norm solution
+                Rb = LA.solve_triangular(self.R, col_scale(rhs, self.left_scaling), trans='T')  # R.T \ rhs
+                return col_scale(np.dot(self.Q, Rb), self.right_scaling)  # minimal norm solution
             else:
                 # Normal case: solve overdetermined system W*x=rhs with W=Q*R
-                Qb = np.dot(self.Q.T, rhs * self.left_scaling)
-                return LA.solve_triangular(self.R, Qb) * self.right_scaling
+                Qb = np.dot(self.Q.T, col_scale(rhs, self.left_scaling))
+                return col_scale(LA.solve_triangular(self.R, Qb), self.right_scaling)
         else:
             logging.warning("model.solve_geom_system not using factorisation")
             W, left_scaling, right_scaling = self.interpolation_matrix()
-            return LA.lstsq(W, rhs * left_scaling)[0] * right_scaling
+            return col_scale(LA.lstsq(W, col_scale(rhs * left_scaling))[0], right_scaling)
 
     def interpolate_mini_models_svd(self, verbose=False, make_full_rank=False, min_sing_val=1e-6, sing_val_frac=1.0, max_jac_cond=1e8,
                                     get_chg_J=False):
@@ -314,31 +316,26 @@ def interpolate_mini_models_svd(self, verbose=False, make_full_rank=False, min_s
         norm_J_error = 0.0
         linalg_resid = 0.0
 
-
         if make_full_rank:
             # Remove old full-rank components of Jacobian
             Y = self.xpt_directions(include_kopt=False).T
             Qy, Ry = LA.qr(Y, mode='full')  # Qy is (n,n), Ry is (n,npt-1)=(n,p)
             Qhat = Qy[:, :Y.shape[1]]
             self.model_jac = np.dot(self.model_jac, np.dot(Qhat, Qhat.T))
-        for m1 in range(self.m()):
-            g_old = self.model_jac[m1, :].copy()
-            rhs = self.fval_v[fval_row_idx, m1]  # length (npt)
-            try:
-                dg = self.solve_geom_system(rhs)
-            except LA.LinAlgError:
-                return False, None, None, None, None  # flag error
-            except ValueError:
-                return False, None, None, None, None  # flag error (e.g. inf or NaN encountered)
-
-            if not np.all(np.isfinite(dg)):  # another check for inf or NaN
-                return False, None, None, None, None  # flag error
 
-            self.model_jac[m1, :] = dg[1:]
-            self.model_const[m1] = dg[0] - np.dot(self.model_jac[m1, :], xopt)  # shift base to xbase
-            if verbose or get_chg_J:
-                norm_J_error += sumsq(dg[1:] - g_old)
-                linalg_resid += sumsq(np.dot(W, dg) - rhs)
+        rhs = self.fval_v[fval_row_idx, :]  # size npt * m
+        try:
+            dg = self.solve_geom_system(rhs)  # size (n+1)*m
+        except LA.LinAlgError:
+            return False, None, None, None, None  # flag error
+        except ValueError:
+            return False, None, None, None, None  # flag error (e.g. inf or NaN encountered)
+        J_old = self.model_jac.copy()
+        self.model_jac = dg[1:,:].T
+        self.model_const = dg[0,:] - np.dot(self.model_jac, xopt)  # shift base to xbase
+        if verbose or get_chg_J:
+            norm_J_error = np.linalg.norm(self.model_jac - J_old, ord='fro')**2
+            linalg_resid = np.linalg.norm(W.dot(dg) - rhs)**2
 
         if make_full_rank:
             try:
@@ -368,20 +365,29 @@ def build_full_model(self):
 
         # Apply scaling based on convention for objective - this code uses sumsq(rvec) not 0.5*sumsq(rvec)
         g = 2.0 * np.dot(J.T, r)  # n-vector
-        hess = Hessian(self.n(), vals=2.0 * np.dot(J.T, J))
-        return g, hess
+        H = 2.0 * np.dot(J.T, J)
+        return g, H
 
-    def lagrange_gradient(self, k, factorise_first=True):
-        assert 0 <= k < self.npt(), "Invalid index %g" % k
+    def lagrange_gradient(self, k=None, factorise_first=True):
         if factorise_first:
             self.factorise_geom_system()
 
-        rhs = np.zeros((self.npt(),))
-        rhs[k] = 1.0
+        if k is not None:
+            assert 0 <= k < self.npt(), "Invalid index %g" % k
+            rhs = np.zeros((self.npt(),))
+            rhs[k] = 1.0
+        else:
+            rhs = np.eye(self.npt())  # find all Lagrange polynomials
         soln = self.solve_geom_system(rhs)
-        c = soln[0]
-        g = soln[1:]
-        return c, g  # constant, gradient [all based at xopt]
+
+        if k is not None:
+            c = soln[0]
+            g = soln[1:]
+            return c, g  # constant, gradient [all based at xopt]
+        else:
+            cs = soln[0, :]
+            gs = soln[1:, :]
+            return cs, gs  # constant terms in each entry and gradient terms in each col [all based at xopt]
 
     def poisedness_constant(self, delta, xbase=None, xbase_in_abs_coords=True):
         # Calculate the poisedness constant of the current interpolation set in B(xbase, delta)
@@ -391,8 +397,13 @@ def poisedness_constant(self, delta, xbase=None, xbase_in_abs_coords=True):
             xbase = self.xopt()
         elif xbase_in_abs_coords:
             xbase = xbase - self.xbase  # shift to correct position
+        # Calculate all Lagrange polynomials at once
+        self.factorise_geom_system()
+        rhs = np.eye(self.npt())  # values to interpolate
+        soln = self.solve_geom_system(rhs)
         for k in range(self.npt()):
-            c, g = self.lagrange_gradient(k, factorise_first=True)
+            # Extract Lagrange poly from soln matrix (based at xopt)
+            c = soln[0,k]; g = soln[1:, k]
             newc = c + np.dot(g, xbase - self.xopt())  # based at xbase
             # Solve problem: bounds are sl <= x <= su, and ||x-xopt|| <= delta
             xmax = trsbox_geometry(xbase, newc, g, self.sl, self.su, delta)

dfols/params.py

@@ -39,7 +39,7 @@ def __init__(self, n, npt, maxfun, objfun_has_noise=False):
         self.params["general.safety_step_thresh"] = 0.5  # safety step called if ||d|| <= thresh * rho
         self.params["general.check_objfun_for_overflow"] = True
         # Initialisation
-        self.params["init.random_initial_directions"] = True
+        self.params["init.random_initial_directions"] = True if npt > (n+1)*(n+2)//2 else False  # only if needed
         self.params["init.run_in_parallel"] = False  # only available for random directions at the moment
         self.params["init.random_directions_make_orthogonal"] = True  # although random > orthogonal, avoid for init
         # Interpolation

dfols/solver.py

@@ -256,7 +256,7 @@ def solve_main(objfun, x0, args, xl, xu, npt, rhobeg, rhoend, maxfun, nruns_so_f
 
 
         # Trust region step
-        d, gopt, hq, gnew, crvmin = control.trust_region_step()
+        d, gopt, H, gnew, crvmin = control.trust_region_step()
         logging.debug("Trust region step is d = " + str(d))
         xnew = control.model.xopt() + d
         dnorm = min(LA.norm(d), control.delta)
@@ -353,7 +353,7 @@ def solve_main(objfun, x0, args, xl, xu, npt, rhobeg, rhoend, maxfun, nruns_so_f
                 diagnostic_info.update_iter_type(ITER_SAFETY)
                 diagnostic_info.update_slow_iter(-1)
 
-            if not control.done_with_current_rho(xnew, gnew, crvmin, hq, current_iter):
+            if not control.done_with_current_rho(xnew, gnew, crvmin, H, current_iter):
                 distsq = (10.0 * control.rho) ** 2
                 number_of_samples = max(nsamples(control.delta, control.rho, current_iter, nruns_so_far), 1)
                 update_delta = True  # we do reduce delta for safety steps
@@ -497,7 +497,7 @@ def solve_main(objfun, x0, args, xl, xu, npt, rhobeg, rhoend, maxfun, nruns_so_f
                 break  # quit
 
             # Estimate f in order to compute 'actual reduction'
-            ratio, exit_info = control.calculate_ratio(current_iter, rvec_list[:num_samples_run, :], d, gopt, hq)
+            ratio, exit_info = control.calculate_ratio(current_iter, rvec_list[:num_samples_run, :], d, gopt, H)
             if exit_info is not None:
                 if exit_info.able_to_do_restart() and params("restarts.use_restarts") and params(
                         "restarts.use_soft_restarts"):
@@ -824,11 +824,14 @@ def solve(objfun, x0, args=(), bounds=None, npt=None, rhobeg=None, rhoend=1e-8,
     if bounds is None:
         xl = None
         xu = None
-        scaling_within_bounds = False
     else:
         assert len(bounds) == 2, "bounds must be a 2-tuple of (lower, upper), where both are arrays of size(x0)"
-        xl = bounds[0].astype(np.float)
-        xu = bounds[1].astype(np.float)
+        xl = bounds[0].astype(np.float) if bounds[0] is not None else None
+        xu = bounds[1].astype(np.float) if bounds[1] is not None else None
+
+    if (xl is None or xu is None) and scaling_within_bounds:
+        scaling_within_bounds = False
+        warnings.warn("Ignoring scaling_within_bounds=True for unconstrained problem/1-sided bounds", RuntimeWarning)
 
     if xl is None:
         xl = -1e20 * np.ones((n,))  # unconstrained

dfols/tests/test_model.py

@@ -280,11 +280,11 @@ def runTest(self):
                                       rosenbrock(x2)-rosenbrock(model.xbase), thresh=1e-10), 'Wrong x2 (no constant)')
         self.assertTrue(array_compare(model.model_value(x2 - model.xbase, d_based_at_xopt=False, with_const_term=False),
                                 rosenbrock(x2) - rosenbrock(model.xbase), thresh=1e-10), 'Wrong x2 (no constant v2)')
-        g, hess = model.build_full_model()
+        g, H = model.build_full_model()
         r = rosenbrock(x1)
         J = model.model_jac
         self.assertTrue(array_compare(g, 2.0*np.dot(J.T, r), thresh=1e-10), 'Bad gradient')
-        self.assertTrue(array_compare(hess.as_full(), 2.0*np.dot(J.T, J)), 'Bad Hessian')
+        self.assertTrue(array_compare(H, 2.0*np.dot(J.T, J)), 'Bad Hessian')
 
 
 class TestGeomSystem(unittest.TestCase):
@@ -426,7 +426,7 @@ def runTest(self):
         self.assertTrue(array_compare(A_for_interp, A_after_scaling), 'Interp matrix 1')
         # For reference: model based around model.xbase
         interp_ok, interp_error, norm_J_error, linalg_resid, ls_interp_cond_num = model.interpolate_mini_models_svd()
-        J_true = np.linalg.lstsq(A_for_interp, b)[0]
+        J_true = np.linalg.lstsq(A_for_interp, b, rcond=None)[0]
         self.assertTrue(interp_ok, 'Interpolation failed')
         # print(model.model_const, model.model_jac)
         # print(J_true[0] - np.dot(J_true[1:], model.xopt()), J_true[1:])
@@ -468,7 +468,7 @@ def runTest(self):
         # For reference: model based around model.xbase
         interp_ok, interp_error, norm_J_error, linalg_resid, ls_interp_cond_num = model.interpolate_mini_models_svd()
         self.assertTrue(interp_ok, 'Interpolation failed')
-        J_true = np.linalg.lstsq(A_for_interp, b)[0]
+        J_true = np.linalg.lstsq(A_for_interp, b, rcond=None)[0]
         g = J_true[1:]
         c = J_true[0] - np.dot(g, model.xopt())  # centred at xbase
         self.assertTrue(array_compare(model.model_const, np.array([c])), 'Wrong constant term')

dfols/tests/test_params.py

@@ -34,9 +34,9 @@ def runTest(self):
         npt = n + 1
         maxfun = 50 * (n + 1)
         p = ParameterList(n, npt, maxfun)
-        self.assertTrue(p("init.random_initial_directions"), 'Bad init dirns/access')
-        p("init.random_initial_directions", False)  # set to False
         self.assertFalse(p("init.random_initial_directions"), 'Bad init dirns/access')
+        p("init.random_initial_directions", True)  # set to True
+        self.assertTrue(p("init.random_initial_directions"), 'Bad init dirns/access')
 
 
 class TestFail(unittest.TestCase):

dfols/tests/test_solver.py

@@ -87,7 +87,7 @@ def runTest(self):
         print(soln.x)
         self.assertTrue(array_compare(soln.x, xmin, thresh=1e-2), "Wrong xmin")
         self.assertTrue(abs(soln.f - fmin) < 1e-4, "Wrong fmin")
-        self.assertTrue(array_compare(soln.jacobian, jacmin, thresh=1e-3), "Wrong jacobian")
+        self.assertTrue(array_compare(soln.jacobian, jacmin, thresh=2e-2), "Wrong jacobian")
 
 
 class TestLinear(unittest.TestCase):
@@ -98,7 +98,7 @@ def runTest(self):
         A = np.random.rand(m, n)
         b = np.random.rand(m)
         objfun = lambda x: np.dot(A, x) - b
-        xmin = np.linalg.lstsq(A, b)[0]
+        xmin = np.linalg.lstsq(A, b, rcond=None)[0]
         fmin = np.dot(objfun(xmin), objfun(xmin))
         x0 = np.zeros((n,))
         soln = dfols.solve(objfun, x0)