Merge pull request #22 from numericalalgorithmsgroup/dev_regularized

Nonsmooth regularizers

Author: lindonroberts

README.rst

@@ -34,13 +34,15 @@ See manual.pdf or `here <https://numericalalgorithmsgroup.github.io/dfols/>`_.
 
 Citation
 --------
-If you use DFO-LS in a paper, please cite:
+The development of DFO-LS is outlined over several publications:
 
-Cartis, C., Fiala, J., Marteau, B. and Roberts, L., `Improving the Flexibility and Robustness of Model-Based Derivative-Free Optimization Solvers <https://doi.org/10.1145/3338517>`_, *ACM Transactions on Mathematical Software*, 45:3 (2019), pp. 32:1-32:41.
+1. C. Cartis, J. Fiala, B. Marteau and L. Roberts, `Improving the Flexibility and Robustness of Model-Based Derivative-Free Optimization Solvers <https://doi.org/10.1145/3338517>`_, *ACM Transactions on Mathematical Software*, 45:3 (2019), pp. 32:1-32:41 [`preprint <https://arxiv.org/abs/1804.00154>`_] . 
+2. M. Hough, and L. Roberts, `Model-Based Derivative-Free Methods for Convex-Constrained Optimization <https://doi.org/10.1137/21M1460971>`_, *SIAM Journal on Optimization*, 21:4 (2022), pp. 2552-2579 [`preprint <https://arxiv.org/abs/2111.05443>`_].
+3. Y. Liu, K. H. Lam and L. Roberts, `Black-box Optimization Algorithms for Regularized Least-squares Problems <http://arxiv.org/abs/2407.14915>`_, *arXiv preprint arXiv:arXiv:2407.14915*, 2024.
 
-If you use DFO-LS for problems with constraints, including bound constraints, please also cite:
-
-Hough, M. and Roberts, L., `Model-Based Derivative-Free Methods for Convex-Constrained Optimization <https://doi.org/10.1137/21M1460971>`_, *SIAM Journal on Optimization*, 21:4 (2022), pp. 2552-2579.
+If you use DFO-LS in a paper, please cite [1]. 
+If your problem has constraints, including bound constraints, please cite [1,2].
+If your problem includes a regularizer, please cite [1,3].
 
 Requirements
 ------------
@@ -70,27 +72,13 @@ For easy installation, use `pip <http://www.pip-installer.org/>`_ as root:
 
 .. code-block:: bash
 
-    $ [sudo] pip install DFO-LS
-
-or alternatively *easy_install*:
-
-.. code-block:: bash
-
-    $ [sudo] easy_install DFO-LS
-
-If you do not have root privileges or you want to install DFO-LS for your private use, you can use:
-
-.. code-block:: bash
-
-    $ pip install --user DFO-LS
-
-which will install DFO-LS in your home directory.
+    $ pip install DFO-LS
 
 Note that if an older install of DFO-LS is present on your system you can use:
 
 .. code-block:: bash
 
-    $ [sudo] pip install --upgrade DFO-LS
+    $ pip install --upgrade DFO-LS
 
 to upgrade DFO-LS to the latest version.
 
@@ -107,22 +95,14 @@ DFO-LS is written in pure Python and requires no compilation. It can be installe
 
  .. code-block:: bash
 
-    $ [sudo] pip install .
-
-If you do not have root privileges or you want to install DFO-LS for your private use, you can use:
-
- .. code-block:: bash
-
-    $ pip install --user .
-
-instead.
+    $ pip install .
 
 To upgrade DFO-LS to the latest version, navigate to the top-level directory (i.e. the one containing :code:`pyproject.toml`) and rerun the installation using :code:`pip`, as above:
 
  .. code-block:: bash
 
     $ git pull
-    $ [sudo] pip install .  # with admin privileges
+    $ pip install .
 
 Testing
 -------
@@ -145,7 +125,7 @@ If DFO-LS was installed using *pip* you can uninstall as follows:
 
  .. code-block:: bash
 
-    $ [sudo] pip uninstall DFO-LS
+    $ pip uninstall DFO-LS
 
 If DFO-LS was installed manually you have to remove the installed files by hand (located in your python site-packages directory).
 

dfols/__init__.py

@@ -39,7 +39,7 @@
 from __future__ import absolute_import, division, print_function, unicode_literals
 
 # DFO-LS version
-__version__ = '1.4.1'
+__version__ = '1.5.0'
 
 # Main solver & exit flags
 from .solver import *

dfols/controller.py

@@ -36,16 +36,16 @@
 import scipy.linalg as LA
 
 from .trust_region import trsbox_geometry
-from .util import sumsq, dykstra
+from .util import sumsq, dykstra, remove_scaling
 
 __all__ = ['Model']
 
 module_logger = logging.getLogger(__name__) 
 
 
 class Model(object):
-    def __init__(self, npt, x0, r0, xl, xu, projections, r0_nsamples, n=None, m=None, abs_tol=1e-12, rel_tol=1e-20, precondition=True,
-                 do_logging=True):
+    def __init__(self, npt, x0, r0, xl, xu, projections, r0_nsamples, h=None, argsh=(), n=None, m=None, abs_tol=1e-12, rel_tol=1e-20, precondition=True,
+                 do_logging=True, scaling_changes=None):
         if n is None:
             n = len(x0)
         if m is None:
@@ -56,11 +56,15 @@ def __init__(self, npt, x0, r0, xl, xu, projections, r0_nsamples, n=None, m=None
         assert xu.shape == (n,), "xu has wrong shape (got %s, expect (%g,))" % (str(xu.shape), n)
         assert r0.shape == (m,), "r0 has wrong shape (got %s, expect (%g,))" % (str(r0.shape), m)
         self.do_logging = do_logging
+        self.scaling_changes = scaling_changes
         self.dim = n
         self.resid_dim = m
         self.num_pts = npt
         self.npt_so_far = 1  # number of points added so far (with function values)
 
+        self.h = h
+        self.argsh = argsh
+
         # Initialise to blank some useful stuff
         # Interpolation points
         self.xbase = x0.copy()
@@ -72,12 +76,15 @@ def __init__(self, npt, x0, r0, xl, xu, projections, r0_nsamples, n=None, m=None
         # Function values
         self.fval_v = np.inf * np.ones((npt, m))  # residuals for each xpt
         self.fval_v[0, :] = r0.copy()
-        self.fval = np.inf * np.ones((npt, ))  # overall objective value for each xpt
-        self.fval[0] = sumsq(r0)
+        
+        self.objval = np.inf * np.ones((npt, ))  # overall objective value for each xpt
+        self.objval[0] = sumsq(r0)
+        if h is not None:
+            self.objval[0] += h(remove_scaling(x0, self.scaling_changes), *argsh)
         self.kopt = 0  # index of current iterate (should be best value so far)
         self.nsamples = np.zeros((npt,), dtype=int)  # number of samples used to evaluate objective at each point
         self.nsamples[0] = r0_nsamples
-        self.fbeg = self.fval[0]  # f(x0), saved to check for sufficient reduction
+        self.objbeg = self.objval[0]  # f(x0), saved to check for sufficient reduction
 
         # Termination criteria
         self.abs_tol = abs_tol
@@ -90,7 +97,7 @@ def __init__(self, npt, x0, r0, xl, xu, projections, r0_nsamples, n=None, m=None
         # Saved point (in absolute coordinates) - always check this value before quitting solver
         self.xsave = None
         self.rsave = None
-        self.fsave = None
+        self.objsave = None
         self.jacsave = None
         self.nsamples_save = None
 
@@ -118,8 +125,8 @@ def xopt(self, abs_coordinates=False):
     def ropt(self):
         return self.fval_v[self.kopt, :]  # residuals for current iterate
 
-    def fopt(self):
-        return self.fval[self.kopt]
+    def objopt(self):
+        return self.objval[self.kopt]
 
     def xpt(self, k, abs_coordinates=False):
         assert 0 <= k < self.npt(), "Invalid index %g" % k
@@ -135,9 +142,9 @@ def rvec(self, k):
         assert 0 <= k < self.npt(), "Invalid index %g" % k
         return self.fval_v[k, :]
 
-    def fval(self, k):
+    def objval(self, k):
         assert 0 <= k < self.npt(), "Invalid index %g" % k
-        return self.fval[k]
+        return self.objval[k]
 
     def as_absolute_coordinates(self, x, full_dykstra=False):
         # If x were an interpolation point, get the absolute coordinates of x
@@ -177,18 +184,20 @@ def change_point(self, k, x, rvec, allow_kopt_update=True):
 
         self.points[k, :] = x.copy()
         self.fval_v[k, :] = rvec.copy()
-        self.fval[k] = sumsq(rvec)
+        self.objval[k] = sumsq(rvec)
+        if self.h is not None:
+            self.objval[k] += self.h(remove_scaling(self.xbase + x, self.scaling_changes), *self.argsh)
         self.nsamples[k] = 1
         self.factorisation_current = False
 
-        if allow_kopt_update and self.fval[k] < self.fopt():
+        if allow_kopt_update and self.objval[k] < self.objopt():
             self.kopt = k
         return
 
     def swap_points(self, k1, k2):
         self.points[[k1, k2], :] = self.points[[k2, k1], :]
         self.fval_v[[k1, k2], :] = self.fval_v[[k2, k1], :]
-        self.fval[[k1, k2]] = self.fval[[k2, k1]]
+        self.objval[[k1, k2]] = self.objval[[k2, k1]]
         if self.kopt == k1:
             self.kopt = k2
         elif self.kopt == k2:
@@ -201,22 +210,27 @@ def add_new_sample(self, k, rvec_extra):
         assert 0 <= k < self.npt(), "Invalid index %g" % k
         t = float(self.nsamples[k]) / float(self.nsamples[k] + 1)
         self.fval_v[k, :] = t * self.fval_v[k, :] + (1 - t) * rvec_extra
-        self.fval[k] = sumsq(self.fval_v[k, :])
+        # NOTE: how to sample when we have h? still at xpt(k), then add h(xpt(k)). Modify test if incorrect!
+        self.objval[k] = sumsq(self.fval_v[k, :])
+        if self.h is not None:
+            self.objval[k] += self.h(remove_scaling(self.xbase + self.points[k, :], self.scaling_changes), *self.argsh)
         self.nsamples[k] += 1
 
-        self.kopt = np.argmin(self.fval[:self.npt()])  # make sure kopt is always the best value we have
+        self.kopt = np.argmin(self.objval[:self.npt()])  # make sure kopt is always the best value we have
         return
 
     def add_new_point(self, x, rvec):
         self.points = np.append(self.points, x.reshape((1, self.n())), axis=0)  # append row to xpt
         self.fval_v = np.append(self.fval_v, rvec.reshape((1, self.m())), axis=0)  # append row to fval_v
-        f = np.dot(rvec, rvec)
-        self.fval = np.append(self.fval, f)  # append entry to fval
+        obj = sumsq(rvec)
+        if self.h is not None:
+            obj += self.h(remove_scaling(self.xbase + x, self.scaling_changes), *self.argsh)
+        self.objval = np.append(self.objval, obj)  # append entry to fval
         self.nsamples = np.append(self.nsamples, 1)  # add new sample number
         self.num_pts += 1  # make sure npt is updated
         self.npt_so_far += 1
 
-        if f < self.fopt():
+        if obj < self.objopt():
             self.kopt = self.npt() - 1
 
         self.factorisation_current = False
@@ -236,27 +250,30 @@ def shift_base(self, xbase_shift):
         return
 
     def save_point(self, x, rvec, nsamples, x_in_abs_coords=True):
-        f = sumsq(rvec)
-        if self.fsave is None or f <= self.fsave:
-            self.xsave = x.copy() if x_in_abs_coords else self.as_absolute_coordinates(x)
+        xabs = x.copy() if x_in_abs_coords else self.as_absolute_coordinates(x)
+        obj = sumsq(rvec)
+        if self.h is not None:
+            obj += self.h(remove_scaling(xabs, self.scaling_changes), *self.argsh)
+        if self.objsave is None or obj <= self.objsave:
+            self.xsave = xabs
             self.rsave = rvec.copy()
-            self.fsave = f
+            self.objsave = obj
             self.jacsave = self.model_jac.copy()
             self.nsamples_save = nsamples
             return True
         else:
             return False  # this value is worse than what we have already - didn't save
 
     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
-            return self.xopt(abs_coordinates=True).copy(), self.ropt().copy(), self.fopt(), self.model_jac.copy(), self.nsamples[self.kopt]
+        # Return x and objval for optimal point (either from xsave+objsave or kopt)
+        if self.objsave is None or self.objopt() <= self.objsave:  # optimal has changed since xsave+objsave were last set
+            return self.xopt(abs_coordinates=True).copy(), self.ropt().copy(), self.objopt(), self.model_jac.copy(), self.nsamples[self.kopt]
         else:
-            return self.xsave.copy(), self.rsave.copy(), self.fsave, self.jacsave, self.nsamples_save
+            return self.xsave.copy(), self.rsave.copy(), self.objsave, self.jacsave, self.nsamples_save
 
     def min_objective_value(self):
         # Get termination criterion for f small: f <= abs_tol or f <= rel_tol * f0
-        return max(self.abs_tol, self.rel_tol * self.fbeg)
+        return max(self.abs_tol, self.rel_tol * self.objbeg)
 
     def model_value(self, d, d_based_at_xopt=True, with_const_term=False):
         if d_based_at_xopt:
@@ -375,7 +392,7 @@ def interpolate_mini_models_svd(self, verbose=False, make_full_rank=False, min_s
         return True, interp_error, sqrt(norm_J_error), linalg_resid, ls_interp_cond_num  # flag ok
 
     def build_full_model(self):
-        # Build full least squares objective model from mini-models
+        # Build full least squares model from mini-models
         # Centred around xopt
         r = self.model_const + np.dot(self.model_jac, self.xopt())  # constant term (for inexact interpolation)
         J = self.model_jac

dfols/model.py

@@ -82,7 +82,7 @@ def __init__(self, n, npt, maxfun, objfun_has_noise=False):
         self.params["restarts.use_soft_restarts"] = True
         self.params["restarts.soft.num_geom_steps"] = 3
         self.params["restarts.soft.move_xk"] = True
-        self.params["restarts.soft.max_fake_successful_steps"] = maxfun  # number ratio>0 steps below fsave allowed
+        self.params["restarts.soft.max_fake_successful_steps"] = maxfun  # number ratio>0 steps below objsave allowed
         self.params["restarts.hard.use_old_rk"] = True  # recycle r(xk) from previous run?
         self.params["restarts.increase_npt"] = False
         self.params["restarts.increase_npt_amt"] = 1
@@ -109,12 +109,20 @@ def __init__(self, n, npt, maxfun, objfun_has_noise=False):
         self.params["growing.full_rank.min_sing_val"] = 1e-6  # absolute floor on singular values
         self.params["growing.full_rank.svd_max_jac_cond"] = 1e8  # maximum condition number of Jacobian
         self.params["growing.perturb_trust_region_step"] = False  # add random direction onto TRS solution?
+        
         # Dykstra's algorithm
         self.params["dykstra.d_tol"] = 1e-10
         self.params["dykstra.max_iters"] = 100
+        
         # Matrix rank algorithm
         self.params["matrix_rank.r_tol"] = 1e-18
-
+        
+        # Function tolerance when applying S-FISTA method
+        self.params["func_tol.criticality_measure"] = 1e-3
+        self.params["func_tol.tr_step"] = 1-1e-1
+        self.params["func_tol.max_iters"] = 500
+        self.params["sfista.max_iters_scaling"] = 2.0
+        
         self.params_changed = {}
         for p in self.params:
             self.params_changed[p] = False
@@ -268,6 +276,14 @@ def param_type(self, key, npt):
             type_str, nonetype_ok, lower, upper = 'int', False, 0, None
         elif key == "matrix_rank.r_tol":
             type_str, nonetype_ok, lower, upper = 'float', False, 0.0, None
+        elif key == "func_tol.criticality_measure":
+            type_str, nonetype_ok, lower, upper = 'float', False, 0.0, 1.0
+        elif key == "func_tol.tr_step":
+            type_str, nonetype_ok, lower, upper = 'float', False, 0.0, 1.0
+        elif key == "func_tol.max_iters":
+            type_str, nonetype_ok, lower, upper = 'int', False, 0, None
+        elif key == "sfista.max_iters_scaling":
+            type_str, nonetype_ok, lower, upper = 'float', False, 1.0, None
         else:
             assert False, "ParameterList.param_type() has unknown key: %s" % key
         return type_str, nonetype_ok, lower, upper