diff --git a/PEPit/examples/minmax_optimization/minmax_prox_blockcoordinate.py b/PEPit/examples/minmax_optimization/minmax_prox_blockcoordinate.py new file mode 100644 index 00000000..4df89bcd --- /dev/null +++ b/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. + `_ + + 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) \ No newline at end of file diff --git a/PEPit/functions/__init__.py b/PEPit/functions/__init__.py index 75ba72f9..bc211753 100644 --- a/PEPit/functions/__init__.py +++ b/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', diff --git a/PEPit/functions/block_convex_concave_function.py b/PEPit/functions/block_convex_concave_function.py new file mode 100644 index 00000000..a3a5c81e --- /dev/null +++ b/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). + `_ + + """ + + 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) + \ No newline at end of file