Repository: CUQI-DTU/CUQIpy

Issue:
Rename existing solvers
**General description of task**
- Class names of [solvers](https://github.com/CUQI-DTU/CUQIpy/blob/main/cuqi/solver/_solver.py) are not following a consistent convection, e.g., `minimize` and `LS`
- Unit tests for several solvers are missing

**Definition of Done (Feature/change)**
- [x] Rename relevant solvers
  - L_BFGS_B -> ScipyLBFGSB
  - minimize -> ScipyMinimizer
  - maximize -> ScipyMaximizer
  - LS -> ScipyLeastSquares
- [x] Any calls on relevant solvers are updated
- [x] Add unit tests for solvers

**Definition of Done (Mandatory)**
- [x] Documentation added (docstrings on all public methods/classes)
- [x] Unit tests added/updated (and passing!)
- [x] Code reviewed and approved by 2 CUQI-maintainers
- [x] Online documentation is rendered correctly (esp. math) 

**Defintion of Done (Optional)**

Discussion:
@chaozg: Hi @jakobsj, could you review the DoD here?@jakobsj: Thanks @chaozg I think everything makes sense, would like one addition namely to go through existing code base and make sure to update any calls to solvers to use the new names. May be covered by the "Rename" task but better to be explicit so that it is clear it has been done. If you agree, then with that DoD is approved.@chaozg: Thank @jakobsj . All calls to the relevant solvers have been updated and I just made it clear in the DoD.

Solution patch:
diff --git a/cuqi/problem/_problem.py b/cuqi/problem/_problem.py
index 824a2d878..455fe7704 100644
--- a/cuqi/problem/_problem.py
+++ b/cuqi/problem/_problem.py
@@ -771,11 +771,11 @@ def gradfunc(x): return -density.gradient(x)
         if self._check_posterior(self, CMRF, must_have_gradient=True): # Use L-BFGS-B for CMRF prior as it has better performance for this multi-modal posterior
             if disp: print(f"Using scipy.optimize.L_BFGS_B on negative log of {density.__class__.__name__}")
             if disp: print("x0: ones vector")
-            solver = cuqi.solver.L_BFGS_B(func, x0, gradfunc=gradfunc)
+            solver = cuqi.solver.ScipyLBFGSB(func, x0, gradfunc=gradfunc)
         else:
             if disp: print(f"Using scipy.optimize.minimize on negative log of {density.__class__.__name__}")
             if disp: print("x0: ones vector")
-            solver = cuqi.solver.minimize(func, x0, gradfunc=gradfunc)
+            solver = cuqi.solver.ScipyMinimizer(func, x0, gradfunc=gradfunc)
 
         x_MAP, solver_info = solver.solve()
 
diff --git a/cuqi/solver/__init__.py b/cuqi/solver/__init__.py
index 5b9e11a19..c3851ae65 100644
--- a/cuqi/solver/__init__.py
+++ b/cuqi/solver/__init__.py
@@ -1,8 +1,8 @@
 from ._solver import (
-    L_BFGS_B,
-    minimize,
-    maximize,
-    LS,
+    ScipyLBFGSB,
+    ScipyMinimizer,
+    ScipyMaximizer,
+    ScipyLeastSquares,
     CGLS,
     LM,
     PDHG,
diff --git a/cuqi/solver/_solver.py b/cuqi/solver/_solver.py
index 98a52a389..5c7868c2d 100644
--- a/cuqi/solver/_solver.py
+++ b/cuqi/solver/_solver.py
@@ -15,7 +15,7 @@
     has_cholmod = False
 
 
-class L_BFGS_B(object):
+class ScipyLBFGSB(object):
     """Wrapper for :meth:`scipy.optimize.fmin_l_bfgs_b`.
 
     Minimize a function func using the L-BFGS-B algorithm.
@@ -30,14 +30,10 @@ class L_BFGS_B(object):
         Initial guess.
     gradfunc : callable f(x,*args), optional
         The gradient of func. 
-        If None, then the solver approximates the gradient.
+        If None, the solver approximates the gradient with a finite difference scheme.
     kwargs : keyword arguments passed to scipy's L-BFGS-B algorithm. See documentation for scipy.optimize.minimize
-
-    Methods
-    ----------
-    :meth:`solve`: Runs the solver and returns the solution and info about the optimization.
     """
-    def __init__(self,func,x0, gradfunc = None, **kwargs):
+    def __init__(self, func, x0, gradfunc = None, **kwargs):
         self.func= func
         self.x0 = x0
         self.gradfunc = gradfunc
@@ -83,7 +79,7 @@ def solve(self):
                 "nfev": solution[2]['funcalls']}
         return solution[0], info
 
-class minimize(object):
+class ScipyMinimizer(object):
     """Wrapper for :meth:`scipy.optimize.minimize`.
 
     Minimize a function func using scipy's optimize.minimize module.
@@ -115,12 +111,8 @@ class minimize(object):
         ‘trust-krylov’ 
         If not given, chosen to be one of BFGS, L-BFGS-B, SLSQP, depending if the problem has constraints or bounds.
     kwargs : keyword arguments passed to scipy's minimizer. See documentation for scipy.optimize.minimize
-
-    Methods
-    ----------
-    :meth:`solve`: Runs the solver and returns the solution and info about the optimization.
     """
-    def __init__(self,func,x0, gradfunc = None, method = None, **kwargs):
+    def __init__(self, func, x0, gradfunc = '2-point', method = None, **kwargs):
         self.func= func
         self.x0 = x0
         self.method = method
@@ -147,18 +139,20 @@ def solve(self):
         info = {"success": solution['success'],
                 "message": solution['message'],
                 "func": solution['fun'],
-                "grad": solution['jac'],
                 "nit": solution['nit'], 
                 "nfev": solution['nfev']}
+        # if gradfunc is callable, record the gradient in the info dict
+        if 'jac' in solution.keys():
+            info['grad'] = solution['jac']
         if isinstance(self.x0,CUQIarray):
             sol = CUQIarray(solution['x'],geometry=self.x0.geometry)
         else:
             sol = solution['x']
         return sol, info
 
-class maximize(minimize):
-    """Simply calls ::class:: cuqi.solver.minimize with -func."""
-    def __init__(self,func,x0, gradfunc = None, method = None, **kwargs):
+class ScipyMaximizer(ScipyMinimizer):
+    """Simply calls ::class:: cuqi.solver.ScipyMinimizer with -func."""
+    def __init__(self, func, x0, gradfunc = None, method = None, **kwargs):
         def nfunc(*args,**kwargs):
             return -func(*args,**kwargs)
         if gradfunc is not None:
@@ -170,7 +164,7 @@ def ngradfunc(*args,**kwargs):
 
 
 
-class LS(object):
+class ScipyLeastSquares(object):
     """Wrapper for :meth:`scipy.optimize.least_squares`.
 
     Solve nonlinear least-squares problems with bounds:
@@ -189,7 +183,7 @@ class LS(object):
         Initial guess.
     Jac : callable f(x,*args), optional
         The Jacobian of func. 
-        If None, then the solver approximates the Jacobian.
+        If not specified, the solver approximates the Jacobian with a finite difference scheme.
     loss: callable rho(x,*args)
         Determines the loss function
         'linear' : rho(z) = z. Gives a standard least-squares problem.
@@ -203,7 +197,7 @@ class LS(object):
         'dogbox', dogleg algorithm with rectangular trust regions, for small problems with bounds.
         'lm', Levenberg-Marquardt algorithm as implemented in MINPACK. Doesn't handle bounds and sparse Jacobians.
     """
-    def __init__(self, func, x0, jacfun=None, method='trf', loss='linear', tol=1e-6, maxit=1e4):
+    def __init__(self, func, x0, jacfun='2-point', method='trf', loss='linear', tol=1e-6, maxit=1e4):
         self.func = func
         self.x0 = x0
         self.jacfun = jacfun
diff --git a/demos/demo12_gradients.py b/demos/demo12_gradients.py
index afb9e74c4..9940df3ff 100644
--- a/demos/demo12_gradients.py
+++ b/demos/demo12_gradients.py
@@ -33,7 +33,7 @@
     # Solve posterior problem using BFGS
     def f(x): return -TP.posterior.logd(x)
     def gradf(x): return -TP.posterior.gradient(x)
-    solver = cuqi.solver.L_BFGS_B(f, x0, gradf)
+    solver = cuqi.solver.ScipyLBFGSB(f, x0, gradf)
     x_MAP, solution_info = solver.solve()
 print('relative error MAP:', np.linalg.norm(x_MAP-x_true)/np.linalg.norm(x_true))
 
diff --git a/demos/try09.py b/demos/try09.py
index 7051e1ca3..1b6d1c6e5 100644
--- a/demos/try09.py
+++ b/demos/try09.py
@@ -48,28 +48,28 @@ def posterior_logpdf_grad(x):
 print('relative error exact MAP:', np.linalg.norm(x_MAP_exact-x_true)/np.linalg.norm(x_true))
 
 #%% L_BFGS_B MAP
-solver = cuqi.solver.L_BFGS_B(posterior_logpdf, x0, gradfunc = posterior_logpdf_grad)
+solver = cuqi.solver.ScipyLBFGSB(posterior_logpdf, x0, gradfunc = posterior_logpdf_grad)
 x_MAP_LBFGS, info_MAP_LBFGS = solver.solve()
 print('relative error L-BFGS MAP:', np.linalg.norm(x_MAP_LBFGS-x_true)/np.linalg.norm(x_true))
 
-#%% L_BFGS_B MAP using minimize
-solver = cuqi.solver.minimize(posterior_logpdf, x0, gradfunc = posterior_logpdf_grad, method = "L-BFGS-B")
+#%% L_BFGS_B MAP using ScipyMinimizer
+solver = cuqi.solver.ScipyMinimizer(posterior_logpdf, x0, gradfunc = posterior_logpdf_grad, method = "L-BFGS-B")
 x_MAP_LBFGS2, info_MAP_LBFGS2 = solver.solve()
 print('relative error L-BFGS MAP:', np.linalg.norm(x_MAP_LBFGS2-x_true)/np.linalg.norm(x_true))
 
 #%% BFGS MAP without gradient input
-solver = cuqi.solver.minimize(posterior_logpdf, x0)
+solver = cuqi.solver.ScipyMinimizer(posterior_logpdf, x0)
 x_MAP_BFGS, info_MAP_BFGS = solver.solve()
 print('relative error BFGS MAP:', np.linalg.norm(x_MAP_BFGS-x_true)/np.linalg.norm(x_true))
 
 #%% BFGS MAP with gradient input
-solver = cuqi.solver.minimize(posterior_logpdf, x0, gradfunc = posterior_logpdf_grad)
+solver = cuqi.solver.ScipyMinimizer(posterior_logpdf, x0, gradfunc = posterior_logpdf_grad)
 x_MAP_BFGSgrad, info_MAP_BFGSgrad = solver.solve()
 print('relative error BFGS MAP:', np.linalg.norm(x_MAP_BFGSgrad-x_true)/np.linalg.norm(x_true))
 
 
 #%% SLSQP MAP
-solver = cuqi.solver.minimize(posterior_logpdf, x0, method = 'SLSQP')
+solver = cuqi.solver.ScipyMinimizer(posterior_logpdf, x0, method = 'SLSQP')
 x_MAP_SLSQP, info_MAP_SLSQP = solver.solve()
 print('relative error SLSQP MAP:', np.linalg.norm(x_MAP_SLSQP-x_true)/np.linalg.norm(x_true))
 
@@ -95,12 +95,12 @@ def likelihood_logpdf_grad(x):
     return grad
 
 # L_BFGS_B MAP
-solver = cuqi.solver.L_BFGS_B(likelihood_logpdf, x0, gradfunc = likelihood_logpdf_grad)
+solver = cuqi.solver.ScipyLBFGSB(likelihood_logpdf, x0, gradfunc = likelihood_logpdf_grad)
 x_ML_LBFGS, info_ML_LBFGS = solver.solve()
 print('relative error L-BFGS ML:', np.linalg.norm(x_MAP_LBFGS-x_true)/np.linalg.norm(x_true))
 
 # BFGS MAP
-solver = cuqi.solver.minimize(likelihood_logpdf, x0)
+solver = cuqi.solver.ScipyMinimizer(likelihood_logpdf, x0)
 x_ML_BFGS, info_ML_BFGS = solver.solve()
 print('relative error BFGS ML:', np.linalg.norm(x_MAP_BFGS-x_true)/np.linalg.norm(x_true))
 


Test patch:
diff --git a/tests/test_solver.py b/tests/test_solver.py
index 8092ee208..bdbb5ba3e 100644
--- a/tests/test_solver.py
+++ b/tests/test_solver.py
@@ -1,7 +1,7 @@
 import numpy as np
 import scipy as sp
 
-from cuqi.solver import CGLS, LM, FISTA, ADMM, ProximalL1, ProjectNonnegative
+from cuqi.solver import ScipyLBFGSB, ScipyMinimizer, ScipyLeastSquares, CGLS, LM, FISTA, ADMM, ProximalL1, ProjectNonnegative
 from scipy.optimize import lsq_linear
 
 
@@ -20,6 +20,61 @@ def test_CGLS():
     # Compare
     assert np.allclose(sol, ref_sol, rtol=1e-3)
 
+def test_ScipyMinimizer_without_gradient():
+    solver = ScipyMinimizer(sp.optimize.rosen,
+                            np.array([1.3, 0.7, 0.8, 1.9, 1.2]),
+                            method='Nelder-Mead',
+                            tol=1e-6)
+    sol, _ = solver.solve()
+    sol_ref = np.array([1.0, 1.0, 1.0, 1.0, 1.0])
+    assert np.allclose(sol, sol_ref)
+
+def test_ScipyMinimizer_with_gradient():
+    solver = ScipyMinimizer(sp.optimize.rosen,
+                            np.array([1.3, 0.7, 0.8, 1.9, 1.2]),
+                            gradfunc=sp.optimize.rosen_der,
+                            method='BFGS',
+                            tol=1e-6)
+    sol, _ = solver.solve()
+    sol_ref = np.array([1.0, 1.0, 1.0, 1.0, 1.0])
+    assert np.allclose(sol, sol_ref)
+
+def test_ScipyLBFGSB_without_gradient():
+    solver = ScipyLBFGSB(sp.optimize.rosen,
+                         np.array([1.3, 0.7, 0.8, 1.9, 1.2]))
+    sol, _ = solver.solve()
+    sol_ref = np.array([1.0, 1.0, 1.0, 1.0, 1.0])
+    assert np.allclose(sol, sol_ref)
+
+def test_ScipyLBFGSB_with_gradient():
+    solver = ScipyLBFGSB(sp.optimize.rosen,
+                         np.array([1.3, 0.7, 0.8, 1.9, 1.2]),
+                         gradfunc=sp.optimize.rosen_der)
+    sol, _ = solver.solve()
+    sol_ref = np.array([1.0, 1.0, 1.0, 1.0, 1.0])
+    assert np.allclose(sol, sol_ref)
+
+def test_ScipyLeastSquares_without_Jac():
+    def fun_rosenbrock(x):
+        return np.array([10 * (x[1] - x[0]**2), (1 - x[0])])
+    x0 = np.array([2, 2])
+    solver = ScipyLeastSquares(fun_rosenbrock, x0)
+    sol, _ = solver.solve()
+    sol_ref = np.array([1, 1])
+    assert np.allclose(sol, sol_ref)
+
+def test_ScipyLeastSquares_with_Jac():
+    def fun_rosenbrock(x):
+        return np.array([10 * (x[1] - x[0]**2), (1 - x[0])])
+    def jac_rosenbrock(x):
+        return np.array([
+            [-20 * x[0], 10],
+            [-1, 0]])
+    x0 = np.array([2, 2])
+    solver = ScipyLeastSquares(fun_rosenbrock, x0, jacfun=jac_rosenbrock)
+    sol, _ = solver.solve()
+    sol_ref = np.array([1, 1])
+    assert np.allclose(sol, sol_ref)
 
 def test_LM():
     # compare to MATLAB's original code solution