Merge branch 'master' into frugal-resolvent-splitting

PEPit/examples/composite_convex_minimization/__init__.py

@@ -8,6 +8,7 @@
 from .no_lips_in_function_value import wc_no_lips_in_function_value
 from .no_lips_in_bregman_divergence import wc_no_lips_in_bregman_divergence
 from .proximal_gradient import wc_proximal_gradient
+from .proximal_gradient_quadratics import wc_proximal_gradient_quadratics
 from .three_operator_splitting import wc_three_operator_splitting
 
 __all__ = ['accelerated_douglas_rachford_splitting', 'wc_accelerated_douglas_rachford_splitting',
@@ -20,5 +21,6 @@
            'no_lips_in_function_value', 'wc_no_lips_in_function_value',
            'no_lips_in_bregman_divergence', 'wc_no_lips_in_bregman_divergence',
            'proximal_gradient', 'wc_proximal_gradient',
+           'proximal_gradient_quadratics', 'wc_proximal_gradient_quadratics',
            'three_operator_splitting', 'wc_three_operator_splitting',
            ]

PEPit/examples/composite_convex_minimization/proximal_gradient_quadratics.py

@@ -0,0 +1,158 @@
+from PEPit import PEP
+from PEPit.functions import ConvexFunction
+from PEPit.functions import SmoothStronglyConvexQuadraticFunction
+from PEPit.primitive_steps import proximal_step
+
+
+def wc_proximal_gradient_quadratics(L, mu, gamma, n, wrapper="cvxpy", solver=None, verbose=1):
+    """
+    Consider the composite convex minimization problem
+
+    .. math:: F_\\star \\triangleq \\min_x \\{F(x) \\equiv f_1(x) + f_2(x)\\},
+
+    where :math:`f_1` is :math:`L`-smooth, :math:`\\mu`-strongly convex and quadratic,
+    and where :math:`f_2` is closed convex and proper.
+
+    This code computes a worst-case guarantee for the **proximal gradient** method (PGM).
+    That is, it computes the smallest possible :math:`\\tau(n, L, \\mu)` such that the guarantee
+
+    .. math :: \\|x_n - x_\\star\\|^2 \\leqslant \\tau(n, L, \\mu) \\|x_0 - x_\\star\\|^2,
+
+    is valid, where :math:`x_n` is the output of the **proximal gradient**,
+    and where :math:`x_\\star` is a minimizer of :math:`F`.
+    In short, for given values of :math:`n`, :math:`L` and :math:`\\mu`,
+    :math:`\\tau(n, L, \\mu)` is computed as the worst-case value of
+    :math:`\\|x_n - x_\\star\\|^2` when :math:`\\|x_0 - x_\\star\\|^2 \\leqslant 1`.
+
+    **Algorithm**: Proximal gradient is described by
+
+        .. math::
+            \\begin{eqnarray}
+                y_t & = & x_t - \\gamma \\nabla f_1(x_t), \\\\
+                x_{t+1} & = & \\arg\\min_x \\left\\{f_2(x)+\\frac{1}{2\gamma}\|x-y_t\|^2 \\right\\},
+            \\end{eqnarray}
+
+    for :math:`t \in \\{ 0, \\dots, n-1\\}` and where :math:`\\gamma` is a step-size.
+
+    **Theoretical guarantee**: It is well known that a **tight** guarantee for PGM is provided by
+
+    .. math :: \\|x_n - x_\\star\\|^2 \\leqslant \\max\\{(1-L\\gamma)^2,(1-\\mu\\gamma)^2\\}^n \\|x_0 - x_\\star\\|^2,
+
+    which can be found in, e.g., [1, Theorem 3.1]. It is a folk knowledge and the result can be found in many references
+    for gradient descent; see, e.g.,[2, Section 1.4: Theorem 3], [3, Section 5.1] and [4, Section 4.4].
+
+    **References**:
+
+    `[1] A. Taylor, J. Hendrickx, F. Glineur (2018).
+    Exact worst-case convergence rates of the proximal gradient method for composite convex minimization.
+    Journal of Optimization Theory and Applications, 178(2), 455-476.
+    <https://arxiv.org/pdf/1705.04398.pdf>`_
+
+    [2] B. Polyak (1987).
+    Introduction to Optimization.
+    Optimization Software New York.
+
+    `[3] E. Ryu, S. Boyd (2016).
+    A primer on monotone operator methods.
+    Applied and Computational Mathematics 15(1), 3-43.
+    <https://web.stanford.edu/~boyd/papers/pdf/monotone_primer.pdf>`_
+
+    `[4] L. Lessard, B. Recht, A. Packard (2016).
+    Analysis and design of optimization algorithms via integral quadratic constraints.
+    SIAM Journal on Optimization 26(1), 57–95.
+    <https://arxiv.org/pdf/1408.3595.pdf>`_
+
+    Args:
+        L (float): the smoothness parameter.
+        mu (float): the strong convexity parameter.
+        gamma (float): proximal step-size.
+        n (int): number of iterations.
+        wrapper (str): the name of the wrapper to be used.
+        solver (str): the name of the solver the wrapper should use.
+        verbose (int): level of information details to print.
+
+                        - -1: No verbose at all.
+                        - 0: This example's output.
+                        - 1: This example's output + PEPit information.
+                        - 2: This example's output + PEPit information + solver details.
+
+    Returns:
+        pepit_tau (float): worst-case value.
+        theoretical_tau (float): theoretical value.
+
+    Example:
+        >>> pepit_tau, theoretical_tau = wc_proximal_gradient_quadratics(L=1, mu=.1, gamma=1, n=2, wrapper="cvxpy", solver=None, verbose=1)
+        (PEPit) Setting up the problem: size of the Gram matrix: 8x8
+        (PEPit) Setting up the problem: performance measure is the minimum of 1 element(s)
+        (PEPit) Setting up the problem: Adding initial conditions and general constraints ...
+        (PEPit) Setting up the problem: initial conditions and general constraints (1 constraint(s) added)
+        (PEPit) Setting up the problem: interpolation conditions for 2 function(s)
+        			Function 1 : Adding 22 scalar constraint(s) ...
+        			Function 1 : 22 scalar constraint(s) added
+        			Function 2 : Adding 6 scalar constraint(s) ...
+        			Function 2 : 6 scalar constraint(s) added
+        (PEPit) Setting up the problem: additional constraints for 0 function(s)
+        (PEPit) Compiling SDP
+        (PEPit) Calling SDP solver
+        (PEPit) Solver status: optimal (wrapper:cvxpy, solver: MOSEK); optimal value: 0.6561000187100321
+        (PEPit) Primal feasibility check:
+        		The solver found a Gram matrix that is positive semi-definite up to an error of 3.8506403023071055e-09
+        		All the primal scalar constraints are verified up to an error of 5.880885747128195e-09
+        (PEPit) Dual feasibility check:
+        		The solver found a residual matrix that is positive semi-definite
+        		All the dual scalar values associated with inequality constraints are nonnegative
+        (PEPit) The worst-case guarantee proof is perfectly reconstituted up to an error of 1.8758326292663516e-07
+        (PEPit) Final upper bound (dual): 0.6561000176340664 and lower bound (primal example): 0.6561000187100321 
+        (PEPit) Duality gap: absolute: -1.0759656499104153e-09 and relative: -1.6399415016416052e-09
+        *** Example file: worst-case performance of the Proximal Gradient Method in function values***
+        	PEPit guarantee:	 ||x_n - x_*||^2 <= 0.6561 ||x0 - xs||^2
+        	Theoretical guarantee:	 ||x_n - x_*||^2 <= 0.6561 ||x0 - xs||^2
+    
+    """
+
+    # Instantiate PEP
+    problem = PEP()
+
+    # Declare a strongly convex smooth function and a closed convex proper function
+    f1 = problem.declare_function(SmoothStronglyConvexQuadraticFunction, mu=mu, L=L)
+    f2 = problem.declare_function(ConvexFunction)
+    func = f1 + f2
+
+    # Start by defining its unique optimal point xs = x_*
+    xs = func.stationary_point()
+
+    # Then define the starting point x0 of the algorithm
+    x0 = problem.set_initial_point()
+
+    # Set the initial constraint that is the distance between x0 and x^*
+    problem.set_initial_condition((x0 - xs) ** 2 <= 1)
+
+    # Run the proximal gradient method starting from x0
+    x = x0
+    for _ in range(n):
+        y = x - gamma * f1.gradient(x)
+        x, _, _ = proximal_step(y, f2, gamma)
+
+    # Set the performance metric to the distance between x and xs
+    problem.set_performance_metric((x - xs) ** 2)
+
+    # Solve the PEP
+    pepit_verbose = max(verbose, 0)
+    pepit_tau = problem.solve(wrapper=wrapper, solver=solver, verbose=pepit_verbose)
+
+    # Compute theoretical guarantee (for comparison)
+    theoretical_tau = max((1 - mu * gamma) ** 2, (1 - L * gamma) ** 2) ** n
+
+    # Print conclusion if required
+    if verbose != -1:
+        print('*** Example file: worst-case performance of the Proximal Gradient Method in function values***')
+        print('\tPEPit guarantee:\t ||x_n - x_*||^2 <= {:.6} ||x0 - xs||^2'.format(pepit_tau))
+        print('\tTheoretical guarantee:\t ||x_n - x_*||^2 <= {:.6} ||x0 - xs||^2'.format(theoretical_tau))
+
+    # Return the worst-case guarantee of the evaluated method ( and the reference theoretical value)
+    return pepit_tau, theoretical_tau
+
+
+if __name__ == "__main__":
+    pepit_tau, theoretical_tau = wc_proximal_gradient_quadratics(L=1, mu=.1, gamma=1, n=2,
+                                                                 wrapper="cvxpy", solver=None, verbose=1)

PEPit/functions/smooth_strongly_convex_quadratic_function.py

@@ -67,17 +67,44 @@ def __init__(self,
         self.mu = mu
         self.L = L
 
+        # Create the unique stationary point
+        super().stationary_point()
+
+    def stationary_point(self, return_gradient_and_function_value=False, name=None):
+        """
+        Access the unique stationary point that has been created in __init__,
+        as well as its zero gradient and its function value.
+
+        Args:
+            return_gradient_and_function_value (bool): if True, return the triplet point (:class:`Point`),
+                                                       gradient (:class:`Point`), function value (:class:`Expression`).
+                                                       Otherwise, return only the point (:class:`Point`).
+            name (str, optional): name of the object. Unused since no object is created.
+
+        Returns:
+            Point or tuple: the minimizer
+
+        """
+        # Create a new point, null gradient and new function value
+        point, g, f = self.list_of_stationary_points[0]
+
+        # Return the required information
+        if return_gradient_and_function_value:
+            return point, g, f
+        else:
+            return point
+
     def set_value_constraint_i(self,
                                xi, gi, fi):
         """
         Set the value of the function.
 
         """
         # Select one stationary point
-        xs = self.list_of_stationary_points[0][0]
+        xs, _, fs = self.list_of_stationary_points[0]
 
         # Value constraint
-        constraint = (fi == 0.5 * (xi - xs) * gi)
+        constraint = (fi - fs == 0.5 * (xi - xs) * gi)
 
         return constraint
 
@@ -101,10 +128,6 @@ def add_class_constraints(self):
         Formulates the list of interpolation constraints for self (smooth strongly convex quadratic function);
         see [1, Theorem 3.9].
         """
-        # Create a stationary point is none exists
-        if self.list_of_stationary_points == list():
-            self.stationary_point()
-
         # Add the quadratic interpolation constraint
         self.add_constraints_from_one_list_of_points(list_of_points=self.list_of_points,
                                                      constraint_name="value",

docs/source/examples/b.rst

@@ -11,6 +11,11 @@ Proximal gradient
 .. autofunction:: PEPit.examples.composite_convex_minimization.wc_proximal_gradient
 
 
+Proximal gradient on quadratics
+-------------------------------
+.. autofunction:: PEPit.examples.composite_convex_minimization.wc_proximal_gradient_quadratics
+
+
 Accelerated proximal gradient
 -----------------------------
 .. autofunction:: PEPit.examples.composite_convex_minimization.wc_accelerated_proximal_gradient

docs/source/whatsnew/0.3.3.rst

@@ -1,9 +1,10 @@
-What's new in PEPit 0.3.2
+What's new in PEPit 0.3.3
 =========================
 
-For users:
-----------
+- Silver step-sizes have been added in the `PEPit.examples.unconstrained_convex_minimization <https://pepit.readthedocs.io/en/0.3.2/examples/a.html>`_ module.
 
-- Silver step-sizes have been added in the PEPit.examples.unconstrained_convex_minimization module.
+- An monotone inclusion example has been added in the `PEPit.examples.monotone_inclusions_variational_inequalities <https://pepit.readthedocs.io/en/0.3.2/examples/e.html>`_ module. It covers parameterized frugal resolvent splittings, including the Malitsky-Tam algorithm, Douglas-Rachford, the Ryu Three Operator Splitting, and a set of novel block splitting methods. This example compares a fixed step size for this class with one optimized using the dual of the PEP.
 
-- An monotone inclusion example has been added which covers parameterized frugal resolvent splittings, including the Malitsky-Tam algorithm, Douglas-Rachford, the Ryu Three Operator Splitting, and a set of novel block splitting methods. This example compares a fixed step size for this class with one optimized using the dual of the PEP.
+- A fix has been added in the class :class:`SmoothStronglyConvexQuadraticFunction`. Prior to that fix, using this class without stationary point in a PEP solved with the direct interface to MOSEK was problematic due to the late creation of a stationary point. After this fix, a stationary point is automatically created when instantiating this class of functions.
+
+- Another modification has been made to the class :class:`SmoothStronglyConvexQuadraticFunction`. Prior to that, the minimum value was assumed to be 0. This is not the case anymore.

setup.py

@@ -3,7 +3,7 @@
 with open("README.md", "r", encoding="utf-8") as fh:
     long_description = fh.read()
 
-version = "0.3.2"
+version = "0.3.3"
 
 setuptools.setup(
     name="PEPit",

tests/test_examples.py

@@ -40,6 +40,7 @@
 from PEPit.examples.composite_convex_minimization import wc_no_lips_in_bregman_divergence
 from PEPit.examples.composite_convex_minimization import wc_no_lips_in_function_value
 from PEPit.examples.composite_convex_minimization import wc_proximal_gradient
+from PEPit.examples.composite_convex_minimization import wc_proximal_gradient_quadratics
 from PEPit.examples.composite_convex_minimization import wc_three_operator_splitting
 from PEPit.examples.nonconvex_optimization import wc_gradient_descent as wc_gradient_descent_non_convex
 from PEPit.examples.nonconvex_optimization import wc_no_lips_1
@@ -451,6 +452,12 @@ def test_proximal_gradient(self):
         wc, theory = wc_proximal_gradient(L=L, mu=mu, gamma=gamma, n=n, wrapper=self.wrapper, verbose=self.verbose)
         self.assertAlmostEqual(wc, theory, delta=self.relative_precision * theory)
 
+    def test_proximal_gradient_quadratics(self):
+        L, mu, gamma, n = 1, .1, 1, 2
+
+        wc, theory = wc_proximal_gradient_quadratics(L=L, mu=mu, gamma=gamma, n=n, wrapper=self.wrapper, verbose=self.verbose)
+        self.assertAlmostEqual(wc, theory, delta=self.relative_precision * theory)
+
     def test_three_operator_splitting(self):
         mu, L1, L3, alpha, theta = 0.1, 10, 1, 1, 1
         n_list = range(1, 4)

tests/test_functions_and_operators.py

@@ -63,6 +63,7 @@ def setUp(self):
         self.operator10 = SymmetricLinearOperator(mu=.1, L=1, name="op10")
         self.operator11 = SkewSymmetricLinearOperator(L=1, name="op11")
 
+        self.point0 = self.func12.list_of_stationary_points[0][0]
         self.point1 = Point(is_leaf=True, decomposition_dict=None, name="pt1")
         self.point2 = Point(is_leaf=True, decomposition_dict=None, name="pt2")
 
@@ -93,6 +94,7 @@ def setUp(self):
         ]
 
         self.all_points = [
+            self.point0,
             self.point1,
             self.point2,
         ]
@@ -136,8 +138,8 @@ def test_is_instance(self):
         self.assertIsInstance(self.operator10, SymmetricLinearOperator)
         self.assertIsInstance(self.operator11, SkewSymmetricLinearOperator)
 
-        self.assertIsInstance(self.point1, Point)
-        self.assertIsInstance(self.point2, Point)
+        for point in self.all_points:
+            self.assertIsInstance(point, Point)
 
         self.assertIsInstance(self.new_operator, Function)
         self.assertIsInstance(self.new_point, Point)
@@ -172,7 +174,7 @@ def test_add_constraints(self):
         self.assertEqual(len(self.new_operator.list_of_class_constraints), 0)
 
         # Call oracles on combination of functions
-        for point in self.all_points:
+        for point in self.all_points[1:]:
             self.new_operator.oracle(point)
         self.new_operator.oracle(self.new_point)
         self.new_operator.stationary_point()
@@ -182,11 +184,12 @@ def test_add_constraints(self):
 
         # Call stationary point on all functions
         for function in self.all_functions_and_operators:
-            function.stationary_point()
+            if not function.list_of_stationary_points:
+                function.stationary_point()
 
         # Set the number of points each function has been evaluated on.
         # Note self.new_operator has been evaluated on 1 point less than the other functions.
-        num_points_eval = len(self.all_points) + 3
+        num_points_eval = len(self.all_points[1:]) + 3
 
         # Add function class constraints
         for function in self.all_functions_and_operators:
@@ -213,6 +216,7 @@ def test_add_constraints(self):
         self.assertEqual(len(self.func10.list_of_class_constraints), num_points_eval * (num_points_eval - 1))
         self.assertEqual(len(self.func11.list_of_class_constraints), num_points_eval * (num_points_eval - 1))
         self.assertEqual(len(self.func12.list_of_class_constraints), num_points_eval * (num_points_eval + 1) / 2)
+
         self.assertEqual(len(self.operator1.list_of_class_constraints), num_points_eval * (num_points_eval - 1) / 2)
         self.assertEqual(len(self.operator2.list_of_class_constraints), num_points_eval * (num_points_eval - 1))
         self.assertEqual(len(self.operator3.list_of_class_constraints), num_points_eval * (num_points_eval - 1) / 2)