Minmax block

PEPit/examples/minmax_optimization/minmax_prox_blockcoordinate.py

@@ -0,0 +1,150 @@
+from PEPit import PEP
+from PEPit.functions import BlockConvexConcaveFunction
+from PEPit.point import Point
+from PEPit.expression import Expression
+import numpy as np
+
+def wc_proximal_point(gamma, n, wrapper="cvxpy", solver=None, verbose=1):
+    """
+    Consider the minimization problem
+
+    .. math:: f_\\star \\triangleq \\min_x f(x),
+
+    where :math:`f` is closed, proper, and convex (and potentially non-smooth).
+
+    This code computes a worst-case guarantee for the **proximal point method** with step-size :math:`\\gamma`.
+    That is, it computes the smallest possible :math:`\\tau(n,\\gamma)` such that the guarantee
+
+    .. math:: f(x_n) - f_\\star \\leqslant \\tau(n, \\gamma)  \\|x_0 - x_\\star\\|^2
+
+    is valid, where :math:`x_n` is the output of the proximal point method, and where :math:`x_\\star` is a
+    minimizer of :math:`f`.
+
+    In short, for given values of :math:`n` and :math:`\\gamma`,
+    :math:`\\tau(n,\\gamma)` is computed as the worst-case value of :math:`f(x_n)-f_\\star`
+    when :math:`\\|x_0 - x_\\star\\|^2 \\leqslant 1`.
+
+    **Algorithm**:
+
+    The proximal point method is described by
+
+        .. math:: x_{t+1} = \\arg\\min_x \\left\\{f(x)+\\frac{1}{2\gamma}\\|x-x_t\\|^2 \\right\\},
+
+    where :math:`\\gamma` is a step-size.
+
+    **Theoretical guarantee**:
+
+    The **tight** theoretical guarantee can be found in [1, Theorem 4.1]:
+
+        .. math:: f(x_n)-f_\\star \\leqslant \\frac{\\|x_0-x_\\star\\|^2}{4\\gamma n},
+
+    where tightness is obtained on, e.g., one-dimensional linear problems on the positive orthant.
+
+    **References**:
+
+    `[1] A. Taylor, J. Hendrickx, F. Glineur (2017).
+    Exact worst-case performance of first-order methods for composite convex optimization.
+    SIAM Journal on Optimization, 27(3):1283–1313.
+    <https://arxiv.org/pdf/1512.07516.pdf>`_
+
+    Args:
+        gamma (float): 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_point(gamma=3, n=4, wrapper="cvxpy", solver=None, verbose=1)
+        (PEPit) Setting up the problem: size of the Gram matrix: 6x6
+        (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 1 function(s)
+        			Function 1 : Adding 20 scalar constraint(s) ...
+        			Function 1 : 20 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.020833335685730252
+        (PEPit) Primal feasibility check:
+        		The solver found a Gram matrix that is positive semi-definite up to an error of 3.626659005644299e-09
+        		All the primal scalar constraints are verified up to an error of 1.1386158081487519e-08
+        (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 3.0464297827437203e-08
+        (PEPit) Final upper bound (dual): 0.020833337068527292 and lower bound (primal example): 0.020833335685730252 
+        (PEPit) Duality gap: absolute: 1.382797040067052e-09 and relative: 6.637425042856655e-08
+        *** Example file: worst-case performance of proximal point method ***
+        	PEPit guarantee:	 f(x_n)-f_* <= 0.0208333 ||x_0 - x_*||^2
+        	Theoretical guarantee:	 f(x_n)-f_* <= 0.0208333 ||x_0 - x_*||^2
+
+    """
+
+    # Instantiate PEP
+    problem = PEP()
+
+    # Declare a convex function
+
+    block_partition = problem.declare_block_partition(d=2)
+    func = problem.declare_function(function_class=BlockConvexConcaveFunction, partition=block_partition)
+    # Start by defining its unique optimal point xs = x_* and corresponding function value fs = f_*
+    zs = func.stationary_point()
+
+    # Then define the starting point x0 of the algorithm
+    z0 = problem.set_initial_point()
+
+    # Set the initial constraint that is the distance between x0 and x^*
+    problem.set_initial_condition((z0 - zs) ** 2 <= 1)
+
+    # Run n steps of the proximal point method
+    z = z0
+    z_prev = z
+
+    for _ in range(n):
+
+
+         # Compute x from the docstring equation.
+         gz = Point()
+         fz = Expression()
+         z_prev = z
+         z = z - gamma * block_partition.get_block(gz, 0)
+         z = z + gamma * block_partition.get_block(gz, 1)
+
+         # Add point to Function f.
+         func.add_point((z, gz, fz))
+
+
+    # Set the performance metric to the final distance to optimum in function values
+    problem.set_performance_metric((z-z_prev)**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 = 1 / (gamma * n)
+    theoretical_tau = np.power((1.0 - (1.0/n)), n-1)/n
+
+    # Print conclusion if required
+    if verbose != -1:
+        print('*** Example file: worst-case performance of proximal point method ***')
+        print('\tPEPit guarantee:\t ||z_n - z_n-1||^2 <= {:.6} ||z_0 - z_*||^2'.format(pepit_tau))
+        print('\tTheoretical guarantee:\t||z_n - z_n-1||^2<= {:.6} ||z_0 - z_*||^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_point(gamma=2, n=5, wrapper="cvxpy", solver=None, verbose=1)

PEPit/functions/__init__.py

@@ -1,3 +1,4 @@
+from .block_convex_concave_function import BlockConvexConcaveFunction
 from .block_smooth_convex_function_cheap import BlockSmoothConvexFunctionCheap
 from .block_smooth_convex_function_expensive import BlockSmoothConvexFunctionExpensive
 from .convex_function import ConvexFunction
@@ -15,7 +16,8 @@
 from .smooth_strongly_convex_quadratic_function import SmoothStronglyConvexQuadraticFunction
 from .strongly_convex_function import StronglyConvexFunction
 
-__all__ = ['block_smooth_convex_function_cheap', 'BlockSmoothConvexFunctionCheap',
+__all__ = ['block_convex_concave_function', 'BlockConvexConcaveFunction',
+           'block_smooth_convex_function_cheap', 'BlockSmoothConvexFunctionCheap',
            'block_smooth_convex_function_expensive', 'BlockSmoothConvexFunctionExpensive',
            'convex_function', 'ConvexFunction',
            'convex_indicator', 'ConvexIndicatorFunction',

PEPit/functions/block_convex_concave_function.py

@@ -0,0 +1,116 @@
+import numpy as np
+
+from PEPit.function import Function
+from PEPit.block_partition import BlockPartition
+
+
+class BlockConvexConcaveFunction(Function):
+    """
+    The :class:`BlockSmoothConvexFunction` class overwrites the `add_class_constraints` method of :class:`Function`,
+    by implementing necessary constraints for interpolation of the class of smooth convex functions by blocks.
+
+    Attributes:
+        partition (BlockPartition): partitioning of the variables (in blocks).
+        L (list): smoothness parameters (one per block).
+
+    Smooth convex functions by blocks are characterized by a list of parameters :math:`L_i` (one per block),
+    hence can be instantiated as
+
+    Example:
+        >>> from PEPit import PEP
+        >>> from PEPit.functions import BlockSmoothConvexFunction
+        >>> problem = PEP()
+        >>> partition = problem.declare_block_partition(d=3)
+        >>> L = [1, 4, 10]
+        >>> func = problem.declare_function(function_class=BlockSmoothConvexFunction, partition=partition, L=L)
+
+    References:
+
+    `[1] Z. Shi, R. Liu (2016).
+    Better worst-case complexity analysis of the block coordinate descent method for large scale machine learning.
+    In 2017 16th IEEE International Conference on Machine Learning and Applications (ICMLA).
+    <https://arxiv.org/pdf/1608.04826.pdf>`_
+
+    """
+
+    def __init__(self,
+                 partition,
+                 is_leaf=True,
+                 decomposition_dict=None,
+                 reuse_gradient=True,
+                 name=None):
+        """
+
+        Args:
+            partition (BlockPartition): a :class:`BlockPartition`.
+            is_leaf (bool): True if self is defined from scratch.
+                            False if self is defined as linear combination of leaf.
+            decomposition_dict (dict): Decomposition of self as linear combination of leaf :class:`Function` objects.
+                                       Keys are :class:`Function` objects and values are their associated coefficients.
+            reuse_gradient (bool): If True, the same subgradient is returned
+                                   when one requires it several times on the same :class:`Point`.
+                                   If False, a new subgradient is computed each time one is required.
+            name (str): name of the object. None by default. Can be updated later through the method `set_name`.
+
+        Note:
+            Smooth convex functions by blocks are necessarily differentiable, hence `reuse_gradient` is set to True.
+
+        """
+        super().__init__(is_leaf=is_leaf,
+                         decomposition_dict=decomposition_dict,
+                         reuse_gradient=True,
+                         name=name,
+                         )
+
+        # Store partition and L
+        assert isinstance(partition, BlockPartition)
+        assert partition.get_nb_blocks() == 2
+        self.partition = partition
+
+
+    def add_class_constraints(self):
+        """
+        Formulates the list of necessary constraints for interpolation for self (block smooth convex function);
+        see [1, Lemma 1.1].
+
+        """
+        # Set function ID
+        function_id = self.get_name()
+        if function_id is None:
+            function_id = "Function_{}".format(self.counter)
+
+        # Set tables_of_constraints attributes
+        self.tables_of_constraints["convexity_concavity_block_{}".format(0)] = [[]]*len(self.list_of_points)
+
+        # Browse list of points and create interpolation constraints
+        for i, point_i in enumerate(self.list_of_points):
+
+            zi, gi, fi = point_i
+            zi_id = zi.get_name()
+            if zi_id is None:
+                zi_id = "Point_{}".format(i)
+
+            for j, point_j in enumerate(self.list_of_points):
+
+                zj, gj, fj = point_j
+                zj_id = zj.get_name()
+                if zj_id is None:
+                    zj_id = "Point_{}".format(j)
+
+                if point_i == point_j:
+                    self.tables_of_constraints["convexity_concavity_block_{}".format(0)][i].append(0)
+
+                else:
+                    gj_x = self.partition.get_block(gj, 0)
+                    xi = self.partition.get_block(zi, 0)
+                    xj = self.partition.get_block(zj, 0)
+
+                    gi_y = self.partition.get_block(gi, 1)
+                    yi = self.partition.get_block(zi, 1)
+                    yj = self.partition.get_block(zj, 1)
+                    constraint = fi >= fj + gj_x * (xi - xj) + gi_y * (yi - yj)
+                    constraint.set_name("IC_{}_convexity_concavity_block_{}({}, {})".format(function_id, 0,
+                                                                                  zi_id, zj_id))
+                    self.tables_of_constraints["convexity_concavity_block_{}".format(0)][i].append(constraint)
+                    self.list_of_class_constraints.append(constraint)
+