Repository: Quantum-TII/qibo (python)

Issue:
Result from optimizers
Hi!
I was working with optimizers implemented in Qibo and I think it would be better if the optimizers returns more than only the value of the function and optimal parameters for example why the minimisation process has stoped, number of function evaluation etc...
Methods like "Powell" have also the possibility in SciPy to set bounds (not in dictionary "options") but this is not possible in Qibo.

Let me know what do you think about it. Thank you all!

Pull request:
Standarizing optimizers input and output
Closes #357. This PR propagates the extra flags from scipy to our optimizer and vqe/qaoa minimizer methods, and returns an extra object which holds the original result output from the underlying optimizer library.

Patch:
diff --git a/examples/aavqe/main.py b/examples/aavqe/main.py
index 89062f372..a4d7935cf 100644
--- a/examples/aavqe/main.py
+++ b/examples/aavqe/main.py
@@ -36,9 +36,9 @@ def AAVQE(nqubits, layers, maxsteps, T_max, initial_parameters, easy_hamiltonian
         print('s =',s)
         hamiltonian =  (1-s)*easy_hamiltonian + s*problem_hamiltonian
         vqe = models.VQE(circuit, hamiltonian)
-        energy, params = vqe.minimize(initial_parameters, method='Nelder-Mead',
-                                      options={'maxfev': maxsteps},
-                                      compile=False)
+        energy, params, _ = vqe.minimize(initial_parameters, method='Nelder-Mead',
+                                         options={'maxfev': maxsteps},
+                                         compile=False)
         initial_parameters = params
     return energy, params
 
diff --git a/examples/adiabatic/optimize.py b/examples/adiabatic/optimize.py
index 9b9bded29..08e046342 100644
--- a/examples/adiabatic/optimize.py
+++ b/examples/adiabatic/optimize.py
@@ -53,9 +53,9 @@ def main(nqubits, hfield, params, dt, solver, method, maxiter, save):
     evolution = models.AdiabaticEvolution(h0, h1, spolynomial, dt=dt,
                                           solver=solver)
     options = {"maxiter": maxiter, "disp": True}
-    energy, parameters = evolution.minimize(params, method=method,
-                                            options=options,
-                                            messages=True)
+    energy, parameters, _ = evolution.minimize(params, method=method,
+                                               options=options,
+                                               messages=True)
 
     print("\nBest energy found:", energy)
     print("Final parameters:", parameters)
diff --git a/examples/adiabatic3sat/main.py b/examples/adiabatic3sat/main.py
index 07e7a9384..5800339b4 100644
--- a/examples/adiabatic3sat/main.py
+++ b/examples/adiabatic3sat/main.py
@@ -87,8 +87,8 @@ def main(nqubits, instance, T, dt, solver, plot, trotter, params,
             options = {"nepochs": maxiter}
         else:
             options = {"maxiter": maxiter, "disp": True}
-        energy, params = evolve.minimize(params, method=method,
-                                         options=options)
+        energy, params, _ = evolve.minimize(params, method=method,
+                                            options=options)
         T = params[-1]
 
     # Perform evolution
@@ -126,4 +126,4 @@ if __name__ == "__main__":
     args = vars(parser.parse_args())
     if args["params"] is not None:
         args["params"] = [float(x) for x in args["params"].split(",")]
-    main(**args)
\ No newline at end of file
+    main(**args)
diff --git a/examples/benchmarks/qaoa.py b/examples/benchmarks/qaoa.py
index 1efde3d3a..145fc8f10 100644
--- a/examples/benchmarks/qaoa.py
+++ b/examples/benchmarks/qaoa.py
@@ -36,8 +36,8 @@ def main(nqubits, nangles, trotter=False, solver="exp",
 
     start_time = time.time()
     options = {'disp': True, 'maxiter': maxiter}
-    best, params = qaoa.minimize(initial_parameters, method=method,
-                                 options=options)
+    best, params, _ = qaoa.minimize(initial_parameters, method=method,
+                                    options=options)
     minimization_time = time.time() - start_time
 
     epsilon = np.log10(1/np.abs(best - target))
diff --git a/examples/benchmarks/vqe.py b/examples/benchmarks/vqe.py
index 1a13a4b39..1eb5d1bbb 100644
--- a/examples/benchmarks/vqe.py
+++ b/examples/benchmarks/vqe.py
@@ -67,8 +67,8 @@ def main(nqubits, nlayers, varlayer=False, method="Powell", maxiter=None):
 
     start_time = time.time()
     options = {'disp': True, 'maxiter': maxiter}
-    best, params = vqe.minimize(initial_parameters, method=method,
-                                options=options, compile=False)
+    best, params, _ = vqe.minimize(initial_parameters, method=method,
+                                   options=options, compile=False)
     minimization_time = time.time() - start_time
     epsilon = np.log10(1/np.abs(best-target))
     print("Found state =", best)
diff --git a/examples/qPDF/qPDF.ipynb b/examples/qPDF/qPDF.ipynb
index dc4823642..80862eea8 100644
--- a/examples/qPDF/qPDF.ipynb
+++ b/examples/qPDF/qPDF.ipynb
@@ -207,7 +207,7 @@
     "\n",
     "np.random.seed(10)\n",
     "params = np.random.rand(mypdf.nparams)\n",
-    "_, params = optimize(loss, params, method='cma')\n",
+    "_, params, _ = optimize(loss, params, method='cma')\n",
     "```"
    ]
   },
diff --git a/src/qibo/evolution.py b/src/qibo/evolution.py
index 781aa744c..d3d0302e2 100644
--- a/src/qibo/evolution.py
+++ b/src/qibo/evolution.py
@@ -326,11 +326,11 @@ class AdiabaticEvolution(StateEvolution):
         else:
             loss = lambda p, ae, h1, msg, hist: self._loss(p, ae, h1, msg, hist).numpy()
 
-        result, parameters = optimizers.optimize(loss, initial_parameters,
+        result, parameters, extra = optimizers.optimize(loss, initial_parameters,
                                                  args=(self, self.h1, self.opt_messages, self.opt_history),
                                                  method=method, options=options)
         if isinstance(parameters, K.tensor_types) and not len(parameters.shape): # pragma: no cover
             # some optimizers like ``Powell`` return number instead of list
             parameters = [parameters]
         self.set_parameters(parameters)
-        return result, parameters
+        return result, parameters, extra
diff --git a/src/qibo/models.py b/src/qibo/models.py
index 4a9586679..f46bba1ca 100644
--- a/src/qibo/models.py
+++ b/src/qibo/models.py
@@ -172,7 +172,9 @@ class VQE(object):
         self.circuit = circuit
         self.hamiltonian = hamiltonian
 
-    def minimize(self, initial_state, method='Powell', options=None, compile=False, processes=None):
+    def minimize(self, initial_state, method='Powell', jac=None, hess=None,
+                 hessp=None, bounds=None, constraints=(), tol=None, callback=None,
+                 options=None, compile=False, processes=None):
         """Search for parameters which minimizes the hamiltonian expectation.
 
         Args:
@@ -181,6 +183,14 @@ class VQE(object):
             method (str): the desired minimization method.
                 See :meth:`qibo.optimizers.optimize` for available optimization
                 methods.
+            jac (dict): Method for computing the gradient vector for scipy optimizers.
+            hess (dict): Method for computing the hessian matrix for scipy optimizers.
+            hessp (callable): Hessian of objective function times an arbitrary
+                vector for scipy optimizers.
+            bounds (sequence or Bounds): Bounds on variables for scipy optimizers.
+            constraints (dict): Constraints definition for scipy optimizers.
+            tol (float): Tolerance of termination for scipy optimizers.
+            callback (callable): Called after each iteration for scipy optimizers.
             options (dict): a dictionary with options for the different optimizers.
             compile (bool): whether the TensorFlow graph should be compiled.
             processes (int): number of processes when using the paralle BFGS method.
@@ -188,6 +198,10 @@ class VQE(object):
         Return:
             The final expectation value.
             The corresponding best parameters.
+            The optimization result object. For scipy methods it
+                returns the ``OptimizeResult``, for ``'cma'`` the
+                ``CMAEvolutionStrategy.result``, and for ``'sgd'``
+                the options used during the optimization.
         """
         def _loss(params, circuit, hamiltonian):
             circuit.set_parameters(params)
@@ -212,12 +226,14 @@ class VQE(object):
             loss = _loss
         else:
             loss = lambda p, c, h: _loss(p, c, h).numpy()
-        result, parameters = self.optimizers.optimize(loss, initial_state,
-                                                      args=(self.circuit, self.hamiltonian),
-                                                      method=method, options=options,
-                                                      compile=compile, processes=processes)
+        result, parameters, extra = self.optimizers.optimize(loss, initial_state,
+                                                             args=(self.circuit, self.hamiltonian),
+                                                             method=method, jac=jac, hess=hess, hessp=hessp,
+                                                             bounds=bounds, constraints=constraints,
+                                                             tol=tol, callback=callback, options=options,
+                                                             compile=compile, processes=processes)
         self.circuit.set_parameters(parameters)
-        return result, parameters
+        return result, parameters, extra
 
 
 class QAOA(object):
@@ -359,7 +375,9 @@ class QAOA(object):
             return self.state_cls.plus_state(self.nqubits).tensor
         return StateCircuit.get_initial_state(self, state)
 
-    def minimize(self, initial_p, initial_state=None, method='Powell', options=None):
+    def minimize(self, initial_p, initial_state=None, method='Powell',
+                 jac=None, hess=None, hessp=None, bounds=None, constraints=(),
+                 tol=None, callback=None, options=None, compile=False, processes=None):
         """Optimizes the variational parameters of the QAOA.
 
         Args:
@@ -368,11 +386,25 @@ class QAOA(object):
             method (str): the desired minimization method.
                 See :meth:`qibo.optimizers.optimize` for available optimization
                 methods.
+            jac (dict): Method for computing the gradient vector for scipy optimizers.
+            hess (dict): Method for computing the hessian matrix for scipy optimizers.
+            hessp (callable): Hessian of objective function times an arbitrary
+                vector for scipy optimizers.
+            bounds (sequence or Bounds): Bounds on variables for scipy optimizers.
+            constraints (dict): Constraints definition for scipy optimizers.
+            tol (float): Tolerance of termination for scipy optimizers.
+            callback (callable): Called after each iteration for scipy optimizers.
             options (dict): a dictionary with options for the different optimizers.
+            compile (bool): whether the TensorFlow graph should be compiled.
+            processes (int): number of processes when using the paralle BFGS method.
 
         Return:
             The final energy (expectation value of the ``hamiltonian``).
             The corresponding best parameters.
+            The optimization result object. For scipy methods it
+                returns the ``OptimizeResult``, for ``'cma'`` the
+                ``CMAEvolutionStrategy.result``, and for ``'sgd'``
+                the options used during the optimization.
         """
         if len(initial_p) % 2 != 0:
             raise_error(ValueError, "Initial guess for the parameters must "
@@ -390,7 +422,10 @@ class QAOA(object):
         else:
             loss = lambda p, c, h: _loss(p, c, h).numpy()
 
-        result, parameters = self.optimizers.optimize(loss, initial_p, args=(self, self.hamiltonian),
-                                                      method=method, options=options)
+        result, parameters, extra = self.optimizers.optimize(loss, initial_p, args=(self, self.hamiltonian),
+                                                             method=method, jac=jac, hess=hess, hessp=hessp,
+                                                             bounds=bounds, constraints=constraints,
+                                                             tol=tol, callback=callback, options=options,
+                                                             compile=compile, processes=processes)
         self.set_parameters(parameters)
-        return result, parameters
+        return result, parameters, extra
diff --git a/src/qibo/optimizers.py b/src/qibo/optimizers.py
index a412de643..8f8145ecb 100644
--- a/src/qibo/optimizers.py
+++ b/src/qibo/optimizers.py
@@ -1,5 +1,6 @@
 def optimize(loss, initial_parameters, args=(), method='Powell',
-             options=None, compile=False, processes=None):
+             jac=None, hess=None, hessp=None, bounds=None, constraints=(),
+             tol=None, callback=None, options=None, compile=False, processes=None):
     """Main optimization method. Selects one of the following optimizers:
         - :meth:`qibo.optimizers.cma`
         - :meth:`qibo.optimizers.newtonian`
@@ -15,12 +16,27 @@ def optimize(loss, initial_parameters, args=(), method='Powell',
         method (str): Name of optimizer to use. Can be ``'cma'``, ``'sgd'`` or
             one of the Newtonian methods supported by
             :meth:`qibo.optimizers.newtonian` and ``'parallel_L-BFGS-B'``.
+        jac (dict): Method for computing the gradient vector for scipy optimizers.
+        hess (dict): Method for computing the hessian matrix for scipy optimizers.
+        hessp (callable): Hessian of objective function times an arbitrary
+            vector for scipy optimizers.
+        bounds (sequence or Bounds): Bounds on variables for scipy optimizers.
+        constraints (dict): Constraints definition for scipy optimizers.
+        tol (float): Tolerance of termination for scipy optimizers.
+        callback (callable): Called after each iteration for scipy optimizers.
         options (dict): Dictionary with options. See the specific optimizer
             bellow for a list of the supported options.
         compile (bool): If ``True`` the Tensorflow optimization graph is compiled.
             This is relevant only for the ``'sgd'`` optimizer.
         processes (int): number of processes when using the parallel BFGS method.
 
+    Returns:
+        loss (float): final best loss value.
+        xbest (float): best parameters obtained by the optimizer.
+        extra: optimizer specific return object containing. For scipy methods it
+            returns the ``OptimizeResult``, for ``'cma'`` the ``CMAEvolutionStrategy.result``,
+            and for ``'sgd'`` the options used during the optimization.
+
     Example:
         ::
 
@@ -40,7 +56,7 @@ def optimize(loss, initial_parameters, args=(), method='Powell',
 
             # optimize using random initial variational parameters
             initial_parameters = np.random.uniform(0, 2, 1)
-            best, params = optimize(myloss, initial_parameters, args=(circuit))
+            best, params, extra = optimize(myloss, initial_parameters, args=(circuit))
 
             # set parameters to circuit
             circuit.set_parameters(params)
@@ -50,7 +66,9 @@ def optimize(loss, initial_parameters, args=(), method='Powell',
     elif method == "sgd":
         return sgd(loss, initial_parameters, args, options, compile)
     else:
-        return newtonian(loss, initial_parameters, args, method, options, processes)
+        return newtonian(loss, initial_parameters, args, method,
+                         jac, hess, hessp, bounds, constraints, tol,
+                         callback, options, processes)
 
 
 def cma(loss, initial_parameters, args=(), options=None):
@@ -68,10 +86,12 @@ def cma(loss, initial_parameters, args=(), options=None):
     """
     import cma
     r = cma.fmin2(loss, initial_parameters, 1.7, options=options, args=args)
-    return r[1].result.fbest, r[1].result.xbest
+    return r[1].result.fbest, r[1].result.xbest, r
 
 
-def newtonian(loss, initial_parameters, args=(), method='Powell', options=None, processes=None):
+def newtonian(loss, initial_parameters, args=(), method='Powell',
+              jac=None, hess=None, hessp=None, bounds=None, constraints=(),
+              tol=None, callback=None, options=None, processes=None):
     """Newtonian optimization approaches based on ``scipy.optimize.minimize``.
 
     For more details check the `scipy documentation <https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html>`_.
@@ -94,6 +114,14 @@ def newtonian(loss, initial_parameters, args=(), method='Powell', options=None,
         args (tuple): optional arguments for the loss function.
         method (str): Name of method supported by ``scipy.optimize.minimize`` and ``'parallel_L-BFGS-B'`` for
             a parallel version of L-BFGS-B algorithm.
+        jac (dict): Method for computing the gradient vector for scipy optimizers.
+        hess (dict): Method for computing the hessian matrix for scipy optimizers.
+        hessp (callable): Hessian of objective function times an arbitrary
+            vector for scipy optimizers.
+        bounds (sequence or Bounds): Bounds on variables for scipy optimizers.
+        constraints (dict): Constraints definition for scipy optimizers.
+        tol (float): Tolerance of termination for scipy optimizers.
+        callback (callable): Called after each iteration for scipy optimizers.
         options (dict): Dictionary with options accepted by
             ``scipy.optimize.minimize``.
         processes (int): number of processes when using the parallel BFGS method.
@@ -101,12 +129,15 @@ def newtonian(loss, initial_parameters, args=(), method='Powell', options=None,
     if method == 'parallel_L-BFGS-B':
         from qibo.parallel import _check_parallel_configuration
         _check_parallel_configuration(processes)
-        o = ParallelBFGS(loss, args=args, options=options, processes=processes)
+        o = ParallelBFGS(loss, args=args, processes=processes,
+                         bounds=bounds, callback=callback, options=options)
         m = o.run(initial_parameters)
     else:
         from scipy.optimize import minimize
-        m = minimize(loss, initial_parameters, args=args, method=method, options=options)
-    return m.fun, m.x
+        m = minimize(loss, initial_parameters, args=args, method=method,
+                     jac=jac, hess=hess, hessp=hessp, bounds=bounds, constraints=constraints,
+                     tol=tol, callback=callback, options=options)
+    return m.fun, m.x, m
 
 
 def sgd(loss, initial_parameters, args=(), options=None, compile=False):
@@ -171,7 +202,7 @@ def sgd(loss, initial_parameters, args=(), options=None, compile=False):
         if e % sgd_options["nmessage"] == 1:
             log.info('ite %d : loss %f', e, l.numpy())
 
-    return loss(vparams, *args).numpy(), vparams.numpy()
+    return loss(vparams, *args).numpy(), vparams.numpy(), sgd_options
 
 
 from qibo.parallel import ParallelResources, _executor


Test patch:
diff --git a/src/qibo/tests/test_evolution.py b/src/qibo/tests/test_evolution.py
index 36c48132c..6d3a7ea5a 100644
--- a/src/qibo/tests/test_evolution.py
+++ b/src/qibo/tests/test_evolution.py
@@ -323,7 +323,7 @@ def test_scheduling_optimization(method, options, messages, trotter, filename):
     h1 = hamiltonians.TFIM(3, trotter=trotter)
     sp = lambda t, p: (1 - p) * np.sqrt(t) + p * t
     adevp = models.AdiabaticEvolution(h0, h1, sp, dt=1e-1)
-    best, params = adevp.minimize([0.5, 1], method=method, options=options,
-                                  messages=messages)
+    best, params, _ = adevp.minimize([0.5, 1], method=method, options=options,
+                                     messages=messages)
     if filename is not None:
         utils.assert_regression_fixture(params, filename)
diff --git a/src/qibo/tests/test_variational.py b/src/qibo/tests/test_variational.py
index 6126c7d22..d71e4d061 100644
--- a/src/qibo/tests/test_variational.py
+++ b/src/qibo/tests/test_variational.py
@@ -46,7 +46,7 @@ def test_vqc(method, options, compile, filename):
     data = np.random.normal(0, 1, size=2**nqubits)
 
     # perform optimization
-    best, params = optimize(myloss, x0, args=(c, data), method=method,
+    best, params, _ = optimize(myloss, x0, args=(c, data), method=method,
                             options=options, compile=compile)
     if filename is not None:
         utils.assert_regression_fixture(params, filename)
@@ -100,8 +100,8 @@ def test_vqe(method, options, compile, filename):
     np.random.seed(0)
     initial_parameters = np.random.uniform(0, 2*np.pi, 2*nqubits*layers + nqubits)
     v = VQE(circuit, hamiltonian)
-    best, params = v.minimize(initial_parameters, method=method,
-                              options=options, compile=compile)
+    best, params, _ = v.minimize(initial_parameters, method=method,
+                                 options=options, compile=compile)
     if method == "cma":
         # remove `outcmaes` folder
         import shutil
@@ -134,12 +134,12 @@ def test_vqe_custom_gates_errors():
     v = VQE(circuit, hamiltonian)
     # compile with custom gates
     with pytest.raises(RuntimeError):
-        best, params = v.minimize(initial_parameters, method="BFGS",
-                                  options={'maxiter': 1}, compile=True)
+        best, params, _ = v.minimize(initial_parameters, method="BFGS",
+                                     options={'maxiter': 1}, compile=True)
     # use SGD with custom gates
     with pytest.raises(RuntimeError):
-        best, params = v.minimize(initial_parameters, method="sgd",
-                                  compile=False)
+        best, params, _ = v.minimize(initial_parameters, method="sgd",
+                                     compile=False)
     qibo.set_backend(original_backend)
 
 
@@ -252,6 +252,6 @@ def test_qaoa_optimization(method, options, trotter, filename):
     h = hamiltonians.XXZ(3, trotter=trotter)
     qaoa = models.QAOA(h)
     initial_p = [0.05, 0.06, 0.07, 0.08]
-    best, params = qaoa.minimize(initial_p, method=method, options=options)
+    best, params, _ = qaoa.minimize(initial_p, method=method, options=options)
     if filename is not None:
         utils.assert_regression_fixture(params, filename)