Merge pull request #118 from PerformanceEstimation/feature/radius_for_convex_indicator

Feature/radius for convex indicator

Author: bgoujaud

PEPit/examples/composite_convex_minimization/frank_wolfe.py

@@ -4,24 +4,25 @@
 from PEPit.primitive_steps import linear_optimization_step
 
 
-def wc_frank_wolfe(L, D, n, wrapper="cvxpy", solver=None, verbose=1):
+def wc_frank_wolfe(L, D, R, center, 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 and convex
-    and where :math:`f_2` is a convex indicator function on :math:`\\mathcal{D}` of diameter at most :math:`D`.
+    and where :math:`f_2` is a convex indicator function on :math:`\\mathcal{D}`
+    of diameter at most :math:`D` and radius at most :math:`R` around `center`.
 
     This code computes a worst-case guarantee for the **conditional gradient** method, aka **Frank-Wolfe** method.
     That is, it computes the smallest possible :math:`\\tau(n, L)` such that the guarantee
 
-    .. math :: F(x_n) - F(x_\\star) \\leqslant \\tau(n, L) D^2,
+    .. math :: F(x_n) - F(x_\\star) \\leqslant \\tau(n, L, D, R),
 
     is valid, where x_n is the output of the **conditional gradient** method,
     and where :math:`x_\\star` is a minimizer of :math:`F`.
-    In short, for given values of :math:`n` and :math:`L`, :math:`\\tau(n, L)` is computed as the worst-case value of
-    :math:`F(x_n) - F(x_\\star)` when :math:`D \\leqslant 1`.
+    In short, for given values of :math:`n` and :math:`L`, :math:`\\tau(n, L, D, R)`
+    is computed as the worst-case value of :math:`F(x_n) - F(x_\\star)`.
 
     **Algorithm**:
 
@@ -36,7 +37,7 @@ def wc_frank_wolfe(L, D, n, wrapper="cvxpy", solver=None, verbose=1):
 
     **Theoretical guarantee**:
 
-    An **upper** guarantee obtained in [2, Theorem 1] is
+    An **upper** guarantee obtained in [2, Theorem 1] when R = infinity is
 
         .. math :: F(x_n) - F(x_\\star) \\leqslant \\frac{2L D^2}{n+2}.
 
@@ -54,7 +55,9 @@ def wc_frank_wolfe(L, D, n, wrapper="cvxpy", solver=None, verbose=1):
 
     Args:
         L (float): the smoothness parameter.
-        D (float): diameter of :math:`f_2`.
+        D (float): diameter of :math:`\\mathcal{D}`.
+        R (float): radius of :math:`\\mathcal{D}`.
+        center (Point): center of :math:`\\mathcal{D}`. If None, the radius constraint must be observed to one center.
         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.
@@ -70,7 +73,8 @@ def wc_frank_wolfe(L, D, n, wrapper="cvxpy", solver=None, verbose=1):
         theoretical_tau (float): theoretical value.
 
     Example:
-        >>> pepit_tau, theoretical_tau = wc_frank_wolfe(L=1, D=1, n=10, wrapper="cvxpy", solver=None, verbose=1)
+        >>> pepit_tau, theoretical_tau = wc_frank_wolfe(L=1, D=1, R=float('inf'), center=None, n=10,
+        >>>                                             wrapper="cvxpy", solver=None, verbose=1)
         (PEPit) Setting up the problem: size of the Gram matrix: 26x26
         (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 ...
@@ -102,9 +106,9 @@ def wc_frank_wolfe(L, D, n, wrapper="cvxpy", solver=None, verbose=1):
     # Instantiate PEP
     problem = PEP()
 
-    # Declare a smooth convex function and a convex indicator of rayon D
+    # Declare a smooth convex function and a convex indicator of diameter D
     func1 = problem.declare_function(function_class=SmoothConvexFunction, L=L)
-    func2 = problem.declare_function(function_class=ConvexIndicatorFunction, D=D)
+    func2 = problem.declare_function(function_class=ConvexIndicatorFunction, D=D, R=R, center=center)
     # Define the function to optimize as the sum of func1 and func2
     func = func1 + func2
 
@@ -135,7 +139,7 @@ def wc_frank_wolfe(L, D, n, wrapper="cvxpy", solver=None, verbose=1):
     pepit_tau = problem.solve(wrapper=wrapper, solver=solver, verbose=pepit_verbose)
 
     # Compute theoretical guarantee (for comparison)
-    # when theta = 1
+    # when theta = 1, and R = infinity
     theoretical_tau = 2 * L * D ** 2 / (n + 2)
 
     # Print conclusion if required
@@ -149,4 +153,5 @@ def wc_frank_wolfe(L, D, n, wrapper="cvxpy", solver=None, verbose=1):
 
 
 if __name__ == "__main__":
-    pepit_tau, theoretical_tau = wc_frank_wolfe(L=1, D=1, n=10, wrapper="cvxpy", solver=None, verbose=1)
+    pepit_tau, theoretical_tau = wc_frank_wolfe(L=1, D=1, R=float('inf'), center=None, n=10,
+                                                wrapper="cvxpy", solver=None, verbose=1)

PEPit/functions/convex_indicator.py

@@ -1,7 +1,9 @@
 import numpy as np
 
+from PEPit.point import Point
 from PEPit.function import Function
 
+
 class ConvexIndicatorFunction(Function):
     """
     The :class:`ConvexIndicatorFunction` class overwrites the `add_class_constraints` method of :class:`Function`,
@@ -10,8 +12,9 @@ class ConvexIndicatorFunction(Function):
     Attributes:
         D (float): upper bound on the diameter of the feasible set, possibly set to np.inf
         R (float): upper bound on the radius of the feasible set, possibly set to np.inf
-        center (Point): Center of the feasible set spanned by the radius constraint, possibly set to None.
-    Convex indicator functions are characterized by a parameter `D` (or `R`), hence can be instantiated as
+        center (Point): Center of the feasible set spanned by the radius constraint.
+                        If set to None, there exists such a point but it remains undefined.
+    Convex indicator functions are characterized by a parameter `D` and/or `R`, hence can be instantiated as
 
     Example:
         >>> from PEPit import PEP
@@ -67,8 +70,11 @@ def __init__(self,
         self.D = D
         # Store the radius R in an attribute
         self.R = R
-        # Store the center in an attribute 
-        self.center = center
+        # Store the center in an attribute
+        if center is None and self.R != np.inf:
+            self.center = Point()
+        else:
+            self.center = center
 
     @staticmethod
     def set_value_constraint_i(xi, gi, fi):
@@ -104,14 +110,17 @@ def set_diameter_constraint_i_j(self,
         # Diameter constraint
         constraint = ((xi - xj) ** 2 <= self.D ** 2)
 
-        # Radius constraint 
-        if self.R < np.inf:
-            # No self.center provided centers the ball on the origin
-            if self.center is None: 
-                constraint = ((xi)**2 <= self.R ** 2)
-            # Centering the ball on self.center
-            else: 
-                constraint = ((self.center - xi)**2 <= self.R ** 2)
+        return constraint
+
+    def set_radius_constraint_i(self,
+                                xi, gi, fi,
+                                ):
+        """
+        Formulate the constraints bounding the radius of the support of self around self.center.
+
+        """
+        # Radius constraint
+        constraint = ((xi - self.center)**2 <= self.R ** 2)
 
         return constraint
 
@@ -130,9 +139,16 @@ def add_class_constraints(self):
                                                       constraint_name="convexity",
                                                       set_class_constraint_i_j=self.set_convexity_constraint_i_j,
                                                       )
-        if (self.D != np.inf) or (self.R != np.inf):
+
+        if self.D != np.inf:
             self.add_constraints_from_two_lists_of_points(list_of_points_1=self.list_of_points,
                                                           list_of_points_2=self.list_of_points,
                                                           constraint_name="diameter",
                                                           set_class_constraint_i_j=self.set_diameter_constraint_i_j,
                                                           )
+
+        if self.R != np.inf:
+            self.add_constraints_from_one_list_of_points(list_of_points=self.list_of_points,
+                                                         constraint_name="radius",
+                                                         set_class_constraint_i=self.set_radius_constraint_i,
+                                                         )

docs/source/whatsnew/0.4.0.rst

@@ -14,3 +14,5 @@ What's new in PEPit 0.4.0
 - New examples: (i) Refined block-coordinate descent (refined interpolation techniques--but computationally more expensive) (ii) introduction of quadratic Lojasiewicz inequalities, examples on gradient descent (with different ways to improse Lojasiewicz inequalities, from naive ones to more advanced SDP-representable ones) (iii) optimistic gradient with more advanced SDP-representable monotonocity/Lipschitz characterization, (iv) difference-of-convex algorithm (DCA, also known as the convex-concave procedure),  (v) online learning settings (online gradient descent, online Frank-Wolfe, follow the leader, follow the regularized leader).
 
 - Corrected example: SGD's closed-form solution has been updated.
+
+- The support of functions of the class :class:`ConvexIndicatorFunction` can be constrained through a radius in addition to a diameter.