Modular model building and model updating

pounders/py/formquad.py

@@ -1,13 +1,14 @@
 import numpy as np
 import scipy.linalg
 
+from .bmpts import bmpts
 from .phi2eval import phi2eval
 
 # from .flipFirstRow import flipFirstRow
 # from .flipSignQ import flipSignQ
 
 
-def formquad(X, F, delta, xk_in, np_max, Pars, vf):
+def formquad(X, F, delta, xk_in, np_max, Pars, vf, H_flag=False, Old_H=[]):
     """
     formquad(X, F, delta, xk_in, np_max, Pars, vf) -> [Mdir, np, valid, G, H, Mind]
     Computes the parameters for m quadratics
@@ -22,13 +23,17 @@ def formquad(X, F, delta, xk_in, np_max, Pars, vf):
     X       [dbl] [nf-by-n] Locations of evaluated points
     F       [dbl] [nf-by-m] Function values of evaluated points
     delta   [dbl] Positive trust region radius
-    xk_in    [int] Index in (X and F) of the current center
-    np_max   [int] Max # interpolation points (>=n+1) (.5*(n+1)*(n+2))
+    xk_in   [int] Index in (X and F) of the current center
+    np_max  [int] Max # interpolation points (>=n+1) (.5*(n+1)*(n+2))
     Pars[0] [dbl] delta multiplier for checking validity
     Pars[1] [dbl] delta multiplier for all interpolation points
     Pars[2] [dbl] Pivot threshold for validity
     Pars[3] [dbl] Pivot threshold for additional points (.001)
     vf      [log] Flag indicating you just want to check model validity
+    H_flag  [int] Flag indicating type of Hessians to return
+                  0: Minimum change in Frobenius norm Hessians
+                  1: Minimum-norm Hessians
+
     --OUTPUTS----------------------------------------------------------------
     Mdir    [dbl] [(n-np+1)-by-n]  Unit directions to improve model
     np      [int] Number of interpolation points (=length(Mind))
@@ -93,7 +98,7 @@ def formquad(X, F, delta, xk_in, np_max, Pars, vf):
         elif aff == 1 and mp < n:  # Need to evaluate more points, then recall
             Mdir = Q[:, mp:n].T  # Output Geometry directions
             G = np.empty(shape=(0, 0))
-            H = np.empty(shape=(0, 0))
+            # H = np.empty(shape=(0, 0))
             return [Mdir, mp, valid, G, H, Mind]
         elif aff == 0:  # Output model-improving directions
             Mdir = Q[:, mp:n].T  # Will be empty if mp=n
@@ -166,4 +171,37 @@ def formquad(X, F, delta, xk_in, np_max, Pars, vf):
     H = H / (delta**2)
     G = G / delta
 
+    if H_flag:
+        H = Old_H + H
+
     return [Mdir, mp, valid, G, H, Mind]
+
+
+def formquad_model_improvement(x_k, Cres, Gres, Hres, Mdir, mp, Low, Upp, delta, Model, combinemodels):
+    n = len(x_k)
+
+    # Update for modelimp; Cres unchanged b/c xk_in unchanged
+    G, H = combinemodels(Cres, Gres, Hres)
+    # Evaluate model-improving points to pick best one
+    # May eventually want to normalize Mdir first for infty norm
+    # Plus directions
+    [Mdir1, mp1] = bmpts(x_k, Mdir[0 : n - mp, :], Low, Upp, delta, Model["Par"][2])
+    Res = np.zeros((n - mp1, 1))
+    for i in range(n - mp1):
+        D = Mdir1[i, :]
+        Res[i, 0] = D @ (G + 0.5 * H @ D.T)
+    b = np.argmin(Res[: n - mp1, 0:1])
+    a1 = np.min(Res[: n - mp1, 0:1])
+    Xsp = Mdir1[b, :]
+    # Minus directions
+    [Mdir1, mp2] = bmpts(x_k, -Mdir[0 : n - mp, :], Low, Upp, delta, Model["Par"][2])
+    Res = np.zeros((n - mp2, 1))
+    for i in range(n - mp2):
+        D = Mdir1[i, :]
+        Res[i, 0] = D @ (G + 0.5 * H @ D.T)
+    b = np.argmin(Res[: n - mp2, 0:1])
+    a2 = np.min(Res[: n - mp2, 0:1])
+    if a2 < a1:
+        Xsp = Mdir1[b, :]
+
+    return Xsp

pounders/py/pounders.py

@@ -5,7 +5,7 @@
 from .bmpts import bmpts
 from .bqmin import bqmin
 from .checkinputss import checkinputss
-from .formquad import formquad
+from .formquad import formquad, formquad_model_improvement
 from .prepare_outputs_before_return import prepare_outputs_before_return
 
 
@@ -109,9 +109,12 @@ def pounders(Ffun, X_0, n, nf_max, g_tol, delta_0, m, Low, Upp, Prior=None, Opti
 
     if Model is None:
         Model = {}
+        Model["H_flag"] = True
         Model["Par"] = _default_model_par_values(n)
         Model["np_max"] = _default_model_np_max(n)
     else:
+        if "H_flag" not in Model:
+            Model["H_flag"] = True
         if "Par" not in Model:
             Model["Par"] = _default_model_par_values(n)
         if "np_max" not in Model:
@@ -190,10 +193,13 @@ def pounders(Ffun, X_0, n, nf_max, g_tol, delta_0, m, Low, Upp, Prior=None, Opti
     Hres = np.zeros((n, n, m))
     ng = np.nan  # Needed for early termination, e.g., if a model is never built
     while nf + 1 < nf_max:
-        #  1a. Compute the interpolation set.
+        #  1a. Compute displacements
         D = X[: nf + 1] - X[xk_in]
         Res[: nf + 1, :] = (F[: nf + 1, :] - Cres) - np.diagonal(0.5 * D @ (np.tensordot(D, Hres, axes=1))).T
-        [Mdir, mp, valid, Gres, Hresdel, Mind] = formquad(X[0 : nf + 1, :], Res[0 : nf + 1, :], delta, xk_in, Model["np_max"], Model["Par"], 0)
+
+        # Request directions for points to be evaluated
+        [Mdir, mp, valid, Gres, Hres, Mind] = formquad(X[0 : nf + 1, :], Res[0 : nf + 1, :], delta, xk_in, Model["np_max"], Model["Par"], False, Model["H_flag"], Hres)
+
         if mp < n:
             [Mdir, mp] = bmpts(X[xk_in], Mdir[0 : n - mp, :], Low, Upp, delta, Model["Par"][2])
             for i in range(int(min(n - mp, nf_max - (nf + 1)))):
@@ -210,14 +216,12 @@ def pounders(Ffun, X_0, n, nf_max, g_tol, delta_0, m, Low, Upp, Prior=None, Opti
                 Res[nf, :] = (F[nf, :] - Cres) - 0.5 * D @ np.tensordot(D.T, Hres, 1)
             if nf + 1 >= nf_max:
                 break
-            [_, mp, valid, Gres, Hresdel, Mind] = formquad(X[0 : nf + 1, :], Res[0 : nf + 1, :], delta, xk_in, Model["np_max"], Model["Par"], False)
+
+            [_, mp, valid, Gres, Hres, Mind] = formquad(X[0 : nf + 1, :], Res[0 : nf + 1, :], delta, xk_in, Model["np_max"], Model["Par"], False, Model["H_flag"], Hres)
             if mp < n:
                 X, F, hF, flag = prepare_outputs_before_return(X, F, hF, nf, -5)
                 return X, F, hF, flag, xk_in
 
-        #  1b. Update the quadratic model
-        Cres = F[xk_in]
-        Hres = Hres + Hresdel
         c = hF[xk_in]
         G, H = combinemodels(Cres, Gres, Hres)
         ind_Lownotbinding = (X[xk_in] > Low) * (G.T > 0)
@@ -253,14 +257,14 @@ def pounders(Ffun, X_0, n, nf_max, g_tol, delta_0, m, Low, Upp, Prior=None, Opti
                     F[nf] = Ffun(X[nf])
                     if np.any(np.isnan(F[nf])):
                         X, F, hF, flag = prepare_outputs_before_return(X, F, hF, nf, -3)
-                        return X, F, flag, xk_in
+                        return X, F, hF, flag, xk_in
                     hF[nf] = hfun(F[nf])
                     if printf:
                         print("%4i   Critical point  %11.5e\n" % (nf, hF[nf]))
                 if nf + 1 >= nf_max:
                     break
                 # Recalculate gradient based on a MFN model
-                [_, _, valid, Gres, Hres, Mind] = formquad(X[: nf + 1, :], F[: nf + 1, :], delta, xk_in, Model["np_max"], Model["Par"], 0)
+                [_, _, valid, Gres, Hres, _] = formquad(X[: nf + 1, :], F[: nf + 1, :], delta, xk_in, Model["np_max"], Model["Par"], 0)
                 G, H = combinemodels(Cres, Gres, Hres)
                 ind_Lownotbinding = (X[xk_in] > Low) * (G.T > 0)
                 ind_Uppnotbinding = (X[xk_in] < Upp) * (G.T < 0)
@@ -336,34 +340,17 @@ def pounders(Ffun, X_0, n, nf_max, g_tol, delta_0, m, Low, Upp, Prior=None, Opti
         # 5. Evaluate a model-improving point if necessary
         if not valid and (nf + 1 < nf_max) and (rho < eta_1):  # Implies xk_in, delta unchanged
             # Need to check because model may be valid after Xsp evaluation
-            [Mdir, mp, valid, _, _, _] = formquad(X[: nf + 1, :], F[: nf + 1, :], delta, xk_in, Model["np_max"], Model["Par"], 1)
+            [Mdir, mp, valid, _, _, _] = formquad(X[: nf + 1, :], F[: nf + 1, :], delta, xk_in, Model["np_max"], Model["Par"], True)
+
             if not valid:  # ! One strategy for choosing model-improving point:
                 # Update model (exists because delta & xk_in unchanged)
                 D = X[: nf + 1] - X[xk_in]
                 Res[: nf + 1, :] = (F[: nf + 1, :] - Cres) - np.diagonal(0.5 * D @ (np.tensordot(D, Hres, axes=1))).T
-                [_, _, valid, Gres, Hresdel, Mind] = formquad(X[: nf + 1, :], Res[: nf + 1, :], delta, xk_in, Model["np_max"], Model["Par"], False)
-                Hres = Hres + Hresdel
-                # Update for modelimp; Cres unchanged b/c xk_in unchanged
-                G, H = combinemodels(Cres, Gres, Hres)
-                # Evaluate model-improving points to pick best one
-                # May eventually want to normalize Mdir first for infty norm
-                # Plus directions
-                [Mdir1, mp1] = bmpts(X[xk_in], Mdir[0 : n - mp, :], Low, Upp, delta, Model["Par"][2])
-                for i in range(n - mp1):
-                    D = Mdir1[i, :]
-                    Res[i, 0] = D @ (G + 0.5 * H @ D.T)
-                b = np.argmin(Res[: n - mp1, 0:1])
-                a1 = np.min(Res[: n - mp1, 0:1])
-                Xsp = Mdir1[b, :]
-                # Minus directions
-                [Mdir1, mp2] = bmpts(X[xk_in], -Mdir[0 : n - mp, :], Low, Upp, delta, Model["Par"][2])
-                for i in range(n - mp2):
-                    D = Mdir1[i, :]
-                    Res[i, 0] = D @ (G + 0.5 * H @ D.T)
-                b = np.argmin(Res[: n - mp2, 0:1])
-                a2 = np.min(Res[: n - mp2, 0:1])
-                if a2 < a1:
-                    Xsp = Mdir1[b, :]
+
+                [_, _, valid, Gres, Hres, Mind] = formquad(X[: nf + 1, :], Res[: nf + 1, :], delta, xk_in, Model["np_max"], Model["Par"], False, Model["H_flag"], Hres)
+
+                Xsp = formquad_model_improvement(X[xk_in], Cres, Gres, Hres, Mdir, mp, Low, Upp, delta, Model, combinemodels)
+
                 nf += 1
                 X[nf] = np.minimum(Upp, np.maximum(Low, X[xk_in] + Xsp))  # Temp safeguard
                 F[nf] = Ffun(X[nf])
@@ -383,6 +370,7 @@ def pounders(Ffun, X_0, n, nf_max, g_tol, delta_0, m, Low, Upp, Prior=None, Opti
                     # Don't actually use
                     for j in range(m):
                         Gres[:, j] = Gres[:, j] + Hres[:, :, j] @ D.T
+
     if printf:
         print("Number of function evals exceeded")
     flag = ng