fix: fix initialization bug in primal dual solvers + add non persiste…

…nt change in initialization variable

doc/tutorials/simple.rst

@@ -66,9 +66,9 @@ And finally solve the problem :
 
    >>> x0 = [0., 0., 0., 0.]
    >>> ret = solvers.solve([f2, f1], x0, solver, atol=1e-5, verbosity='HIGH')
+   INFO: Forward-backward method
        func evaluation: 0.000000e+00
        norm_l2 evaluation: 1.260000e+02
-   INFO: Forward-backward method
    Iteration 1 of forward_backward:
        func evaluation: 0.000000e+00
        norm_l2 evaluation: 1.400000e+01

pyunlocbox/solvers.py

@@ -54,7 +54,7 @@
 
 
 def solve(functions, x0, solver=None, atol=None, dtol=None, rtol=1e-3,
-          xtol=None, maxit=200, verbosity='LOW'):
+          xtol=None, maxit=200, verbosity='LOW', inplace=False):
     r"""
     Solve an optimization problem whose objective function is the sum of some
     convex functions.
@@ -79,10 +79,7 @@ def solve(functions, x0, solver=None, atol=None, dtol=None, rtol=1e-3,
         documentation of the considered solver.
     x0 : array_like
         Starting point of the algorithm, :math:`x_0 \in \mathbb{R}^{n \times
-        N}`. Note that if you pass a numpy array it will be modified in place
-        during execution to save memory. It will then contain the solution. Be
-        careful to pass data of the type (int, float32, float64) you want your
-        computations to use.
+        N}`.
     solver : solver class instance, optional
         The solver algorithm. It is an object who must inherit from
         :class:`pyunlocbox.solvers.solver` and implement the :meth:`_pre`,
@@ -111,6 +108,11 @@ def solve(functions, x0, solver=None, atol=None, dtol=None, rtol=1e-3,
         convergence, ``'HIGH'`` for info at all solving steps, ``'ALL'`` for
         all possible outputs, including at each steps of the proximal operators
         computation. Default is ``'LOW'``.
+    inplace : bool, optional
+        If True and x0 is a numpy array, then x0 will be modified in place
+        during execution to save memory. It will then contain the solution. Be
+        careful to pass data of the type (int, float32, float64) you want your
+        computations to use.
 
     Returns
     -------
@@ -143,28 +145,28 @@ def solve(functions, x0, solver=None, atol=None, dtol=None, rtol=1e-3,
     >>> ret = solvers.solve([f], x0, atol=1e-2, verbosity='ALL')
     INFO: Dummy objective function added.
     INFO: Selected solver: forward_backward
-        norm_l2 evaluation: 1.260000e+02
-        dummy evaluation: 0.000000e+00
     INFO: Forward-backward method
+        dummy evaluation: 0.000000e+00
+        norm_l2 evaluation: 1.260000e+02
     Iteration 1 of forward_backward:
-        norm_l2 evaluation: 1.400000e+01
         dummy evaluation: 0.000000e+00
+        norm_l2 evaluation: 1.400000e+01
         objective = 1.40e+01
     Iteration 2 of forward_backward:
-        norm_l2 evaluation: 2.963739e-01
         dummy evaluation: 0.000000e+00
+        norm_l2 evaluation: 2.963739e-01
         objective = 2.96e-01
     Iteration 3 of forward_backward:
-        norm_l2 evaluation: 7.902529e-02
         dummy evaluation: 0.000000e+00
+        norm_l2 evaluation: 7.902529e-02
         objective = 7.90e-02
     Iteration 4 of forward_backward:
-        norm_l2 evaluation: 5.752265e-02
         dummy evaluation: 0.000000e+00
+        norm_l2 evaluation: 5.752265e-02
         objective = 5.75e-02
     Iteration 5 of forward_backward:
-        norm_l2 evaluation: 5.142032e-03
         dummy evaluation: 0.000000e+00
+        norm_l2 evaluation: 5.142032e-03
         objective = 5.14e-03
     Solution found after 5 iterations:
         objective function f(sol) = 5.142032e-03
@@ -198,6 +200,9 @@ def solve(functions, x0, solver=None, atol=None, dtol=None, rtol=1e-3,
     [0.05752265..., 0], [0.00514203..., 0]]
 
     """
+    # to prevent any modification of the input
+    if not(inplace):
+        x0 = x0.copy()
 
     if verbosity not in ['NONE', 'LOW', 'HIGH', 'ALL']:
         raise ValueError('Verbosity should be either NONE, LOW, HIGH or ALL.')
@@ -243,12 +248,14 @@ def solve(functions, x0, solver=None, atol=None, dtol=None, rtol=1e-3,
     tstart = time.time()
     crit = None
     niter = 0
-    objective = [[f.eval(x0) for f in functions]]
     rtol_only_zeros = True
 
     # Solver specific initialization.
     solver.pre(functions, x0)
 
+    # Evaluate the objective function at the begining
+    objective = [solver.objective(x0)]
+
     while not crit:
 
         niter += 1
@@ -263,7 +270,7 @@ def solve(functions, x0, solver=None, atol=None, dtol=None, rtol=1e-3,
         # Solver iterative algorithm.
         solver.algo(objective, niter)
 
-        objective.append([f.eval(solver.sol) for f in functions])
+        objective.append(solver.objective(solver.sol))
         current = np.sum(objective[-1])
         last = np.sum(objective[-2])
 
@@ -425,6 +432,19 @@ def post(self):
     def _post(self):
         raise NotImplementedError("Class user should define this method.")
 
+    def objective(self, x):
+        """
+        Return the objective function at x.
+
+        Necessitate `solver._pre(...)` to be run first.
+        """
+        return self._objective(x)
+
+    def _objective(self, x):
+        obj_smooth = [f.eval(x) for f in self.smooth_funs]
+        obj_nonsmooth = [f.eval(x) for f in self.non_smooth_funs]
+        return obj_nonsmooth + obj_smooth
+
 
 class gradient_descent(solver):
     r"""
@@ -780,14 +800,21 @@ def __init__(self, L=None, Lt=None, d0=None, *args, **kwargs):
     def _pre(self, functions, x0):
         # Dual variable.
         if self.d0 is None:
-            self.dual_sol = self.L(x0)
+            # The copy is necessary in case `L = lambda x: x`.
+            self.dual_sol = self.L(np.asarray(x0).copy())
         else:
             self.dual_sol = self.d0
 
     def _post(self):
         self.d0 = None
         del self.dual_sol
 
+    def _objective(self, x):
+        obj_smooth = [f.eval(x) for f in self.smooth_funs]
+        obj_nonsmooth = [self.non_smooth_funs[0].eval(x),
+                         self.non_smooth_funs[1].eval(self.L(x))]
+        return obj_nonsmooth + obj_smooth
+
 
 class mlfbf(primal_dual):
     r"""
@@ -831,7 +858,7 @@ class mlfbf(primal_dual):
     >>> solver = solvers.mlfbf(L=L, step=max_step/2.)
     >>> ret = solvers.solve([f, g, h], x0, solver, maxit=1000, rtol=0)
     Solution found after 1000 iterations:
-        objective function f(sol) = 1.833865e+05
+        objective function f(sol) = 1.839060e+05
         stopping criterion: MAXIT
     >>> ret['sol']
     array([1., 1., 1.])

pyunlocbox/tests/test_functions.py

@@ -600,5 +600,22 @@ def test_independent_problems(self):
                     res[:, iN] = f.grad(X[:, iN])
                 nptest.assert_array_almost_equal(res, f.grad(X))
 
+    def test_prox_star(self):
+        n = 10
+        x = 3 * np.random.randn(n, 1)
+        f = functions.norm_l1()
+        f2 = functions.dummy()
+        f2.prox = lambda x, T: functions._prox_star(f, x, T)
+        gamma = np.random.rand()
+
+        p1 = f.prox(x, gamma)
+        p2 = functions._prox_star(f2, x, gamma)
+
+        np.testing.assert_array_almost_equal(p1, p2)
+
+        p1 = f.prox(x, gamma) - x
+        p2 = -gamma * f2.prox(x / gamma, 1 / gamma)
+        np.testing.assert_array_almost_equal(p1, p2)
+
 
 suite = unittest.TestLoader().loadTestsFromTestCase(TestCase)

pyunlocbox/tests/test_solvers.py

@@ -297,6 +297,59 @@ def x0(): return np.zeros(len(x))
         # Sanity check
         self.assertRaises(ValueError, solver.pre, [f, g], x0())
 
+        # Make a second test where the solution is calculated by hand
+        n = 10
+        y = np.random.rand(n) * 2
+        z = np.random.rand(n)
+        c = 1
+
+        delta = (y - z - c)**2 + 4 * (1 + y * z - z * c)
+        sol = 0.5 * ((y - z - c) + np.sqrt(delta))
+
+        class mlog(functions.func):
+            def __init__(self, z):
+                super().__init__()
+                self.z = z
+
+            def _eval(self, x):
+                return -np.sum(np.log(x + self.z))
+
+            def _prox(self, x, T):
+                delta = (x - self.z)**2 + 4 * (T + x * self.z)
+                sol = 0.5 * (x - self.z + np.sqrt(delta))
+                return sol
+
+        f = functions.norm_l1(lambda_=c)
+        g = mlog(z=z)
+        h = functions.norm_l2(lambda_=0.5, y=y)
+
+        mu = 1 + 1
+        step = 1 / mu / 2
+
+        solver = solvers.mlfbf(step=step)
+        ret = solvers.solve([f, g, h],
+                            y.copy(),
+                            solver,
+                            maxit=200,
+                            rtol=0,
+                            verbosity="NONE")
+
+        nptest.assert_allclose(ret["sol"], sol, atol=1e-10)
+
+        # Make a final test where the function g can not be evaluate
+        # on the primal variables
+        y = np.random.rand(3)
+        y_2 = L.dot(y)
+        L = np.array([[5, 9, 3], [7, 8, 5], [4, 4, 9], [0, 1, 7]])
+        x0 = np.zeros(len(y))
+        f = functions.norm_l1(y=y)
+        g = functions.norm_l2(lambda_=0.5, y=y_2)
+        h = functions.norm_l2(y=y, lambda_=0.5)
+        max_step = 1 / (1 + np.linalg.norm(L, 2))
+        solver = solvers.mlfbf(L=L, step=max_step / 2.)
+        ret = solvers.solve([f, g, h], x0, solver, maxit=1000, rtol=0)
+        np.testing.assert_allclose(ret["sol"], y)
+
     def test_projection_based(self):
         """
         Test the projection-based solver with arbitrarily selected functions.
@@ -354,7 +407,12 @@ def test_solver_comparison(self):
             ret = solvers.solve([f1, f2], x0, solver, **params)
             nptest.assert_allclose(ret['sol'], sol)
             self.assertEqual(ret['niter'], niter)
-            self.assertIs(ret['sol'], x0)  # The initial value was modified.
+            # The initial value not was modified.
+            np.testing.assert_array_equal(np.zeros(len(y)), x0)
+            ret = solvers.solve([f1, f2], x0, solver,
+                                inplace=True, **params)
+            # The initial value was modified.
+            self.assertIs(ret['sol'], x0)
 
     def test_primal_dual_solver_comparison(self):
         """
@@ -366,33 +424,45 @@ def test_primal_dual_solver_comparison(self):
         """
 
         # Convex functions.
-        y = np.array([294, 390, 361])
-        sol = [1., 1., 1.]
-        L = np.array([[5, 9, 3], [7, 8, 5], [4, 4, 9], [0, 1, 7]])
+        y = np.random.randn(3)
+        L = np.random.randn(4, 3)
+
+        sol = y
+        y2 = L.dot(y)
         f1 = functions.norm_l1(y=y)
-        f2 = functions.norm_l1()
+        f2 = functions.norm_l2(y=y2)
         f3 = functions.dummy()
 
         # Solvers.
         step = 0.5 / (1 + np.linalg.norm(L, 2))
         slvs = []
-        slvs.append(solvers.mlfbf(step=step))
-        slvs.append(solvers.projection_based(step=step))
+        slvs.append(solvers.mlfbf(step=step, L=L))
+        slvs.append(solvers.projection_based(step=step, L=L))
 
         # Compare solutions.
-        params = {'rtol': 0, 'verbosity': 'NONE', 'maxit': 50}
-        niters = [50, 50]
-        for solver, niter in zip(slvs, niters):
+        niter = 1000
+        params = {'rtol': 0, 'verbosity': 'NONE', 'maxit': niter}
+        for solver in slvs:
             x0 = np.zeros(len(y))
 
             if type(solver) is solvers.mlfbf:
                 ret = solvers.solve([f1, f2, f3], x0, solver, **params)
             else:
                 ret = solvers.solve([f1, f2], x0, solver, **params)
-
             nptest.assert_allclose(ret['sol'], sol)
             self.assertEqual(ret['niter'], niter)
-            self.assertIs(ret['sol'], x0)  # The initial value was modified.
+            # The initial value was not modified.
+            nptest.assert_array_equal(x0, np.zeros(len(y)))
+
+            if type(solver) is solvers.mlfbf:
+                ret = solvers.solve([f1, f2, f3], x0, solver,
+                                    inplace=True, **params)
+            else:
+                ret = solvers.solve([f1, f2], x0, solver,
+                                    inplace=True, **params)
+            # The initial value was modified.
+            self.assertIs(ret['sol'], x0)
+            nptest.assert_allclose(ret['sol'], sol)
 
 
 suite = unittest.TestLoader().loadTestsFromTestCase(TestCase)