Use np.ndarray for Hessians, not Hessian class (faster matrix-vector multiplication)

Author: lindonroberts

pybobyqa/controller.py

@@ -35,7 +35,6 @@
 import numpy as np
 import scipy.linalg as LA
 
-from .hessian import *
 from .model import *
 from .trust_region import *
 from .util import *
@@ -266,21 +265,21 @@ def initialise_random_directions(self, number_of_samples, num_directions, params
 
     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, hess = self.model.lagrange_polynomial(knew)  # based at xopt
+            c, g, H = self.model.lagrange_polynomial(knew)  # based at xopt
             # Solve problem: bounds are sl <= xnew <= su, and ||xnew-xopt|| <= adelt
-            xnew = trsbox_geometry(self.model.xopt(), c, g, hess, self.model.sl, self.model.su, adelt)
+            xnew = trsbox_geometry(self.model.xopt(), c, g, H, self.model.sl, self.model.su, 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)
@@ -300,8 +299,8 @@ def geometry_step(self, knew, adelt, number_of_samples, params):
         # Estimate actual reduction to add to diffs vector
         f = np.mean(f_list[:num_samples_run])  # 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
@@ -381,9 +380,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:
@@ -402,7 +401,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
@@ -425,10 +424,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, f_list, d, gopt, hq):
+    def calculate_ratio(self, current_iter, f_list, d, gopt, H):
         exit_info = None
         f = np.mean(f_list)  # 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!

pybobyqa/model.py

@@ -34,7 +34,7 @@
 import numpy as np
 import scipy.linalg as LA
 
-from .hessian import Hessian, to_upper_triangular_vector
+from .hessian import to_upper_triangular_vector
 from .trust_region import trsbox_geometry
 from .util import sumsq, model_value
 
@@ -74,7 +74,7 @@ def __init__(self, npt, x0, f0, xl, xu, f0_nsamples, n=None, abs_tol=-1e20, prec
         # Model information
         self.model_const = 0.0  # constant term for model m(s) = c + J*s
         self.model_grad = np.zeros((n,))  # Jacobian term for model m(s) = c + J*s
-        self.model_hess = Hessian(n)
+        self.model_hess = np.zeros((n,n))
 
         # Saved point (in absolute coordinates) - always check this value before quitting solver
         self.xsave = None
@@ -199,7 +199,7 @@ def shift_base(self, xbase_shift):
         self.factorisation_current = False
 
         # Update model (always centred on xbase)
-        Hx = self.model_hess.vec_mul(xbase_shift)
+        Hx = self.model_hess.dot(xbase_shift)
         self.model_const += np.dot(self.model_grad + 0.5*Hx, xbase_shift)
         self.model_grad += Hx
         return
@@ -209,7 +209,7 @@ def save_point(self, x, f, nsamples, x_in_abs_coords=True):
             self.xsave = x.copy() if x_in_abs_coords else self.as_absolute_coordinates(x)
             self.fsave = f
             self.gradsave = self.model_grad.copy()
-            self.hesssave = self.model_hess.as_full().copy()
+            self.hesssave = self.model_hess.copy()
             self.nsamples_save = nsamples
             return True
         else:
@@ -219,7 +219,7 @@ def get_final_results(self):
         # Return x and fval for optimal point (either from xsave+fsave or kopt)
         if self.fsave is None or self.fopt() <= self.fsave:  # optimal has changed since xsave+fsave were last set
             g, hess = self.build_full_model()  # model based at xopt
-            return self.xopt(abs_coordinates=True).copy(), self.fopt(), g, hess.as_full(), self.nsamples[self.kopt]
+            return self.xopt(abs_coordinates=True).copy(), self.fopt(), g, hess, self.nsamples[self.kopt]
         else:
             return self.xsave, self.fsave, self.gradsave, self.hesssave, self.nsamples_save
 
@@ -231,7 +231,7 @@ def model_value(self, d, d_based_at_xopt=True, with_const_term=False):
         # Model is always centred around xbase
         const = self.model_const if with_const_term else 0.0
         d_to_use = d + self.xopt() if d_based_at_xopt else d
-        Hd = self.model_hess.vec_mul(d_to_use)
+        Hd = self.model_hess.dot(d_to_use)
         return const + np.dot(self.model_grad + 0.5 * Hd, d_to_use)
 
     def interpolation_matrix(self):
@@ -281,7 +281,7 @@ def interpolate_model(self, verbose=False, min_chg_hess=True, get_norm_model_chg
                 # It's good to see which bits are needed for this specifically (here & 1 line below)
                 for t in range(self.npt()-1):
                     dx = self.xpt(fval_row_idx[t]) - self.xopt()
-                    rhs[t] = rhs[t] - 0.5 * np.dot(dx, self.model_hess.vec_mul(dx))  # include old Hessian
+                    rhs[t] = rhs[t] - 0.5 * np.dot(dx, self.model_hess.dot(dx))  # include old Hessian
 
         try:
             coeffs = self.solve_system(rhs)
@@ -296,7 +296,7 @@ def interpolate_model(self, verbose=False, min_chg_hess=True, get_norm_model_chg
         # Old gradient and Hessian (save so can compute changes later)
         if verbose or get_norm_model_chg:
             old_model_grad = self.model_grad.copy()
-            old_model_hess = self.model_hess.as_full()
+            old_model_hess = self.model_hess.copy()
         else:
             old_model_grad = None
             old_model_hess = None
@@ -305,24 +305,21 @@ def interpolate_model(self, verbose=False, min_chg_hess=True, get_norm_model_chg
         self.model_const = self.fopt()  # true in all cases
         if self.npt() == self.n() + 1:
             self.model_grad = coeffs.copy()
-            self.model_hess = Hessian(self.n())  # zeros
+            self.model_hess = np.zeros((self.n(), self.n()))
         elif self.npt() == (self.n() + 1) * (self.n() + 2) // 2:
             self.model_grad = coeffs[:self.n()]
-            self.model_hess = Hessian(self.n(), coeffs[self.n():])  # rest of coeffs are upper triangular part of Hess
+            self.model_hess = build_symmetric_matrix_from_vector(self.n(), coeffs[self.n():])  # rest of coeffs are upper triangular part of Hess
         else:
             self.model_grad = coeffs[self.npt()-1:]  # last n values
-            if min_chg_hess:
-                hess_full = self.model_hess.as_full()
-            else:
-                hess_full = np.zeros((self.n(), self.n()))
+            if not min_chg_hess:
+                self.model_hess = np.zeros((self.n(), self.n()))
             for i in range(self.npt()-1):
                 dx = self.xpt(fval_row_idx[i]) - self.xopt()
-                hess_full += coeffs[i] * np.outer(dx, dx)
-            self.model_hess = Hessian(self.n(), hess_full)
+                self.model_hess += coeffs[i] * np.outer(dx, dx)
 
         # Base model at xbase, not xopt (note negative signs)
         xopt = self.xopt()
-        Hx = self.model_hess.vec_mul(xopt)
+        Hx = self.model_hess.dot(xopt)
         self.model_const += np.dot(-self.model_grad + 0.5*Hx, xopt)
         self.model_grad += -Hx
 
@@ -331,7 +328,7 @@ def interpolate_model(self, verbose=False, min_chg_hess=True, get_norm_model_chg
         norm_chg_hess = 0.0
         if verbose or get_norm_model_chg:
             norm_chg_grad = LA.norm(self.model_grad - old_model_grad)
-            norm_chg_hess = LA.norm(self.model_hess.as_full() - old_model_hess, ord='fro')
+            norm_chg_hess = LA.norm(self.model_hess - old_model_hess, ord='fro')
         if verbose:
             for k in range(self.npt()):
                 f_pred = self.model_value(self.xpt(k), d_based_at_xopt=False, with_const_term=True)
@@ -342,7 +339,7 @@ def interpolate_model(self, verbose=False, min_chg_hess=True, get_norm_model_chg
 
     def build_full_model(self):
         # Make model centred around xopt
-        g = self.model_grad + self.model_hess.vec_mul(self.xopt())
+        g = self.model_grad + self.model_hess.dot(self.xopt())
         return g, self.model_hess
 
     def lagrange_polynomial(self, k, factorise_first=True):
@@ -382,21 +379,20 @@ def lagrange_polynomial(self, k, factorise_first=True):
         c = 1.0 if k==self.kopt else 0.0  # true in all cases
         if self.npt() == self.n() + 1:
             g = coeffs.copy()
-            hess = Hessian(self.n())  # zeros
+            H = np.zeros((self.n(), self.n()))
         elif self.npt() == (self.n() + 1) * (self.n() + 2) // 2:
             g = coeffs[:self.n()]
-            hess = Hessian(self.n(), coeffs[self.n():])  # rest of coeffs are upper triangular part of Hess
+            H = build_symmetric_matrix_from_vector(self.n(), coeffs[self.n():])  # rest of coeffs are upper triangular part of Hess
         else:
             g = coeffs[self.npt() - 1:]  # last n values
             fval_row_idx = np.delete(np.arange(self.npt()), self.kopt)  # indices of all rows except kopt
-            hess_full = np.zeros((self.n(), self.n()))
+            H = np.zeros((self.n(), self.n()))
             for i in range(self.npt() - 1):
                 dx = self.xpt(fval_row_idx[i]) - self.xopt()
-                hess_full += coeffs[i] * np.outer(dx, dx)
-            hess = Hessian(self.n(), hess_full)
+                H += coeffs[i] * np.outer(dx, dx)
 
         # (c, g, hess) currently based around xopt
-        return c, g, hess
+        return c, g, H
 
     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)
@@ -407,14 +403,14 @@ def poisedness_constant(self, delta, xbase=None, xbase_in_abs_coords=True):
         elif xbase_in_abs_coords:
             xbase = xbase - self.xbase  # shift to correct position
         for k in range(self.npt()):
-            c, g, hess = self.lagrange_polynomial(k, factorise_first=True)  # based at self.xopt()
+            c, g, H = self.lagrange_polynomial(k, factorise_first=True)  # based at self.xopt()
             # Switch base of poly from xopt to xbase, as required by trsbox_geometry
             base_chg = self.xopt() - xbase
-            Hx = hess.vec_mul(base_chg)
+            Hx = H.dot(base_chg)
             c += np.dot(-g + 0.5 * Hx, base_chg)
             g += -Hx
-            xmax = trsbox_geometry(xbase, c, g, hess, self.sl, self.su, delta)
-            lmax = abs(c + model_value(g, hess, xmax-xbase))  # evaluate Lagrange poly
+            xmax = trsbox_geometry(xbase, c, g, H, self.sl, self.su, delta)
+            lmax = abs(c + model_value(g, H, xmax-xbase))  # evaluate Lagrange poly
             if overall_max is None or lmax > overall_max:
                 overall_max = lmax
         return overall_max
@@ -456,3 +452,14 @@ def build_interpolation_matrix(Y, approx_delta=1.0):
         right_scaling[p:] = approx_delta
     return A, left_scaling, right_scaling
 
+
+def build_symmetric_matrix_from_vector(n, entries):
+    assert entries.shape == (n*(n+1)//2,), "Entries vector has wrong size, got %g, expect %g (for n=%g)" % (len(entries), n*(n+1)//2, n)
+    A = np.zeros((n, n))
+    ih = -1
+    for j in range(n):  # j = 0, ..., n-1
+        for i in range(j + 1):  # i = 0, ..., j
+            ih += 1
+            A[i, j] = entries[ih]
+            A[j, i] = entries[ih]
+    return A

pybobyqa/solver.py

@@ -222,7 +222,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)
@@ -242,7 +242,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
@@ -378,7 +378,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, f_list[:num_samples_run], d, gopt, hq)
+            ratio, exit_info = control.calculate_ratio(current_iter, f_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"):