From 61c1b52b7b45b7eefd402ff24d49a2a793aed62f Mon Sep 17 00:00:00 2001 From: Peter Barkley Date: Tue, 6 Aug 2024 15:21:21 -0700 Subject: [PATCH 1/6] Add parameterized frugal resolvent splitting example --- .../__init__.py | 4 + .../frugal_resolvent_splitting.py | 277 ++++++++++++++++++ .../reduced_frugal_resolvent_splitting.py | 253 ++++++++++++++++ docs/source/examples/e.rst | 6 + docs/source/whatsnew/0.3.3.rst | 2 + 5 files changed, 542 insertions(+) create mode 100644 PEPit/examples/monotone_inclusions_variational_inequalities/frugal_resolvent_splitting.py create mode 100644 PEPit/examples/monotone_inclusions_variational_inequalities/reduced_frugal_resolvent_splitting.py diff --git a/PEPit/examples/monotone_inclusions_variational_inequalities/__init__.py b/PEPit/examples/monotone_inclusions_variational_inequalities/__init__.py index fe73a65e..eb4c1f3d 100644 --- a/PEPit/examples/monotone_inclusions_variational_inequalities/__init__.py +++ b/PEPit/examples/monotone_inclusions_variational_inequalities/__init__.py @@ -5,6 +5,8 @@ from .past_extragradient import wc_past_extragradient from .proximal_point import wc_proximal_point from .three_operator_splitting import wc_three_operator_splitting +from .frugal_resolvent_splitting import wc_frugal_resolvent_splitting +from .reduced_frugal_resolvent_splitting import wc_reduced_frugal_resolvent_splitting __all__ = ['accelerated_proximal_point', 'wc_accelerated_proximal_point', 'douglas_rachford_splitting', 'wc_douglas_rachford_splitting', @@ -13,4 +15,6 @@ 'past_extragradient', 'wc_past_extragradient', 'proximal_point', 'wc_proximal_point', 'three_operator_splitting', 'wc_three_operator_splitting', + 'frugal_resolvent_splitting', 'wc_frugal_resolvent_splitting', + 'reduced_frugal_resolvent_splitting', 'wc_reduced_frugal_resolvent_splitting' ] diff --git a/PEPit/examples/monotone_inclusions_variational_inequalities/frugal_resolvent_splitting.py b/PEPit/examples/monotone_inclusions_variational_inequalities/frugal_resolvent_splitting.py new file mode 100644 index 00000000..ad50a946 --- /dev/null +++ b/PEPit/examples/monotone_inclusions_variational_inequalities/frugal_resolvent_splitting.py @@ -0,0 +1,277 @@ +from PEPit import PEP, null_point, Constraint +from PEPit.primitive_steps import proximal_step +from PEPit.functions import SmoothStronglyConvexFunction +from numpy import array + +def wc_frugal_resolvent_splitting(L, W, lipschitz_values, mu_values, alpha=1, gamma=0.5, wrapper="cvxpy", solver=None, verbose=1): + """ + Consider the the monotone inclusion problem + + .. math:: \\mathrm{Find}\\, x:\\, 0 \\in \\sum_{i=1}^{n} A_i(x), + + where :math:`A_i` is the subdifferential of an :math:`l_i`-Lipschitz smooth and :math:`\\mu_i`-strongly convex function for all :math:`i \\leq n`. + We denote by :math:`J_{\\alpha A_i}` the resolvent of :math:`\\alpha A_i`. + We denote the lifted vector operator :math:`\\mathbf{A}` as :math:`\\mathbf{A} = [A_1, \\dots, A_n]`, + and use lifted :math:`\\mathbf{x} = [x_1, \\dots, x_n]` and :math:`\\mathbf{v} = [v_1, \\dots, v_n]`. + We denote by :math:`L, W \\in \mathbb{R}^{n \\times n}` the algorithm design matrices, and by :math:`l` and :math:`\\mu` the vectors of Lipschitz and strong convexity constants of the lifted operator :math:`\\mathbf{A}`. + :math:`L` is assumed to be strictly lower diagonal. + + This code computes a worst-case guarantee for any frugal resolvent splitting with design matrices :math:`L, W`. + As shown in [1] and [2], this can include the Malitsky-Tam [3], Ryu Three Operator Splitting [4], Douglas-Rachford [5], or block splitting algorithms [1]. + That is, given two lifted initial points :math:`\\mathbf{v}^{(0)}_t` and :math:`\\mathbf{v}^{(1)}_t` (each of which sums to 0), + this code computes the smallest possible :math:`\\tau(L, W, l, \\mu, \\alpha, \\gamma)` + (a.k.a. "contraction factor") such that the guarantee + + .. math:: \\|\\mathbf{v}^{(0)}_{t+1} - \\mathbf{v}^{(1)}_{t+1}\\|^2 \\leqslant \\tau(L, W, l, \\mu, \\alpha, \\gamma) \\|\\mathbf{v}^{(0)}_{t} - \\mathbf{v}^{(1)}_{t}\\|^2, + + is valid, where :math:`\\mathbf{v}^{(0)}_{t+1}` and :math:`\\mathbf{v}^{(1)}_{t+1}` are obtained after one iteration of the frugal resolvent splitting from respectively :math:`\\mathbf{v}^{(0)}_{t}` and :math:`\\mathbf{v}^{(1)}_{t}`. + + In short, for given values of :math:`L`, :math:`W`, :math:`l`, :math:`\\mu`, :math:`\\alpha` and :math:`\\gamma`, the contraction factor :math:`\\tau(L, W, \\mu, \\alpha, \\theta)` is computed as the worst-case value of + :math:`\\|\\mathbf{v}^{(0)}_{t+1} - \\mathbf{v}^{(1)}_{t+1}\\|^2` when :math:`\\|\\mathbf{v}^{(0)}_{t} - \\mathbf{v}^{(1)}_{t}\\|^2 \\leqslant 1`. + + **Algorithm**: One iteration of the parameterized frugal resolvent splitting is described as follows: + + .. math:: + :nowrap: + + \\begin{eqnarray} + \\mathbf{x}_{t+1} & = & J_{\\alpha \\mathbf{A}} (\\mathbf{L} \\mathbf{x}_{t+1} + \\mathbf{v}_t),\\\\ + \\mathbf{v}_{t+1} & = & \\mathbf{v}_t - \\gamma \\mathbf{W} \\mathbf{x}_{t+1}. + \\end{eqnarray} + + **References**: + + `[1] R. Bassett, P. Barkley (2024). + Optimal Design of Resolvent Splitting Algorithms. arxiv:2407.16159. + `_ + + `[2] M. Tam (2023). Frugal and decentralised resolvent splittings defined by nonexpansive operators. Optimization Letters pp 1–19. `_ + + `[3] Y. Malitsky, M. Tam (2023). Resolvent splitting for sums of monotone operators + with minimal lifting. Mathematical Programming 201(1-2):231–262. `_ + + `[4] E. Ryu (2020). Uniqueness of drs as the 2 operator resolvent-splitting and + impossibility of 3 operator resolvent-splitting. Mathematical Programming 182(1- + 2):233–273. `_ + + `[5] J. Eckstein, D. Bertsekas (1992). On the Douglas—Rachford splitting method and the proximal point algorithm for maximal monotone operators. Mathematical Programming 55:293–318. `_ + + + Args: + L (ndarray): n x n numpy array of resolvent multipliers for step 1. + W (ndarray): n x n numpy array of resolvent multipliers for step 2. + lipschitz_values (array): n Lipschitz parameters for the subdifferentials. + mu_values (array): n convexity parameters for the subdifferentials. + alpha (float): resolvent scaling parameter. + gamma (float): step size parameter. + 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 + + Example: + >>> pepit_tau = wc_frugal_resolvent_splitting( + L=array([[0,0],[2,0]]), + W=array([[1,-1],[-1,1]]), + lipschitz_values=[2, 1000], + mu_values=[1, 0], + alpha=1, + gamma=0.5, + wrapper="cvxpy", + 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 (3 constraint(s) added) + (PEPit) Setting up the problem: interpolation conditions for 2 function(s) + Function 1 : Adding 2 scalar constraint(s) ... + Function 1 : 2 scalar constraint(s) added + Function 2 : Adding 2 scalar constraint(s) ... + Function 2 : 2 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.6941881289170055 + (PEPit) Postprocessing: solver's output is not entirely feasible (smallest eigenvalue of the Gram matrix is: -3.01e-09 < 0). + Small deviation from 0 may simply be due to numerical error. Big ones should be deeply investigated. + In any case, from now the provided values of parameters are based on the projection of the Gram matrix onto the cone of symmetric semi-definite matrix. + (PEPit) Primal feasibility check: + The solver found a Gram matrix that is positive semi-definite up to an error of 3.0070570485427986e-09 + All the primal scalar constraints are verified up to an error of 7.876224117353559e-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 2.165773640630705e-06 + (PEPit) Final upper bound (dual): 0.6941851987362065 and lower bound (primal example): 0.6941881289170055 + (PEPit) Duality gap: absolute: -2.9301807989989825e-06 and relative: -4.2210183046061616e-06 + *** Example file: worst-case performance of parameterized frugal resolvent splitting` *** + PEPit guarantee: ||v_(t+1)^0 - v_(t+1)^1||^2 <= 0.694185 ||v_(t)^0 - v_(t)^1||^2 + `` + >>> comparison() + `` + Contraction factor of different designs with standard step size, optimal step size, and optimal W matrix + Optimized step sizes and W matrix from [4] when n=4 + Design 0.5 step size Optimal step size Optimal W matrix + --------------------------------------------------------------------- + MT 0.858 0.758 0.750 + Full 0.423 0.101 0.077 + Block 0.444 0.088 0.059 + `` + + """ + + # Instantiate PEP + problem = PEP() + + # Declare operators (works for LipschitzStronglyMonotoneOperator as well) + operators = [problem.declare_function(SmoothStronglyConvexFunction, L=l, mu=mu) for l, mu in zip(lipschitz_values, mu_values)] + + # Then define the starting points v0 and v1 + n = W.shape[0] + v0 = [problem.set_initial_point() for _ in range(n)] + v1 = [problem.set_initial_point() for _ in range(n)] + + # Set the initial constraint that is the distance between v0 and v1 + problem.set_initial_condition(sum((v0[i] - v1[i]) ** 2 for i in range(n)) <= 1) + + # Constraint on the lifted starting points so each sums to 0 + v0constraint = Constraint(expression=sum((v0[i] for i in range(n)), start=null_point)**2, equality_or_inequality="equality") + v1constraint = Constraint(expression=sum((v1[i] for i in range(n)), start=null_point)**2, equality_or_inequality="equality") + problem.set_initial_condition(v0constraint) + problem.set_initial_condition(v1constraint) + + # Define the step for each element of the lifted vector + def resolvent(i, x, v, L, alpha): + Lx = sum((L[i, j]*x[j] for j in range(i)), start=null_point) + x, _, _ = proximal_step(v[i] + Lx, operators[i], alpha) + return x + + x0 = [] + x1 = [] + for i in range(n): + x0.append(resolvent(i, x0, v0, L, alpha)) + x1.append(resolvent(i, x1, v1, L, alpha)) + + z0 = [] + z1 = [] + for i in range(n): + z0.append(v0[i] - gamma*W[i,:]@x0) + z1.append(v1[i] - gamma*W[i,:]@x1) + + + # Set the performance metric to the distance between z0 and z1 + problem.set_performance_metric(sum((z0[i] - z1[i]) ** 2 for i in range(n))) + + # Solve the PEP + pepit_verbose = max(verbose, 0) + pepit_tau = problem.solve(wrapper=wrapper, solver=solver, verbose=pepit_verbose) + + # Print conclusion if required + if verbose != -1: + print('*** Example file: worst-case performance of parameterized frugal resolvent splitting` ***') + print('\tPEPit guarantee:\t ||v_(t+1)^0 - v_(t+1)^1||^2 <= {:.6} ||v_(t)^0 - v_(t)^1||^2'.format(pepit_tau)) + + # Return the worst-case guarantee of the evaluated method (and the reference theoretical value) + return pepit_tau + +def comparison(): + # Comparison for 4 operators for Malitsky-Tam, Fully Connected, and 2-Block designs + # with and without optimized step sizes and W matrices + n = 4 + lipschitz_values = [2, 2, 2, 2] + mu_values = [1, 1, 1, 1] + + # Malitsky-Tam [3] + L_MT = array([[0,0,0,0], + [1,0,0,0], + [0,1,0,0], + [1,0,1,0]]) + W_MT = array([[1,-1,0,0], + [-1,2,-1,0], + [0,-1,2,-1], + [0,0,-1,1]]) + W_MT_opt = array([ + [ 1.071, -1.071, -0. , -0. ], + [-1.071, 2.071, -1. , -0. ], + [-0. , -1. , 2.5 , -1.5 ], + [-0. , -0. , -1.5 , 1.5 ]]) + + # Fully Connected + L_full = array([[0,0,0,0], + [2/3,0,0,0], + [2/3,2/3,0,0], + [2/3,2/3,2/3,0]]) + W_full = array([[2, -2/3, -2/3, -2/3], + [-2/3, 2, -2/3, -2/3], + [-2/3, -2/3, 2, -2/3], + [-2/3, -2/3, -2/3, 2]]) + W_full_opt = array([ + [ 2.226, -0.577, -0.714, -0.936], + [-0.577, 1.94 , -0.604, -0.759], + [-0.714, -0.604, 1.976, -0.658], + [-0.936, -0.759, -0.658, 2.353]]) + + # 2-Block [1] + L_block = array([[0,0,0,0], + [0,0,0,0], + [1,1,0,0], + [1,1,0,0]]) + W_block = array([[ 2., 0., -1., -1.], + [ 0., 2., -1., -1.], + [-1., -1., 2., 0.], + [-1., -1., 0., 2.]]) + W_block_opt = array([ + [ 2.2, -0.2, -1. , -1. ], + [-0.2, 2.2, -1. , -1. ], + [-1. , -1. , 2.2, -0.2], + [-1. , -1. , -0.2, 2.2]]) + + print('\nContraction factors of different designs with standard step size, optimal step size, and optimal W matrix') + print('Optimized step sizes and W matrix from [4] when n=4', '\n') + print('Design\t', '0.5 step size\t', 'Optimal step size\t', 'Optimal W matrix\t') + print('---------------------------------------------------------------------') + + # Malitsky-Tam [3] + tau_MT = wc_frugal_resolvent_splitting(L_MT, W_MT, lipschitz_values, mu_values, gamma=0.5, verbose=-1) + tau_MT_opt_step = wc_frugal_resolvent_splitting(L_MT, W_MT, lipschitz_values, mu_values, gamma=1.09, verbose=-1) + tau_MT_opt_W = wc_frugal_resolvent_splitting(L_MT, W_MT_opt, lipschitz_values, mu_values, gamma=1, verbose=-1) + # string format for the output of the function rounding to 3 decimal places with tab separation + print('MT \t {:.3f} \t\t {:.3f} \t\t\t {:.3f}'.format(tau_MT, tau_MT_opt_step, tau_MT_opt_W)) + + # Fully Connected + tau_f = wc_frugal_resolvent_splitting(L_full, W_full, lipschitz_values, mu_values, gamma=0.5, verbose=-1) + tau_f_opt_step = wc_frugal_resolvent_splitting(L_full, W_full, lipschitz_values, mu_values, gamma=1.09, verbose=-1) + tau_f_opt_W = wc_frugal_resolvent_splitting(L_full, W_full_opt, lipschitz_values, mu_values, gamma=1, verbose=-1) + print('Full \t {:.3f} \t\t {:.3f} \t\t\t {:.3f}'.format(tau_f, tau_f_opt_step, tau_f_opt_W)) + + # 2-Block [1] + tau_b = wc_frugal_resolvent_splitting(L_block, W_block, lipschitz_values, mu_values, gamma=0.5, verbose=-1) + tau_b_opt_step = wc_frugal_resolvent_splitting(L_block, W_block, lipschitz_values, mu_values, gamma=1.09, verbose=-1) + tau_b_opt_W = wc_frugal_resolvent_splitting(L_block, W_block_opt, lipschitz_values, mu_values, gamma=1, verbose=-1) + print('Block \t {:.3f} \t\t {:.3f} \t\t\t {:.3f}'.format(tau_b, tau_b_opt_step, tau_b_opt_W)) + return 0 + +if __name__ == "__main__": + # Douglas-Rachford [5] + pepit_tau = wc_frugal_resolvent_splitting( + L=array([[0,0],[2,0]]), + W=array([[1,-1],[-1,1]]), + lipschitz_values=[2, 1000], + mu_values=[1, 0], + alpha=1, + gamma=0.5, + wrapper="cvxpy", + verbose=1) + + # Comparison for 4 operators for Malitsky-Tam, Fully Connected, and 2-Block designs + # with and without optimized step sizes and W matrices from [1] + comparison() + + diff --git a/PEPit/examples/monotone_inclusions_variational_inequalities/reduced_frugal_resolvent_splitting.py b/PEPit/examples/monotone_inclusions_variational_inequalities/reduced_frugal_resolvent_splitting.py new file mode 100644 index 00000000..3c5ae8c3 --- /dev/null +++ b/PEPit/examples/monotone_inclusions_variational_inequalities/reduced_frugal_resolvent_splitting.py @@ -0,0 +1,253 @@ +from PEPit import PEP, null_point, Constraint +from PEPit.primitive_steps import proximal_step +from PEPit.functions import SmoothStronglyConvexFunction +from numpy import array + +def wc_reduced_frugal_resolvent_splitting(L, M, lipschitz_values, mu_values, alpha=1, gamma=0.5, wrapper="cvxpy", solver=None, verbose=1): + """ + Consider the the monotone inclusion problem + + .. math:: \\mathrm{Find}\\, x:\\, 0 \\in \\sum_{i=1}^{n} A_i(x), + + where :math:`A_i` is the subdifferential of an :math:`l_i`-Lipschitz smooth and :math:`\\mu_i`-strongly convex function for all :math:`i \\leq n`. + We denote by :math:`J_{\\alpha A_i}` the resolvent of :math:`\\alpha A_i`. + We denote the lifted vector operator :math:`\\mathbf{A}` as :math:`\\mathbf{A} = [A_1, \\dots, A_n]`, + and use lifted :math:`\\mathbf{x} = [x_1, \\dots, x_n]` and :math:`\\mathbf{w} = [w_1, \\dots, w_d]`. + We denote by :math:`L \\in \mathbb{R}^{n \\times n}` and :math:`M \\in \mathbb{R}^{n-1 \\times n}` the reduced algorithm design matrices, and by :math:`l` and :math:`\\mu` the vectors of Lipschitz and strong convexity constants of the lifted operator :math:`\\mathbf{A}`. + :math:`L` is assumed to be strictly lower diagonal. + + This code computes a worst-case guarantee for any reduced frugal resolvent splitting with design matrices :math:`L, M`. + As shown in [1] and [2], this can include the Malitsky-Tam [3], Ryu Three Operator Splitting [4], Douglas-Rachford [5], or block splitting algorithms [1]. + That is, given two lifted initial points :math:`\\mathbf{w}^{(0)}_t` and :math:`\\mathbf{w}^{(1)}_t` + this code computes the smallest possible :math:`\\tau(L, M, l, \\mu, \\alpha, \\gamma)` + (a.k.a. "contraction factor") such that the guarantee + + .. math:: \\|\\mathbf{w}^{(0)}_{t+1} - \\mathbf{w}^{(1)}_{t+1}\\|^2 \\leqslant \\tau(L, M, l, \\mu, \\alpha, \\gamma) \\|\\mathbf{w}^{(0)}_{t} - \\mathbf{w}^{(1)}_{t}\\|^2, + + is valid, where :math:`\\mathbf{w}^{(0)}_{t+1}` and :math:`\\mathbf{w}^{(1)}_{t+1}` are obtained after one iteration of the frugal resolvent splitting from respectively :math:`\\mathbf{w}^{(0)}_{t}` and :math:`\\mathbf{w}^{(1)}_{t}`. + + In short, for given values of :math:`L`, :math:`M`, :math:`l`, :math:`\\mu`, :math:`\\alpha` and :math:`\\gamma`, the contraction factor :math:`\\tau(L, M, \\mu, \\alpha, \\theta)` is computed as the worst-case value of + :math:`\\|\\mathbf{w}^{(0)}_{t+1} - \\mathbf{w}^{(1)}_{t+1}\\|^2` when :math:`\\|\\mathbf{w}^{(0)}_{t} - \\mathbf{w}^{(1)}_{t}\\|^2 \\leqslant 1`. + + **Algorithm**: One iteration of the parameterized frugal resolvent splitting is described as follows: + + .. math:: + :nowrap: + + \\begin{eqnarray} + \\mathbf{x}_{t+1} & = & J_{\\alpha \\mathbf{A}} (\\mathbf{L} \\mathbf{x}_{t+1} - \\mathbf{M}^T \\mathbf{w}_t),\\\\ + \\mathbf{w}_{t+1} & = & \\mathbf{w}_t + \\gamma \\mathbf{M} \\mathbf{x}_{t+1}. + \\end{eqnarray} + + **References**: + + `[1] R. Bassett, P. Barkley (2024). + Optimal Design of Resolvent Splitting Algorithms. arxiv:2407.16159. + `_ + + `[2] M. Tam (2023). Frugal and decentralised resolvent splittings defined by nonexpansive operators. Optimization Letters pp 1–19. `_ + + `[3] Y. Malitsky, M. Tam (2023). Resolvent splitting for sums of monotone operators + with minimal lifting. Mathematical Programming 201(1-2):231–262. `_ + + `[4] E. Ryu (2020). Uniqueness of DRS as the 2 operator resolvent-splitting and impossibility of 3 operator resolvent-splitting. Mathematical Programming 182(1- + 2):233–273. `_ + + `[5] J. Eckstein, D. Bertsekas (1992). On the Douglas—Rachford splitting method and the proximal point algorithm for maximal monotone operators. Mathematical Programming 55:293–318. `_ + + + Args: + L (ndarray): n x n numpy array of resolvent multipliers for step 1. + M (ndarray): (n-1) x n numpy array of multipliers for steps 1 and 2. + lipschitz_values (array): n Lipschitz parameters for the subdifferentials. + mu_values (array): n convexity parameters for the subdifferentials. + alpha (float): resolvent scaling parameter. + gamma (float): step size parameter. + 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 + + Example: + >>> pepit_tau = wc_reduced_frugal_resolvent_splitting( + L=array([[0,0],[2,0]]), + M=array([[1,-1]]), + lipschitz_values=[2, 1000], + mu_values=[1, 0], + alpha=1, + gamma=0.5, + wrapper="cvxpy", + 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 2 function(s) + Function 1 : Adding 2 scalar constraint(s) ... + Function 1 : 2 scalar constraint(s) added + Function 2 : Adding 2 scalar constraint(s) ... + Function 2 : 2 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.6941669886097721 + (PEPit) Postprocessing: solver's output is not entirely feasible (smallest eigenvalue of the Gram matrix is: -1.56e-09 < 0). + Small deviation from 0 may simply be due to numerical error. Big ones should be deeply investigated. + In any case, from now the provided values of parameters are based on the projection of the Gram matrix onto the cone of symmetric semi-definite matrix. + (PEPit) Primal feasibility check: + The solver found a Gram matrix that is positive semi-definite up to an error of 1.5617660559982366e-09 + All the primal scalar constraints are verified up to an error of 1.1003756295036027e-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.7629368974077815e-08 + (PEPit) Final upper bound (dual): 0.6941669881458048 and lower bound (primal example): 0.6941669886097721 + (PEPit) Duality gap: absolute: -4.639673090167662e-10 and relative: -6.683799671113239e-10 + *** Example file: worst-case performance of parameterized frugal resolvent splitting` *** + PEPit guarantee: ||w_(t+1)^0 - w_(t+1)^1||^2 <= 0.694167 ||w_(t)^0 - w_(t)^1||^2 + `` + >>> comparison() + `` + Contraction factors of different designs with standard step size vs optimal step size + Optimized step sizes from [4] when n=4 + Design 0.5 step size Optimal step size + ------------------------------------------ + MT 0.837 0.603 + Full 0.423 0.101 + Block 0.445 0.067 + `` + + """ + + # Instantiate PEP + problem = PEP() + + # Declare operators (works for LipschitzStronglyMonotoneOperator as well) + operators = [problem.declare_function(SmoothStronglyConvexFunction, L=l, mu=mu) for l, mu in zip(lipschitz_values, mu_values)] + + # Then define the starting points v0 and v1 + n = L.shape[0] + d = M.shape[0] + w0 = [problem.set_initial_point() for _ in range(d)] + w1 = [problem.set_initial_point() for _ in range(d)] + + # Set the initial constraint that is the distance between v0 and v1 + problem.set_initial_condition(sum((w0[i] - w1[i]) ** 2 for i in range(d)) <= 1) + + # Define the step for each element of the lifted vector + def resolvent(i, x, w, L, M, alpha): + Lx = sum((L[i, j]*x[j] for j in range(i)), start=null_point) + x, _, _ = proximal_step(-M.T[i,:]@w + Lx, operators[i], alpha) + return x + + x0 = [] + x1 = [] + for i in range(n): + x0.append(resolvent(i, x0, w0, L, M, alpha)) + x1.append(resolvent(i, x1, w1, L, M, alpha)) + + z0 = [] + z1 = [] + for i in range(d): + z0.append(w0[i] + gamma*M[i,:]@x0) + z1.append(w1[i] + gamma*M[i,:]@x1) + + + # Set the performance metric to the distance between z0 and z1 + problem.set_performance_metric(sum((z0[i] - z1[i]) ** 2 for i in range(d))) + + # Solve the PEP + pepit_verbose = max(verbose, 0) + pepit_tau = problem.solve(wrapper=wrapper, solver=solver, verbose=pepit_verbose) + + # Print conclusion if required + if verbose != -1: + print('*** Example file: worst-case performance of parameterized frugal resolvent splitting` ***') + print('\tPEPit guarantee:\t ||w_(t+1)^0 - w_(t+1)^1||^2 <= {:.6} ||w_(t)^0 - w_(t)^1||^2'.format(pepit_tau)) + + # Return the worst-case guarantee of the evaluated method (and the reference theoretical value) + return pepit_tau + +def comparison(): + # Comparison for 4 operators for Malitsky-Tam, Fully Connected, and 2-Block designs + # with and without optimized step sizes and W matrices + n = 4 + lipschitz_values = [2, 2, 2, 2] + mu_values = [1, 1, 1, 1] + + # Malitsky-Tam [3] + L_MT = array([[0,0,0,0], + [1,0,0,0], + [0,1,0,0], + [1,0,1,0]]) + M_MT = array([ + [-1., 1., 0., 0.], + [ 0., -1., 1., 0.], + [ 0., 0., -1., 1.]]) + + # Fully Connected + L_full = array([[0, 0, 0, 0], + [2/3,0, 0, 0], + [2/3,2/3,0, 0], + [2/3,2/3,2/3,0]]) + M_full = array([ + [-1.155, 1.155, 0. , 0. ], + [-0.667, -0.667, 1.333, 0. ], + [-0.471, -0.471, -0.471, 1.414]]) + + # 2-Block [1] + L_block = array([[0,0,0,0], + [0,0,0,0], + [1,1,0,0], + [1,1,0,0]]) + M_block = array([ + [-1. , 1. , 0. , 0. ], + [-0.707, -0.707, 1.414, 0. ], + [-0.707, -0.707, 0. , 1.414]]) + + print('\nContraction factors of different designs with standard step size vs optimal step size') + print('Optimized step sizes from [4] when n=4', '\n') + print('Design\t', '0.5 step size\t', 'Optimal step size') + print('-------------------------------------------') + + # Malitsky-Tam [3] + tau_MT = wc_reduced_frugal_resolvent_splitting(L_MT, M_MT, lipschitz_values, mu_values, gamma=0.5, verbose=-1) + tau_MT_opt_step = wc_reduced_frugal_resolvent_splitting(L_MT, M_MT, lipschitz_values, mu_values, gamma=1.405, verbose=-1) + # string format for the output of the function rounding to 3 decimal places with tab separation + print('MT \t {:.3f} \t\t {:.3f}'.format(tau_MT, tau_MT_opt_step)) + + # Fully Connected + tau_f = wc_reduced_frugal_resolvent_splitting(L_full, M_full, lipschitz_values, mu_values, gamma=0.5, verbose=-1) + tau_f_opt_step = wc_reduced_frugal_resolvent_splitting(L_full, M_full, lipschitz_values, mu_values, gamma=1.09, verbose=-1) + print('Full \t {:.3f} \t\t {:.3f}'.format(tau_f, tau_f_opt_step)) + + # 2-Block [1] + tau_b = wc_reduced_frugal_resolvent_splitting(L_block, M_block, lipschitz_values, mu_values, gamma=0.5, verbose=-1) + tau_b_opt_step = wc_reduced_frugal_resolvent_splitting(L_block, M_block, lipschitz_values, mu_values, gamma=1.12, verbose=-1) + print('Block \t {:.3f} \t\t {:.3f}'.format(tau_b, tau_b_opt_step)) + return 0 + +if __name__ == "__main__": + # Douglas-Rachford [5] + pepit_tau = wc_reduced_frugal_resolvent_splitting( + L=array([[0,0],[2,0]]), + M=array([[1,-1]]), + lipschitz_values=[2, 1000], + mu_values=[1, 0], + alpha=1, + gamma=0.5, + wrapper="cvxpy", + verbose=1) + + # Comparison for 4 operators for Malitsky-Tam, Fully Connected, and 2-Block designs + # with and without optimized step sizes from [1] + comparison() + + diff --git a/docs/source/examples/e.rst b/docs/source/examples/e.rst index 8e08170b..f3e045b1 100644 --- a/docs/source/examples/e.rst +++ b/docs/source/examples/e.rst @@ -39,3 +39,9 @@ Optimistic gradient Past extragradient ------------------ .. autofunction:: PEPit.examples.monotone_inclusions_variational_inequalities.wc_past_extragradient + + +Frugal Resolvent splittings +--------------------------- +.. autofunction:: PEPit.examples.monotone_inclusions_variational_inequalities.wc_frugal_resolvent_splitting +.. autofunction:: PEPit.examples.monotone_inclusions_variational_inequalities.wc_reduced_frugal_resolvent_splitting \ No newline at end of file diff --git a/docs/source/whatsnew/0.3.3.rst b/docs/source/whatsnew/0.3.3.rst index 17e190fb..04ad4430 100644 --- a/docs/source/whatsnew/0.3.3.rst +++ b/docs/source/whatsnew/0.3.3.rst @@ -5,3 +5,5 @@ For users: ---------- - Silver step-sizes have been added in the PEPit.examples.unconstrained_convex_minimization module. + +- 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. From 1456abc3f4283c2c7aa2edd88f88603117363cc3 Mon Sep 17 00:00:00 2001 From: Peter Barkley Date: Wed, 7 Aug 2024 13:54:41 -0700 Subject: [PATCH 2/6] Small fix to math notation --- .../frugal_resolvent_splitting.py | 2 +- .../reduced_frugal_resolvent_splitting.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/PEPit/examples/monotone_inclusions_variational_inequalities/frugal_resolvent_splitting.py b/PEPit/examples/monotone_inclusions_variational_inequalities/frugal_resolvent_splitting.py index ad50a946..1e454316 100644 --- a/PEPit/examples/monotone_inclusions_variational_inequalities/frugal_resolvent_splitting.py +++ b/PEPit/examples/monotone_inclusions_variational_inequalities/frugal_resolvent_splitting.py @@ -13,7 +13,7 @@ def wc_frugal_resolvent_splitting(L, W, lipschitz_values, mu_values, alpha=1, ga We denote by :math:`J_{\\alpha A_i}` the resolvent of :math:`\\alpha A_i`. We denote the lifted vector operator :math:`\\mathbf{A}` as :math:`\\mathbf{A} = [A_1, \\dots, A_n]`, and use lifted :math:`\\mathbf{x} = [x_1, \\dots, x_n]` and :math:`\\mathbf{v} = [v_1, \\dots, v_n]`. - We denote by :math:`L, W \\in \mathbb{R}^{n \\times n}` the algorithm design matrices, and by :math:`l` and :math:`\\mu` the vectors of Lipschitz and strong convexity constants of the lifted operator :math:`\\mathbf{A}`. + We denote by :math:`L, W \\in \\mathbb{R}^{n \\times n}` the algorithm design matrices, and by :math:`l` and :math:`\\mu` the vectors of Lipschitz and strong convexity constants of the lifted operator :math:`\\mathbf{A}`. :math:`L` is assumed to be strictly lower diagonal. This code computes a worst-case guarantee for any frugal resolvent splitting with design matrices :math:`L, W`. diff --git a/PEPit/examples/monotone_inclusions_variational_inequalities/reduced_frugal_resolvent_splitting.py b/PEPit/examples/monotone_inclusions_variational_inequalities/reduced_frugal_resolvent_splitting.py index 3c5ae8c3..86f215a0 100644 --- a/PEPit/examples/monotone_inclusions_variational_inequalities/reduced_frugal_resolvent_splitting.py +++ b/PEPit/examples/monotone_inclusions_variational_inequalities/reduced_frugal_resolvent_splitting.py @@ -13,7 +13,7 @@ def wc_reduced_frugal_resolvent_splitting(L, M, lipschitz_values, mu_values, alp We denote by :math:`J_{\\alpha A_i}` the resolvent of :math:`\\alpha A_i`. We denote the lifted vector operator :math:`\\mathbf{A}` as :math:`\\mathbf{A} = [A_1, \\dots, A_n]`, and use lifted :math:`\\mathbf{x} = [x_1, \\dots, x_n]` and :math:`\\mathbf{w} = [w_1, \\dots, w_d]`. - We denote by :math:`L \\in \mathbb{R}^{n \\times n}` and :math:`M \\in \mathbb{R}^{n-1 \\times n}` the reduced algorithm design matrices, and by :math:`l` and :math:`\\mu` the vectors of Lipschitz and strong convexity constants of the lifted operator :math:`\\mathbf{A}`. + We denote by :math:`L \\in \\mathbb{R}^{n \\times n}` and :math:`M \\in \\mathbb{R}^{n-1 \\times n}` the reduced algorithm design matrices, and by :math:`l` and :math:`\\mu` the vectors of Lipschitz and strong convexity constants of the lifted operator :math:`\\mathbf{A}`. :math:`L` is assumed to be strictly lower diagonal. This code computes a worst-case guarantee for any reduced frugal resolvent splitting with design matrices :math:`L, M`. From 6adf853e74256dd14113a9981014ab25c79d30b4 Mon Sep 17 00:00:00 2001 From: Peter Barkley Date: Thu, 19 Sep 2024 22:05:42 -0700 Subject: [PATCH 3/6] Add requested changes to frugal resolvent splitting examples --- .../frugal_resolvent_splitting.py | 153 ++++++++++-------- .../reduced_frugal_resolvent_splitting.py | 129 +++++++++------ 2 files changed, 165 insertions(+), 117 deletions(-) diff --git a/PEPit/examples/monotone_inclusions_variational_inequalities/frugal_resolvent_splitting.py b/PEPit/examples/monotone_inclusions_variational_inequalities/frugal_resolvent_splitting.py index 1e454316..a4a0c8e9 100644 --- a/PEPit/examples/monotone_inclusions_variational_inequalities/frugal_resolvent_splitting.py +++ b/PEPit/examples/monotone_inclusions_variational_inequalities/frugal_resolvent_splitting.py @@ -3,18 +3,17 @@ from PEPit.functions import SmoothStronglyConvexFunction from numpy import array -def wc_frugal_resolvent_splitting(L, W, lipschitz_values, mu_values, alpha=1, gamma=0.5, wrapper="cvxpy", solver=None, verbose=1): +def wc_frugal_resolvent_splitting(L, W, l_values, mu_values, alpha=1, gamma=0.5, wrapper="cvxpy", solver=None, verbose=1): """ - Consider the the monotone inclusion problem + Consider the the problem - .. math:: \\mathrm{Find}\\, x:\\, 0 \\in \\sum_{i=1}^{n} A_i(x), + .. math:: \\mathrm{Find}\\, x:\\, 0 = \\sum_{i=1}^{n} A_i(x), - where :math:`A_i` is the subdifferential of an :math:`l_i`-Lipschitz smooth and :math:`\\mu_i`-strongly convex function for all :math:`i \\leq n`. + where :math:`A_i` is the subgradient of an :math:`l_i`-smooth and :math:`\\mu_i`-strongly convex function for all :math:`i \\leq n`. We denote by :math:`J_{\\alpha A_i}` the resolvent of :math:`\\alpha A_i`. We denote the lifted vector operator :math:`\\mathbf{A}` as :math:`\\mathbf{A} = [A_1, \\dots, A_n]`, - and use lifted :math:`\\mathbf{x} = [x_1, \\dots, x_n]` and :math:`\\mathbf{v} = [v_1, \\dots, v_n]`. - We denote by :math:`L, W \\in \\mathbb{R}^{n \\times n}` the algorithm design matrices, and by :math:`l` and :math:`\\mu` the vectors of Lipschitz and strong convexity constants of the lifted operator :math:`\\mathbf{A}`. - :math:`L` is assumed to be strictly lower diagonal. + and use lifted :math:`\\mathbf{x}^T = [x_1, \\dots, x_n]` and :math:`\\mathbf{v}^T = [v_1, \\dots, v_n]`. + We denote by :math:`L, W \\in \\mathbb{R}^{n \\times n}` the algorithm design matrices, and by :math:`l` and :math:`\\mu` the vectors of smoothness and strong convexity constants of the functions associated with the lifted operator :math:`\\mathbf{A}`. This code computes a worst-case guarantee for any frugal resolvent splitting with design matrices :math:`L, W`. As shown in [1] and [2], this can include the Malitsky-Tam [3], Ryu Three Operator Splitting [4], Douglas-Rachford [5], or block splitting algorithms [1]. @@ -26,7 +25,7 @@ def wc_frugal_resolvent_splitting(L, W, lipschitz_values, mu_values, alpha=1, ga is valid, where :math:`\\mathbf{v}^{(0)}_{t+1}` and :math:`\\mathbf{v}^{(1)}_{t+1}` are obtained after one iteration of the frugal resolvent splitting from respectively :math:`\\mathbf{v}^{(0)}_{t}` and :math:`\\mathbf{v}^{(1)}_{t}`. - In short, for given values of :math:`L`, :math:`W`, :math:`l`, :math:`\\mu`, :math:`\\alpha` and :math:`\\gamma`, the contraction factor :math:`\\tau(L, W, \\mu, \\alpha, \\theta)` is computed as the worst-case value of + In short, for given values of :math:`L`, :math:`W`, :math:`l`, :math:`\\mu`, :math:`\\alpha` and :math:`\\gamma`, the contraction factor :math:`\\tau(L, W, l, \\mu, \\alpha, \\gamma)` is computed as the worst-case value of :math:`\\|\\mathbf{v}^{(0)}_{t+1} - \\mathbf{v}^{(1)}_{t+1}\\|^2` when :math:`\\|\\mathbf{v}^{(0)}_{t} - \\mathbf{v}^{(1)}_{t}\\|^2 \\leqslant 1`. **Algorithm**: One iteration of the parameterized frugal resolvent splitting is described as follows: @@ -39,6 +38,10 @@ def wc_frugal_resolvent_splitting(L, W, lipschitz_values, mu_values, alpha=1, ga \\mathbf{v}_{t+1} & = & \\mathbf{v}_t - \\gamma \\mathbf{W} \\mathbf{x}_{t+1}. \\end{eqnarray} + :math:`L` is assumed to be strictly lower triangular to make each resolvent :math:`J_{A_i}` of the algorithm rely only on :math:`\\mathbf{v}_i` and the results of the previous resolvents in that iteration (and not subsequent resolvents). + + :math:`W` is assumed to be positive semi-definite with a nullspace equal to the span of the ones vector, so that :math:`W 1 = 0`, ensuring that the output vector sums to 0 like the input. + **References**: `[1] R. Bassett, P. Barkley (2024). @@ -60,8 +63,8 @@ def wc_frugal_resolvent_splitting(L, W, lipschitz_values, mu_values, alpha=1, ga Args: L (ndarray): n x n numpy array of resolvent multipliers for step 1. W (ndarray): n x n numpy array of resolvent multipliers for step 2. - lipschitz_values (array): n Lipschitz parameters for the subdifferentials. - mu_values (array): n convexity parameters for the subdifferentials. + l_values (array): n smoothness parameters for the functions. + mu_values (array): n strong convexity parameters for the functions. alpha (float): resolvent scaling parameter. gamma (float): step size parameter. wrapper (str): the name of the wrapper to be used. @@ -80,16 +83,16 @@ def wc_frugal_resolvent_splitting(L, W, lipschitz_values, mu_values, alpha=1, ga >>> pepit_tau = wc_frugal_resolvent_splitting( L=array([[0,0],[2,0]]), W=array([[1,-1],[-1,1]]), - lipschitz_values=[2, 1000], + l_values=[2, 1000], mu_values=[1, 0], alpha=1, gamma=0.5, wrapper="cvxpy", verbose=1) - ``(PEPit) Setting up the problem: size of the Gram matrix: 8x8 + (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 (3 constraint(s) added) + (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 2 scalar constraint(s) ... Function 1 : 2 scalar constraint(s) added @@ -98,54 +101,68 @@ def wc_frugal_resolvent_splitting(L, W, lipschitz_values, mu_values, alpha=1, ga (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.6941881289170055 - (PEPit) Postprocessing: solver's output is not entirely feasible (smallest eigenvalue of the Gram matrix is: -3.01e-09 < 0). + (PEPit) Solver status: optimal (wrapper:cvxpy, solver: SCS); optimal value: 0.6941671049622371 + (PEPit) Postprocessing: solver's output is not entirely feasible (smallest eigenvalue of the Gram matrix is: -4.99e-09 < 0). Small deviation from 0 may simply be due to numerical error. Big ones should be deeply investigated. In any case, from now the provided values of parameters are based on the projection of the Gram matrix onto the cone of symmetric semi-definite matrix. (PEPit) Primal feasibility check: - The solver found a Gram matrix that is positive semi-definite up to an error of 3.0070570485427986e-09 - All the primal scalar constraints are verified up to an error of 7.876224117353559e-09 + The solver found a Gram matrix that is positive semi-definite up to an error of 4.987058462382981e-09 + All the primal scalar constraints are verified up to an error of 4.26979977485864e-08 (PEPit) Dual feasibility check: - The solver found a residual matrix that is positive semi-definite + The solver found a residual matrix that is positive semi-definite up to an error of 7.726891430196974e-15 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 2.165773640630705e-06 - (PEPit) Final upper bound (dual): 0.6941851987362065 and lower bound (primal example): 0.6941881289170055 - (PEPit) Duality gap: absolute: -2.9301807989989825e-06 and relative: -4.2210183046061616e-06 - *** Example file: worst-case performance of parameterized frugal resolvent splitting` *** - PEPit guarantee: ||v_(t+1)^0 - v_(t+1)^1||^2 <= 0.694185 ||v_(t)^0 - v_(t)^1||^2 - `` + (PEPit) The worst-case guarantee proof is perfectly reconstituted up to an error of 4.912845922439013e-08 + (PEPit) Final upper bound (dual): 0.6941669727111097 and lower bound (primal example): 0.6941671049622371 + (PEPit) Duality gap: absolute: -1.3225112738268763e-07 and relative: -1.9051771027075984e-07 + *** Example file: worst-case performance of parameterized frugal resolvent splitting *** + PEPit guarantee: ||v_(t+1)^0 - v_(t+1)^1||^2 <= 0.694167 ||v_(t)^0 - v_(t)^1||^2 >>> comparison() - `` - Contraction factor of different designs with standard step size, optimal step size, and optimal W matrix - Optimized step sizes and W matrix from [4] when n=4 - Design 0.5 step size Optimal step size Optimal W matrix - --------------------------------------------------------------------- - MT 0.858 0.758 0.750 - Full 0.423 0.101 0.077 - Block 0.444 0.088 0.059 - `` + --------------------------------------------------------------------------- + Contraction factors of different designs with fixed step size, optimal step size, and optimal W matrix + when n=4 and each function has l=2, mu=1 + --------------------------------------------------------------------------- + Contraction Factors + Design constant step size 0.5 optimized step size optimized W matrix + --------------------------------------------------------------------------- + MT 0.858 0.758 0.750 + Full 0.423 0.101 0.077 + Block 0.444 0.088 0.059 + --------------------------------------------------------------------------- + Optimized step sizes and W matrix found using the dual of the PEP as in [1]. + MT is the Malitsky-Tam algorithm from [3]. + Full is the fully connected algorithm in which L has no zero entries in the lower triangle. + Block is the 2-Block design from [1]. """ + # Store the number of operators + n = W.shape[0] + + # Validate input sizes are consistent + assert L.shape == W.shape == (n,n) + assert len(l_values) == len(mu_values) == n # Instantiate PEP problem = PEP() # Declare operators (works for LipschitzStronglyMonotoneOperator as well) - operators = [problem.declare_function(SmoothStronglyConvexFunction, L=l, mu=mu) for l, mu in zip(lipschitz_values, mu_values)] + operators = [problem.declare_function(SmoothStronglyConvexFunction, L=l, mu=mu) for l, mu in zip(l_values, mu_values)] - # Then define the starting points v0 and v1 - n = W.shape[0] - v0 = [problem.set_initial_point() for _ in range(n)] - v1 = [problem.set_initial_point() for _ in range(n)] + # Define the starting points v0 and v1 for the n-1 independent values + v0 = [problem.set_initial_point() for _ in range(n-1)] + v1 = [problem.set_initial_point() for _ in range(n-1)] + + # Add the n-th value to v so that it is implicitly constrained to equal 0 + v0.append(-sum(v0, start=null_point)) + v1.append(-sum(v1, start=null_point)) # Set the initial constraint that is the distance between v0 and v1 problem.set_initial_condition(sum((v0[i] - v1[i]) ** 2 for i in range(n)) <= 1) # Constraint on the lifted starting points so each sums to 0 - v0constraint = Constraint(expression=sum((v0[i] for i in range(n)), start=null_point)**2, equality_or_inequality="equality") - v1constraint = Constraint(expression=sum((v1[i] for i in range(n)), start=null_point)**2, equality_or_inequality="equality") - problem.set_initial_condition(v0constraint) - problem.set_initial_condition(v1constraint) + # v0constraint = Constraint(expression=sum((v0[i] for i in range(n)), start=null_point)**2, equality_or_inequality="equality") + # v1constraint = Constraint(expression=sum((v1[i] for i in range(n)), start=null_point)**2, equality_or_inequality="equality") + # problem.set_initial_condition(v0constraint) + # problem.set_initial_condition(v1constraint) # Define the step for each element of the lifted vector def resolvent(i, x, v, L, alpha): @@ -159,12 +176,12 @@ def resolvent(i, x, v, L, alpha): x0.append(resolvent(i, x0, v0, L, alpha)) x1.append(resolvent(i, x1, v1, L, alpha)) + # z is the updated version of v z0 = [] z1 = [] for i in range(n): z0.append(v0[i] - gamma*W[i,:]@x0) z1.append(v1[i] - gamma*W[i,:]@x1) - # Set the performance metric to the distance between z0 and z1 problem.set_performance_metric(sum((z0[i] - z1[i]) ** 2 for i in range(n))) @@ -175,7 +192,7 @@ def resolvent(i, x, v, L, alpha): # Print conclusion if required if verbose != -1: - print('*** Example file: worst-case performance of parameterized frugal resolvent splitting` ***') + print('*** Example file: worst-case performance of parameterized frugal resolvent splitting ***') print('\tPEPit guarantee:\t ||v_(t+1)^0 - v_(t+1)^1||^2 <= {:.6} ||v_(t)^0 - v_(t)^1||^2'.format(pepit_tau)) # Return the worst-case guarantee of the evaluated method (and the reference theoretical value) @@ -184,8 +201,7 @@ def resolvent(i, x, v, L, alpha): def comparison(): # Comparison for 4 operators for Malitsky-Tam, Fully Connected, and 2-Block designs # with and without optimized step sizes and W matrices - n = 4 - lipschitz_values = [2, 2, 2, 2] + l_values = [2, 2, 2, 2] mu_values = [1, 1, 1, 1] # Malitsky-Tam [3] @@ -232,38 +248,46 @@ def comparison(): [-0.2, 2.2, -1. , -1. ], [-1. , -1. , 2.2, -0.2], [-1. , -1. , -0.2, 2.2]]) - - print('\nContraction factors of different designs with standard step size, optimal step size, and optimal W matrix') - print('Optimized step sizes and W matrix from [4] when n=4', '\n') - print('Design\t', '0.5 step size\t', 'Optimal step size\t', 'Optimal W matrix\t') - print('---------------------------------------------------------------------') + + print('---------------------------------------------------------------------------') + print('Contraction factors of different designs with fixed step size, optimal step size, and optimal W matrix\n when n=4 and each function has l=2, mu=1') + print('---------------------------------------------------------------------------') + print('\t\t\tContraction Factors\t\t\t') + print('Design\t', 'constant step size 0.5\t', 'optimized step size\t', 'optimized W matrix\t') + print('---------------------------------------------------------------------------') # Malitsky-Tam [3] - tau_MT = wc_frugal_resolvent_splitting(L_MT, W_MT, lipschitz_values, mu_values, gamma=0.5, verbose=-1) - tau_MT_opt_step = wc_frugal_resolvent_splitting(L_MT, W_MT, lipschitz_values, mu_values, gamma=1.09, verbose=-1) - tau_MT_opt_W = wc_frugal_resolvent_splitting(L_MT, W_MT_opt, lipschitz_values, mu_values, gamma=1, verbose=-1) + tau_MT = wc_frugal_resolvent_splitting(L_MT, W_MT, l_values, mu_values, gamma=0.5, verbose=-1) + tau_MT_opt_step = wc_frugal_resolvent_splitting(L_MT, W_MT, l_values, mu_values, gamma=1.09, verbose=-1) + tau_MT_opt_W = wc_frugal_resolvent_splitting(L_MT, W_MT_opt, l_values, mu_values, gamma=1, verbose=-1) # string format for the output of the function rounding to 3 decimal places with tab separation - print('MT \t {:.3f} \t\t {:.3f} \t\t\t {:.3f}'.format(tau_MT, tau_MT_opt_step, tau_MT_opt_W)) + print('MT \t {:.3f} \t\t\t {:.3f} \t\t\t {:.3f}'.format(tau_MT, tau_MT_opt_step, tau_MT_opt_W)) # Fully Connected - tau_f = wc_frugal_resolvent_splitting(L_full, W_full, lipschitz_values, mu_values, gamma=0.5, verbose=-1) - tau_f_opt_step = wc_frugal_resolvent_splitting(L_full, W_full, lipschitz_values, mu_values, gamma=1.09, verbose=-1) - tau_f_opt_W = wc_frugal_resolvent_splitting(L_full, W_full_opt, lipschitz_values, mu_values, gamma=1, verbose=-1) - print('Full \t {:.3f} \t\t {:.3f} \t\t\t {:.3f}'.format(tau_f, tau_f_opt_step, tau_f_opt_W)) + tau_f = wc_frugal_resolvent_splitting(L_full, W_full, l_values, mu_values, gamma=0.5, verbose=-1) + tau_f_opt_step = wc_frugal_resolvent_splitting(L_full, W_full, l_values, mu_values, gamma=1.09, verbose=-1) + tau_f_opt_W = wc_frugal_resolvent_splitting(L_full, W_full_opt, l_values, mu_values, gamma=1, verbose=-1) + print('Full \t {:.3f} \t\t\t {:.3f} \t\t\t {:.3f}'.format(tau_f, tau_f_opt_step, tau_f_opt_W)) # 2-Block [1] - tau_b = wc_frugal_resolvent_splitting(L_block, W_block, lipschitz_values, mu_values, gamma=0.5, verbose=-1) - tau_b_opt_step = wc_frugal_resolvent_splitting(L_block, W_block, lipschitz_values, mu_values, gamma=1.09, verbose=-1) - tau_b_opt_W = wc_frugal_resolvent_splitting(L_block, W_block_opt, lipschitz_values, mu_values, gamma=1, verbose=-1) - print('Block \t {:.3f} \t\t {:.3f} \t\t\t {:.3f}'.format(tau_b, tau_b_opt_step, tau_b_opt_W)) + tau_b = wc_frugal_resolvent_splitting(L_block, W_block, l_values, mu_values, gamma=0.5, verbose=-1) + tau_b_opt_step = wc_frugal_resolvent_splitting(L_block, W_block, l_values, mu_values, gamma=1.09, verbose=-1) + tau_b_opt_W = wc_frugal_resolvent_splitting(L_block, W_block_opt, l_values, mu_values, gamma=1, verbose=-1) + print('Block \t {:.3f} \t\t\t {:.3f} \t\t\t {:.3f}'.format(tau_b, tau_b_opt_step, tau_b_opt_W)) + print('---------------------------------------------------------------------------') + print(''' Optimized step sizes and W matrix found using the dual of the PEP as in [1]. + MT is the Malitsky-Tam algorithm from [3]. + Full is the fully connected algorithm in which L has no zero entries in the lower triangle. + Block is the 2-Block design from [1].''') return 0 if __name__ == "__main__": # Douglas-Rachford [5] + print("\n1. Basic test using Douglas-Rachford matrices\n") pepit_tau = wc_frugal_resolvent_splitting( L=array([[0,0],[2,0]]), W=array([[1,-1],[-1,1]]), - lipschitz_values=[2, 1000], + l_values=[2, 1000], mu_values=[1, 0], alpha=1, gamma=0.5, @@ -272,6 +296,7 @@ def comparison(): # Comparison for 4 operators for Malitsky-Tam, Fully Connected, and 2-Block designs # with and without optimized step sizes and W matrices from [1] + print('\n2. Comparison of various algorithm designs') comparison() diff --git a/PEPit/examples/monotone_inclusions_variational_inequalities/reduced_frugal_resolvent_splitting.py b/PEPit/examples/monotone_inclusions_variational_inequalities/reduced_frugal_resolvent_splitting.py index 86f215a0..a011c1b2 100644 --- a/PEPit/examples/monotone_inclusions_variational_inequalities/reduced_frugal_resolvent_splitting.py +++ b/PEPit/examples/monotone_inclusions_variational_inequalities/reduced_frugal_resolvent_splitting.py @@ -1,20 +1,19 @@ -from PEPit import PEP, null_point, Constraint +from PEPit import PEP, null_point from PEPit.primitive_steps import proximal_step from PEPit.functions import SmoothStronglyConvexFunction from numpy import array -def wc_reduced_frugal_resolvent_splitting(L, M, lipschitz_values, mu_values, alpha=1, gamma=0.5, wrapper="cvxpy", solver=None, verbose=1): +def wc_reduced_frugal_resolvent_splitting(L, M, l_values, mu_values, alpha=1, gamma=0.5, wrapper="cvxpy", solver=None, verbose=1): """ - Consider the the monotone inclusion problem + Consider the the problem - .. math:: \\mathrm{Find}\\, x:\\, 0 \\in \\sum_{i=1}^{n} A_i(x), + .. math:: \\mathrm{Find}\\, x:\\, 0 = \\sum_{i=1}^{n} A_i(x), - where :math:`A_i` is the subdifferential of an :math:`l_i`-Lipschitz smooth and :math:`\\mu_i`-strongly convex function for all :math:`i \\leq n`. + where :math:`A_i` is the subgradient of an :math:`l_i`-smooth and :math:`\\mu_i`-strongly convex function for all :math:`i \\leq n`. We denote by :math:`J_{\\alpha A_i}` the resolvent of :math:`\\alpha A_i`. We denote the lifted vector operator :math:`\\mathbf{A}` as :math:`\\mathbf{A} = [A_1, \\dots, A_n]`, - and use lifted :math:`\\mathbf{x} = [x_1, \\dots, x_n]` and :math:`\\mathbf{w} = [w_1, \\dots, w_d]`. - We denote by :math:`L \\in \\mathbb{R}^{n \\times n}` and :math:`M \\in \\mathbb{R}^{n-1 \\times n}` the reduced algorithm design matrices, and by :math:`l` and :math:`\\mu` the vectors of Lipschitz and strong convexity constants of the lifted operator :math:`\\mathbf{A}`. - :math:`L` is assumed to be strictly lower diagonal. + and use lifted :math:`\\mathbf{x}^T = [x_1, \\dots, x_n]` and :math:`\\mathbf{w}^T = [w_1, \\dots, w_d]`. + We denote by :math:`L \\in \\mathbb{R}^{n \\times n}` and :math:`M \\in \\mathbb{R}^{n-1 \\times n}` the reduced algorithm design matrices, and by :math:`l` and :math:`\\mu` the vectors of smoothness and strong convexity constants of the functions associated with the lifted operator :math:`\\mathbf{A}`. This code computes a worst-case guarantee for any reduced frugal resolvent splitting with design matrices :math:`L, M`. As shown in [1] and [2], this can include the Malitsky-Tam [3], Ryu Three Operator Splitting [4], Douglas-Rachford [5], or block splitting algorithms [1]. @@ -24,12 +23,12 @@ def wc_reduced_frugal_resolvent_splitting(L, M, lipschitz_values, mu_values, alp .. math:: \\|\\mathbf{w}^{(0)}_{t+1} - \\mathbf{w}^{(1)}_{t+1}\\|^2 \\leqslant \\tau(L, M, l, \\mu, \\alpha, \\gamma) \\|\\mathbf{w}^{(0)}_{t} - \\mathbf{w}^{(1)}_{t}\\|^2, - is valid, where :math:`\\mathbf{w}^{(0)}_{t+1}` and :math:`\\mathbf{w}^{(1)}_{t+1}` are obtained after one iteration of the frugal resolvent splitting from respectively :math:`\\mathbf{w}^{(0)}_{t}` and :math:`\\mathbf{w}^{(1)}_{t}`. + is valid, where :math:`\\mathbf{w}^{(0)}_{t+1}` and :math:`\\mathbf{w}^{(1)}_{t+1}` are obtained after one iteration of the reduced frugal resolvent splitting from respectively :math:`\\mathbf{w}^{(0)}_{t}` and :math:`\\mathbf{w}^{(1)}_{t}`. - In short, for given values of :math:`L`, :math:`M`, :math:`l`, :math:`\\mu`, :math:`\\alpha` and :math:`\\gamma`, the contraction factor :math:`\\tau(L, M, \\mu, \\alpha, \\theta)` is computed as the worst-case value of + In short, for given values of :math:`L`, :math:`M`, :math:`l`, :math:`\\mu`, :math:`\\alpha` and :math:`\\gamma`, the contraction factor :math:`\\tau(L, M, l, \\mu, \\alpha, \\gamma)` is computed as the worst-case value of :math:`\\|\\mathbf{w}^{(0)}_{t+1} - \\mathbf{w}^{(1)}_{t+1}\\|^2` when :math:`\\|\\mathbf{w}^{(0)}_{t} - \\mathbf{w}^{(1)}_{t}\\|^2 \\leqslant 1`. - **Algorithm**: One iteration of the parameterized frugal resolvent splitting is described as follows: + **Algorithm**: One iteration of the reduced parameterized frugal resolvent splitting is described as follows: .. math:: :nowrap: @@ -39,6 +38,10 @@ def wc_reduced_frugal_resolvent_splitting(L, M, lipschitz_values, mu_values, alp \\mathbf{w}_{t+1} & = & \\mathbf{w}_t + \\gamma \\mathbf{M} \\mathbf{x}_{t+1}. \\end{eqnarray} + :math:`L` is assumed to be strictly lower triangular to make each resolvent :math:`i` of the algorithm rely only on :math:`\\mathbf{v}_i` and the results of the previous resolvents in that iteration (and not subsequent resolvents). + + :math:`M` is assumed to have nullspace equal to the span of the ones vector, so that :math:`M 1 = 0`. + **References**: `[1] R. Bassett, P. Barkley (2024). @@ -59,8 +62,8 @@ def wc_reduced_frugal_resolvent_splitting(L, M, lipschitz_values, mu_values, alp Args: L (ndarray): n x n numpy array of resolvent multipliers for step 1. M (ndarray): (n-1) x n numpy array of multipliers for steps 1 and 2. - lipschitz_values (array): n Lipschitz parameters for the subdifferentials. - mu_values (array): n convexity parameters for the subdifferentials. + l_values (array): n smoothness parameters for the functions. + mu_values (array): n strong convexity parameters for the functions. alpha (float): resolvent scaling parameter. gamma (float): step size parameter. wrapper (str): the name of the wrapper to be used. @@ -79,13 +82,13 @@ def wc_reduced_frugal_resolvent_splitting(L, M, lipschitz_values, mu_values, alp >>> pepit_tau = wc_reduced_frugal_resolvent_splitting( L=array([[0,0],[2,0]]), M=array([[1,-1]]), - lipschitz_values=[2, 1000], + l_values=[2, 1000], mu_values=[1, 0], alpha=1, gamma=0.5, wrapper="cvxpy", verbose=1) - ``(PEPit) Setting up the problem: size of the Gram matrix: 6x6 + (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) @@ -97,40 +100,52 @@ def wc_reduced_frugal_resolvent_splitting(L, M, lipschitz_values, mu_values, alp (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.6941669886097721 - (PEPit) Postprocessing: solver's output is not entirely feasible (smallest eigenvalue of the Gram matrix is: -1.56e-09 < 0). + (PEPit) Solver status: optimal (wrapper:cvxpy, solver: SCS); optimal value: 0.694166349969032 + (PEPit) Postprocessing: solver's output is not entirely feasible (smallest eigenvalue of the Gram matrix is: -5.11e-08 < 0). Small deviation from 0 may simply be due to numerical error. Big ones should be deeply investigated. In any case, from now the provided values of parameters are based on the projection of the Gram matrix onto the cone of symmetric semi-definite matrix. (PEPit) Primal feasibility check: - The solver found a Gram matrix that is positive semi-definite up to an error of 1.5617660559982366e-09 - All the primal scalar constraints are verified up to an error of 1.1003756295036027e-08 + The solver found a Gram matrix that is positive semi-definite up to an error of 5.105166188821449e-08 + All the primal scalar constraints are verified up to an error of 2.4801379455707817e-07 (PEPit) Dual feasibility check: - The solver found a residual matrix that is positive semi-definite + The solver found a residual matrix that is positive semi-definite up to an error of 1.44343037744619e-15 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.7629368974077815e-08 - (PEPit) Final upper bound (dual): 0.6941669881458048 and lower bound (primal example): 0.6941669886097721 - (PEPit) Duality gap: absolute: -4.639673090167662e-10 and relative: -6.683799671113239e-10 - *** Example file: worst-case performance of parameterized frugal resolvent splitting` *** + (PEPit) The worst-case guarantee proof is perfectly reconstituted up to an error of 1.480013558669696e-07 + (PEPit) Final upper bound (dual): 0.694166983156479 and lower bound (primal example): 0.694166349969032 + (PEPit) Duality gap: absolute: 6.331874469189813e-07 and relative: 9.121552016274325e-07 + *** Example file: worst-case performance of reduced parameterized frugal resolvent splitting *** PEPit guarantee: ||w_(t+1)^0 - w_(t+1)^1||^2 <= 0.694167 ||w_(t)^0 - w_(t)^1||^2 - `` >>> comparison() - `` - Contraction factors of different designs with standard step size vs optimal step size - Optimized step sizes from [4] when n=4 - Design 0.5 step size Optimal step size - ------------------------------------------ - MT 0.837 0.603 - Full 0.423 0.101 - Block 0.445 0.067 - `` + ---------------------------------------------------------------- + Contraction factors of different designs with constant step size 0.5 and optimized step size + when n=4 and each function has l=2, mu=1. + ---------------------------------------------------------------- + Contraction factor with Contraction factor with + Design constant step size 0.5 optimized step size + ---------------------------------------------------------------- + MT 0.837 0.603 + Full 0.423 0.101 + Block 0.445 0.067 + ---------------------------------------------------------------- + Optimized step sizes found using the dual of the PEP as in [1]. + MT is the Malitsky-Tam algorithm from [3]. + Full is the fully connected algorithm in which L has no zero entries in the lower triangle. + Block is the 2-Block design from [1]. """ + # Store the number of operators + n = L.shape[0] + d = M.shape[0] + + # Validate input sizes are consistent + assert L.shape == (n,n) + assert M.shape[1] == len(l_values) == len(mu_values) == n # Instantiate PEP problem = PEP() # Declare operators (works for LipschitzStronglyMonotoneOperator as well) - operators = [problem.declare_function(SmoothStronglyConvexFunction, L=l, mu=mu) for l, mu in zip(lipschitz_values, mu_values)] + operators = [problem.declare_function(SmoothStronglyConvexFunction, L=l, mu=mu) for l, mu in zip(l_values, mu_values)] # Then define the starting points v0 and v1 n = L.shape[0] @@ -153,12 +168,12 @@ def resolvent(i, x, w, L, M, alpha): x0.append(resolvent(i, x0, w0, L, M, alpha)) x1.append(resolvent(i, x1, w1, L, M, alpha)) + # z is the updated version of w z0 = [] z1 = [] for i in range(d): z0.append(w0[i] + gamma*M[i,:]@x0) z1.append(w1[i] + gamma*M[i,:]@x1) - # Set the performance metric to the distance between z0 and z1 problem.set_performance_metric(sum((z0[i] - z1[i]) ** 2 for i in range(d))) @@ -169,7 +184,7 @@ def resolvent(i, x, w, L, M, alpha): # Print conclusion if required if verbose != -1: - print('*** Example file: worst-case performance of parameterized frugal resolvent splitting` ***') + print('*** Example file: worst-case performance of reduced parameterized frugal resolvent splitting ***') print('\tPEPit guarantee:\t ||w_(t+1)^0 - w_(t+1)^1||^2 <= {:.6} ||w_(t)^0 - w_(t)^1||^2'.format(pepit_tau)) # Return the worst-case guarantee of the evaluated method (and the reference theoretical value) @@ -178,8 +193,7 @@ def resolvent(i, x, w, L, M, alpha): def comparison(): # Comparison for 4 operators for Malitsky-Tam, Fully Connected, and 2-Block designs # with and without optimized step sizes and W matrices - n = 4 - lipschitz_values = [2, 2, 2, 2] + l_values = [2, 2, 2, 2] mu_values = [1, 1, 1, 1] # Malitsky-Tam [3] @@ -212,34 +226,42 @@ def comparison(): [-0.707, -0.707, 1.414, 0. ], [-0.707, -0.707, 0. , 1.414]]) - print('\nContraction factors of different designs with standard step size vs optimal step size') - print('Optimized step sizes from [4] when n=4', '\n') - print('Design\t', '0.5 step size\t', 'Optimal step size') - print('-------------------------------------------') + print('----------------------------------------------------------------') + print('Contraction factors of different designs with constant step size 0.5 and optimized step size\n when n=4 and each function has l=2, mu=1.') + print('----------------------------------------------------------------') + print('\t', 'Contraction factor with\t', 'Contraction factor with') + print('Design\t', 'constant step size 0.5\t\t', 'optimized step size') + print('----------------------------------------------------------------') # Malitsky-Tam [3] - tau_MT = wc_reduced_frugal_resolvent_splitting(L_MT, M_MT, lipschitz_values, mu_values, gamma=0.5, verbose=-1) - tau_MT_opt_step = wc_reduced_frugal_resolvent_splitting(L_MT, M_MT, lipschitz_values, mu_values, gamma=1.405, verbose=-1) + tau_MT = wc_reduced_frugal_resolvent_splitting(L_MT, M_MT, l_values, mu_values, gamma=0.5, verbose=-1) + tau_MT_opt_step = wc_reduced_frugal_resolvent_splitting(L_MT, M_MT, l_values, mu_values, gamma=1.405, verbose=-1) # string format for the output of the function rounding to 3 decimal places with tab separation - print('MT \t {:.3f} \t\t {:.3f}'.format(tau_MT, tau_MT_opt_step)) + print('MT \t {:.3f} \t\t\t\t {:.3f}'.format(tau_MT, tau_MT_opt_step)) # Fully Connected - tau_f = wc_reduced_frugal_resolvent_splitting(L_full, M_full, lipschitz_values, mu_values, gamma=0.5, verbose=-1) - tau_f_opt_step = wc_reduced_frugal_resolvent_splitting(L_full, M_full, lipschitz_values, mu_values, gamma=1.09, verbose=-1) - print('Full \t {:.3f} \t\t {:.3f}'.format(tau_f, tau_f_opt_step)) + tau_f = wc_reduced_frugal_resolvent_splitting(L_full, M_full, l_values, mu_values, gamma=0.5, verbose=-1) + tau_f_opt_step = wc_reduced_frugal_resolvent_splitting(L_full, M_full, l_values, mu_values, gamma=1.09, verbose=-1) + print('Full \t {:.3f} \t\t\t\t {:.3f}'.format(tau_f, tau_f_opt_step)) # 2-Block [1] - tau_b = wc_reduced_frugal_resolvent_splitting(L_block, M_block, lipschitz_values, mu_values, gamma=0.5, verbose=-1) - tau_b_opt_step = wc_reduced_frugal_resolvent_splitting(L_block, M_block, lipschitz_values, mu_values, gamma=1.12, verbose=-1) - print('Block \t {:.3f} \t\t {:.3f}'.format(tau_b, tau_b_opt_step)) + tau_b = wc_reduced_frugal_resolvent_splitting(L_block, M_block, l_values, mu_values, gamma=0.5, verbose=-1) + tau_b_opt_step = wc_reduced_frugal_resolvent_splitting(L_block, M_block, l_values, mu_values, gamma=1.12, verbose=-1) + print('Block \t {:.3f} \t\t\t\t {:.3f}'.format(tau_b, tau_b_opt_step)) + print('----------------------------------------------------------------') + print(''' Optimized step sizes found using the dual of the PEP as in [1]. + MT is the Malitsky-Tam algorithm from [3]. + Full is the fully connected algorithm in which L has no zero entries in the lower triangle. + Block is the 2-Block design from [1].''') return 0 if __name__ == "__main__": # Douglas-Rachford [5] + print("\n1. Basic test using Douglas-Rachford matrices\n") pepit_tau = wc_reduced_frugal_resolvent_splitting( L=array([[0,0],[2,0]]), M=array([[1,-1]]), - lipschitz_values=[2, 1000], + l_values=[2, 1000], mu_values=[1, 0], alpha=1, gamma=0.5, @@ -248,6 +270,7 @@ def comparison(): # Comparison for 4 operators for Malitsky-Tam, Fully Connected, and 2-Block designs # with and without optimized step sizes from [1] + print('\n2. Comparison of various algorithm designs') comparison() From 7074b31322dc61f0139ceb5f5d08ae73a0f6a66e Mon Sep 17 00:00:00 2001 From: Peter Barkley Date: Thu, 19 Sep 2024 22:19:42 -0700 Subject: [PATCH 4/6] Remove extra import --- .../frugal_resolvent_splitting.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PEPit/examples/monotone_inclusions_variational_inequalities/frugal_resolvent_splitting.py b/PEPit/examples/monotone_inclusions_variational_inequalities/frugal_resolvent_splitting.py index a4a0c8e9..00c5e30f 100644 --- a/PEPit/examples/monotone_inclusions_variational_inequalities/frugal_resolvent_splitting.py +++ b/PEPit/examples/monotone_inclusions_variational_inequalities/frugal_resolvent_splitting.py @@ -1,4 +1,4 @@ -from PEPit import PEP, null_point, Constraint +from PEPit import PEP, null_point from PEPit.primitive_steps import proximal_step from PEPit.functions import SmoothStronglyConvexFunction from numpy import array From c87609e97e95d2a236fc9939a6a67a4d8c4be19e Mon Sep 17 00:00:00 2001 From: Peter Barkley Date: Thu, 19 Sep 2024 22:44:13 -0700 Subject: [PATCH 5/6] Remove extraneous comments --- .../frugal_resolvent_splitting.py | 6 ------ 1 file changed, 6 deletions(-) diff --git a/PEPit/examples/monotone_inclusions_variational_inequalities/frugal_resolvent_splitting.py b/PEPit/examples/monotone_inclusions_variational_inequalities/frugal_resolvent_splitting.py index 00c5e30f..f3a072fc 100644 --- a/PEPit/examples/monotone_inclusions_variational_inequalities/frugal_resolvent_splitting.py +++ b/PEPit/examples/monotone_inclusions_variational_inequalities/frugal_resolvent_splitting.py @@ -158,12 +158,6 @@ def wc_frugal_resolvent_splitting(L, W, l_values, mu_values, alpha=1, gamma=0.5, # Set the initial constraint that is the distance between v0 and v1 problem.set_initial_condition(sum((v0[i] - v1[i]) ** 2 for i in range(n)) <= 1) - # Constraint on the lifted starting points so each sums to 0 - # v0constraint = Constraint(expression=sum((v0[i] for i in range(n)), start=null_point)**2, equality_or_inequality="equality") - # v1constraint = Constraint(expression=sum((v1[i] for i in range(n)), start=null_point)**2, equality_or_inequality="equality") - # problem.set_initial_condition(v0constraint) - # problem.set_initial_condition(v1constraint) - # Define the step for each element of the lifted vector def resolvent(i, x, v, L, alpha): Lx = sum((L[i, j]*x[j] for j in range(i)), start=null_point) From 210d0384e98be5eefe2172c46e2e9765b0adfad5 Mon Sep 17 00:00:00 2001 From: Peter Barkley Date: Mon, 23 Sep 2024 10:29:44 -0700 Subject: [PATCH 6/6] Modify FRS examples to accept arbitrary problem with monotone operators. Add FRS tests using Douglas-Rachford matrices. --- .../frugal_resolvent_splitting.py | 92 +++++++------- .../reduced_frugal_resolvent_splitting.py | 115 ++++++++++-------- .../__init__.py | 3 + .../frugal_resolvent_splitting.py | 80 ++++++++++++ tests/test_uselessly_complexified_examples.py | 20 ++- 5 files changed, 213 insertions(+), 97 deletions(-) create mode 100644 tests/additional_complexified_examples_tests/frugal_resolvent_splitting.py diff --git a/PEPit/examples/monotone_inclusions_variational_inequalities/frugal_resolvent_splitting.py b/PEPit/examples/monotone_inclusions_variational_inequalities/frugal_resolvent_splitting.py index f3a072fc..f4fb0b44 100644 --- a/PEPit/examples/monotone_inclusions_variational_inequalities/frugal_resolvent_splitting.py +++ b/PEPit/examples/monotone_inclusions_variational_inequalities/frugal_resolvent_splitting.py @@ -1,31 +1,32 @@ from PEPit import PEP, null_point from PEPit.primitive_steps import proximal_step from PEPit.functions import SmoothStronglyConvexFunction +from PEPit.operators import MonotoneOperator from numpy import array -def wc_frugal_resolvent_splitting(L, W, l_values, mu_values, alpha=1, gamma=0.5, wrapper="cvxpy", solver=None, verbose=1): +def wc_frugal_resolvent_splitting(L, W, problem, alpha=1, gamma=0.5, wrapper="cvxpy", solver=None, verbose=1): """ Consider the the problem - .. math:: \\mathrm{Find}\\, x:\\, 0 = \\sum_{i=1}^{n} A_i(x), + .. math:: \\mathrm{Find}\\, x:\\, 0 \\in \\sum_{i=1}^{n} A_i(x), - where :math:`A_i` is the subgradient of an :math:`l_i`-smooth and :math:`\\mu_i`-strongly convex function for all :math:`i \\leq n`. + where :math:`A_i` is a maximal monotone operator for all :math:`i \\leq n`. We denote by :math:`J_{\\alpha A_i}` the resolvent of :math:`\\alpha A_i`. We denote the lifted vector operator :math:`\\mathbf{A}` as :math:`\\mathbf{A} = [A_1, \\dots, A_n]`, and use lifted :math:`\\mathbf{x}^T = [x_1, \\dots, x_n]` and :math:`\\mathbf{v}^T = [v_1, \\dots, v_n]`. - We denote by :math:`L, W \\in \\mathbb{R}^{n \\times n}` the algorithm design matrices, and by :math:`l` and :math:`\\mu` the vectors of smoothness and strong convexity constants of the functions associated with the lifted operator :math:`\\mathbf{A}`. + We denote by :math:`L, W \\in \\mathbb{R}^{n \\times n}` the algorithm design matrices. This code computes a worst-case guarantee for any frugal resolvent splitting with design matrices :math:`L, W`. - As shown in [1] and [2], this can include the Malitsky-Tam [3], Ryu Three Operator Splitting [4], Douglas-Rachford [5], or block splitting algorithms [1]. + As shown in [1] and [2], this includes the block splitting algorithms of [1], and (via linear transformation of the iterates), the Malitsky-Tam [3], Ryu Three Operator Splitting [4], and Douglas-Rachford [5] algorithms. That is, given two lifted initial points :math:`\\mathbf{v}^{(0)}_t` and :math:`\\mathbf{v}^{(1)}_t` (each of which sums to 0), - this code computes the smallest possible :math:`\\tau(L, W, l, \\mu, \\alpha, \\gamma)` + this code computes the smallest possible :math:`\\tau(L, W, \\alpha, \\gamma)` (a.k.a. "contraction factor") such that the guarantee - .. math:: \\|\\mathbf{v}^{(0)}_{t+1} - \\mathbf{v}^{(1)}_{t+1}\\|^2 \\leqslant \\tau(L, W, l, \\mu, \\alpha, \\gamma) \\|\\mathbf{v}^{(0)}_{t} - \\mathbf{v}^{(1)}_{t}\\|^2, + .. math:: \\|\\mathbf{v}^{(0)}_{t+1} - \\mathbf{v}^{(1)}_{t+1}\\|^2 \\leqslant \\tau(L, W, \\alpha, \\gamma) \\|\\mathbf{v}^{(0)}_{t} - \\mathbf{v}^{(1)}_{t}\\|^2, is valid, where :math:`\\mathbf{v}^{(0)}_{t+1}` and :math:`\\mathbf{v}^{(1)}_{t+1}` are obtained after one iteration of the frugal resolvent splitting from respectively :math:`\\mathbf{v}^{(0)}_{t}` and :math:`\\mathbf{v}^{(1)}_{t}`. - In short, for given values of :math:`L`, :math:`W`, :math:`l`, :math:`\\mu`, :math:`\\alpha` and :math:`\\gamma`, the contraction factor :math:`\\tau(L, W, l, \\mu, \\alpha, \\gamma)` is computed as the worst-case value of + In short, for given values of :math:`L`, :math:`W`, :math:`\\alpha` and :math:`\\gamma`, the contraction factor :math:`\\tau(L, W, \\alpha, \\gamma)` is computed as the worst-case value of :math:`\\|\\mathbf{v}^{(0)}_{t+1} - \\mathbf{v}^{(1)}_{t+1}\\|^2` when :math:`\\|\\mathbf{v}^{(0)}_{t} - \\mathbf{v}^{(1)}_{t}\\|^2 \\leqslant 1`. **Algorithm**: One iteration of the parameterized frugal resolvent splitting is described as follows: @@ -38,7 +39,7 @@ def wc_frugal_resolvent_splitting(L, W, l_values, mu_values, alpha=1, gamma=0.5, \\mathbf{v}_{t+1} & = & \\mathbf{v}_t - \\gamma \\mathbf{W} \\mathbf{x}_{t+1}. \\end{eqnarray} - :math:`L` is assumed to be strictly lower triangular to make each resolvent :math:`J_{A_i}` of the algorithm rely only on :math:`\\mathbf{v}_i` and the results of the previous resolvents in that iteration (and not subsequent resolvents). + :math:`L` is assumed to be strictly lower triangular to make each resolvent :math:`J_{A_i}` of the algorithm rely only on the value of :math:`v_i` and the results of the previous resolvents in that iteration (and not subsequent resolvents). :math:`W` is assumed to be positive semi-definite with a nullspace equal to the span of the ones vector, so that :math:`W 1 = 0`, ensuring that the output vector sums to 0 like the input. @@ -63,8 +64,7 @@ def wc_frugal_resolvent_splitting(L, W, l_values, mu_values, alpha=1, gamma=0.5, Args: L (ndarray): n x n numpy array of resolvent multipliers for step 1. W (ndarray): n x n numpy array of resolvent multipliers for step 2. - l_values (array): n smoothness parameters for the functions. - mu_values (array): n strong convexity parameters for the functions. + problem (PEP): PEP problem with exactly n maximal monotone operators. alpha (float): resolvent scaling parameter. gamma (float): step size parameter. wrapper (str): the name of the wrapper to be used. @@ -80,15 +80,13 @@ def wc_frugal_resolvent_splitting(L, W, l_values, mu_values, alpha=1, gamma=0.5, pepit_tau (float): worst-case value Example: + >>> problem = PEP() + >>> problem.declare_function(SmoothStronglyConvexFunction, L=2, mu=1) + >>> problem.declare_function(MonotoneOperator) >>> pepit_tau = wc_frugal_resolvent_splitting( L=array([[0,0],[2,0]]), W=array([[1,-1],[-1,1]]), - l_values=[2, 1000], - mu_values=[1, 0], - alpha=1, - gamma=0.5, - wrapper="cvxpy", - verbose=1) + problem=problem) (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 ... @@ -119,7 +117,7 @@ def wc_frugal_resolvent_splitting(L, W, l_values, mu_values, alpha=1, gamma=0.5, >>> comparison() --------------------------------------------------------------------------- Contraction factors of different designs with fixed step size, optimal step size, and optimal W matrix - when n=4 and each function has l=2, mu=1 + with 4 smooth strongly convex functions having l=2, mu=1 --------------------------------------------------------------------------- Contraction Factors Design constant step size 0.5 optimized step size optimized W matrix @@ -137,15 +135,12 @@ def wc_frugal_resolvent_splitting(L, W, l_values, mu_values, alpha=1, gamma=0.5, # Store the number of operators n = W.shape[0] + # Get problem operators + operators = problem.list_of_functions + # Validate input sizes are consistent assert L.shape == W.shape == (n,n) - assert len(l_values) == len(mu_values) == n - - # Instantiate PEP - problem = PEP() - - # Declare operators (works for LipschitzStronglyMonotoneOperator as well) - operators = [problem.declare_function(SmoothStronglyConvexFunction, L=l, mu=mu) for l, mu in zip(l_values, mu_values)] + assert len(operators) == n # Define the starting points v0 and v1 for the n-1 independent values v0 = [problem.set_initial_point() for _ in range(n-1)] @@ -244,30 +239,36 @@ def comparison(): [-1. , -1. , -0.2, 2.2]]) print('---------------------------------------------------------------------------') - print('Contraction factors of different designs with fixed step size, optimal step size, and optimal W matrix\n when n=4 and each function has l=2, mu=1') + print('Contraction factors of different designs with fixed step size, optimal step size, and optimal W matrix\n with 4 smooth strongly convex functions having l=2, mu=1') print('---------------------------------------------------------------------------') print('\t\t\tContraction Factors\t\t\t') print('Design\t', 'constant step size 0.5\t', 'optimized step size\t', 'optimized W matrix\t') print('---------------------------------------------------------------------------') # Malitsky-Tam [3] - tau_MT = wc_frugal_resolvent_splitting(L_MT, W_MT, l_values, mu_values, gamma=0.5, verbose=-1) - tau_MT_opt_step = wc_frugal_resolvent_splitting(L_MT, W_MT, l_values, mu_values, gamma=1.09, verbose=-1) - tau_MT_opt_W = wc_frugal_resolvent_splitting(L_MT, W_MT_opt, l_values, mu_values, gamma=1, verbose=-1) + taus = [] + for W, gamma in [(W_MT, 0.5), (W_MT, 1.09), (W_MT_opt, 1)]: + problem = PEP() + operators = [problem.declare_function(SmoothStronglyConvexFunction, L=l, mu=mu) for l, mu in zip(l_values, mu_values)] + taus.append(wc_frugal_resolvent_splitting(L_MT, W, problem, gamma=gamma, verbose=-1)) # string format for the output of the function rounding to 3 decimal places with tab separation - print('MT \t {:.3f} \t\t\t {:.3f} \t\t\t {:.3f}'.format(tau_MT, tau_MT_opt_step, tau_MT_opt_W)) + print('MT \t {:.3f} \t\t\t {:.3f} \t\t\t {:.3f}'.format(*taus)) # Fully Connected - tau_f = wc_frugal_resolvent_splitting(L_full, W_full, l_values, mu_values, gamma=0.5, verbose=-1) - tau_f_opt_step = wc_frugal_resolvent_splitting(L_full, W_full, l_values, mu_values, gamma=1.09, verbose=-1) - tau_f_opt_W = wc_frugal_resolvent_splitting(L_full, W_full_opt, l_values, mu_values, gamma=1, verbose=-1) - print('Full \t {:.3f} \t\t\t {:.3f} \t\t\t {:.3f}'.format(tau_f, tau_f_opt_step, tau_f_opt_W)) + taus = [] + for W, gamma in [(W_full, 0.5), (W_full, 1.09), (W_full_opt, 1)]: + problem = PEP() + operators = [problem.declare_function(SmoothStronglyConvexFunction, L=l, mu=mu) for l, mu in zip(l_values, mu_values)] + taus.append(wc_frugal_resolvent_splitting(L_full, W, problem, gamma=gamma, verbose=-1)) + print('Full \t {:.3f} \t\t\t {:.3f} \t\t\t {:.3f}'.format(*taus)) # 2-Block [1] - tau_b = wc_frugal_resolvent_splitting(L_block, W_block, l_values, mu_values, gamma=0.5, verbose=-1) - tau_b_opt_step = wc_frugal_resolvent_splitting(L_block, W_block, l_values, mu_values, gamma=1.09, verbose=-1) - tau_b_opt_W = wc_frugal_resolvent_splitting(L_block, W_block_opt, l_values, mu_values, gamma=1, verbose=-1) - print('Block \t {:.3f} \t\t\t {:.3f} \t\t\t {:.3f}'.format(tau_b, tau_b_opt_step, tau_b_opt_W)) + taus = [] + for W, gamma in [(W_block, 0.5), (W_block, 1.09), (W_block_opt, 1)]: + problem = PEP() + operators = [problem.declare_function(SmoothStronglyConvexFunction, L=l, mu=mu) for l, mu in zip(l_values, mu_values)] + taus.append(wc_frugal_resolvent_splitting(L_block, W, problem, gamma=gamma, verbose=-1)) + print('Block \t {:.3f} \t\t\t {:.3f} \t\t\t {:.3f}'.format(*taus)) print('---------------------------------------------------------------------------') print(''' Optimized step sizes and W matrix found using the dual of the PEP as in [1]. MT is the Malitsky-Tam algorithm from [3]. @@ -278,15 +279,18 @@ def comparison(): if __name__ == "__main__": # Douglas-Rachford [5] print("\n1. Basic test using Douglas-Rachford matrices\n") + + # Instantiate PEP + problem = PEP() + + # Declare operators + problem.declare_function(SmoothStronglyConvexFunction, L=2, mu=1) + problem.declare_function(MonotoneOperator) + pepit_tau = wc_frugal_resolvent_splitting( L=array([[0,0],[2,0]]), W=array([[1,-1],[-1,1]]), - l_values=[2, 1000], - mu_values=[1, 0], - alpha=1, - gamma=0.5, - wrapper="cvxpy", - verbose=1) + problem=problem) # Comparison for 4 operators for Malitsky-Tam, Fully Connected, and 2-Block designs # with and without optimized step sizes and W matrices from [1] diff --git a/PEPit/examples/monotone_inclusions_variational_inequalities/reduced_frugal_resolvent_splitting.py b/PEPit/examples/monotone_inclusions_variational_inequalities/reduced_frugal_resolvent_splitting.py index a011c1b2..e2a7f8a4 100644 --- a/PEPit/examples/monotone_inclusions_variational_inequalities/reduced_frugal_resolvent_splitting.py +++ b/PEPit/examples/monotone_inclusions_variational_inequalities/reduced_frugal_resolvent_splitting.py @@ -1,31 +1,32 @@ from PEPit import PEP, null_point from PEPit.primitive_steps import proximal_step from PEPit.functions import SmoothStronglyConvexFunction +from PEPit.operators import MonotoneOperator from numpy import array -def wc_reduced_frugal_resolvent_splitting(L, M, l_values, mu_values, alpha=1, gamma=0.5, wrapper="cvxpy", solver=None, verbose=1): +def wc_reduced_frugal_resolvent_splitting(L, M, problem, alpha=1, gamma=0.5, wrapper="cvxpy", solver=None, verbose=1): """ Consider the the problem - .. math:: \\mathrm{Find}\\, x:\\, 0 = \\sum_{i=1}^{n} A_i(x), + .. math:: \\mathrm{Find}\\, x:\\, 0 \\in \\sum_{i=1}^{n} A_i(x), - where :math:`A_i` is the subgradient of an :math:`l_i`-smooth and :math:`\\mu_i`-strongly convex function for all :math:`i \\leq n`. + where :math:`A_i` is a maximal monotone operator for all :math:`i \\leq n`. We denote by :math:`J_{\\alpha A_i}` the resolvent of :math:`\\alpha A_i`. We denote the lifted vector operator :math:`\\mathbf{A}` as :math:`\\mathbf{A} = [A_1, \\dots, A_n]`, and use lifted :math:`\\mathbf{x}^T = [x_1, \\dots, x_n]` and :math:`\\mathbf{w}^T = [w_1, \\dots, w_d]`. - We denote by :math:`L \\in \\mathbb{R}^{n \\times n}` and :math:`M \\in \\mathbb{R}^{n-1 \\times n}` the reduced algorithm design matrices, and by :math:`l` and :math:`\\mu` the vectors of smoothness and strong convexity constants of the functions associated with the lifted operator :math:`\\mathbf{A}`. + We denote by :math:`L \\in \\mathbb{R}^{n \\times n}` and :math:`M \\in \\mathbb{R}^{n-1 \\times n}` the reduced algorithm design matrices. This code computes a worst-case guarantee for any reduced frugal resolvent splitting with design matrices :math:`L, M`. - As shown in [1] and [2], this can include the Malitsky-Tam [3], Ryu Three Operator Splitting [4], Douglas-Rachford [5], or block splitting algorithms [1]. + As shown in [1] and [2], this can include the Malitsky-Tam [3], Ryu Three Operator Splitting [4], Douglas-Rachford [5], and the reduced version of the block splitting algorithms in [1]. That is, given two lifted initial points :math:`\\mathbf{w}^{(0)}_t` and :math:`\\mathbf{w}^{(1)}_t` - this code computes the smallest possible :math:`\\tau(L, M, l, \\mu, \\alpha, \\gamma)` + this code computes the smallest possible :math:`\\tau(L, M, \\alpha, \\gamma)` (a.k.a. "contraction factor") such that the guarantee - .. math:: \\|\\mathbf{w}^{(0)}_{t+1} - \\mathbf{w}^{(1)}_{t+1}\\|^2 \\leqslant \\tau(L, M, l, \\mu, \\alpha, \\gamma) \\|\\mathbf{w}^{(0)}_{t} - \\mathbf{w}^{(1)}_{t}\\|^2, + .. math:: \\|\\mathbf{w}^{(0)}_{t+1} - \\mathbf{w}^{(1)}_{t+1}\\|^2 \\leqslant \\tau(L, M, \\alpha, \\gamma) \\|\\mathbf{w}^{(0)}_{t} - \\mathbf{w}^{(1)}_{t}\\|^2, is valid, where :math:`\\mathbf{w}^{(0)}_{t+1}` and :math:`\\mathbf{w}^{(1)}_{t+1}` are obtained after one iteration of the reduced frugal resolvent splitting from respectively :math:`\\mathbf{w}^{(0)}_{t}` and :math:`\\mathbf{w}^{(1)}_{t}`. - In short, for given values of :math:`L`, :math:`M`, :math:`l`, :math:`\\mu`, :math:`\\alpha` and :math:`\\gamma`, the contraction factor :math:`\\tau(L, M, l, \\mu, \\alpha, \\gamma)` is computed as the worst-case value of + In short, for given values of :math:`L`, :math:`M`, :math:`\\alpha` and :math:`\\gamma`, the contraction factor :math:`\\tau(L, M, \\alpha, \\gamma)` is computed as the worst-case value of :math:`\\|\\mathbf{w}^{(0)}_{t+1} - \\mathbf{w}^{(1)}_{t+1}\\|^2` when :math:`\\|\\mathbf{w}^{(0)}_{t} - \\mathbf{w}^{(1)}_{t}\\|^2 \\leqslant 1`. **Algorithm**: One iteration of the reduced parameterized frugal resolvent splitting is described as follows: @@ -38,7 +39,7 @@ def wc_reduced_frugal_resolvent_splitting(L, M, l_values, mu_values, alpha=1, ga \\mathbf{w}_{t+1} & = & \\mathbf{w}_t + \\gamma \\mathbf{M} \\mathbf{x}_{t+1}. \\end{eqnarray} - :math:`L` is assumed to be strictly lower triangular to make each resolvent :math:`i` of the algorithm rely only on :math:`\\mathbf{v}_i` and the results of the previous resolvents in that iteration (and not subsequent resolvents). + :math:`L` is assumed to be strictly lower triangular to make each resolvent :math:`i` of the algorithm rely only on :math:`\\mathbf{w}_t` and the results of the previous resolvents in that iteration (and not subsequent resolvents). :math:`M` is assumed to have nullspace equal to the span of the ones vector, so that :math:`M 1 = 0`. @@ -62,8 +63,7 @@ def wc_reduced_frugal_resolvent_splitting(L, M, l_values, mu_values, alpha=1, ga Args: L (ndarray): n x n numpy array of resolvent multipliers for step 1. M (ndarray): (n-1) x n numpy array of multipliers for steps 1 and 2. - l_values (array): n smoothness parameters for the functions. - mu_values (array): n strong convexity parameters for the functions. + problem (PEP): PEP problem with exactly n maximal monotone operators. alpha (float): resolvent scaling parameter. gamma (float): step size parameter. wrapper (str): the name of the wrapper to be used. @@ -79,15 +79,13 @@ def wc_reduced_frugal_resolvent_splitting(L, M, l_values, mu_values, alpha=1, ga pepit_tau (float): worst-case value Example: + >>> problem = PEP() + >>> problem.declare_function(SmoothStronglyConvexFunction, L=2, mu=1) + >>> problem.declare_function(MonotoneOperator) >>> pepit_tau = wc_reduced_frugal_resolvent_splitting( L=array([[0,0],[2,0]]), M=array([[1,-1]]), - l_values=[2, 1000], - mu_values=[1, 0], - alpha=1, - gamma=0.5, - wrapper="cvxpy", - verbose=1) + problem=problem) (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 ... @@ -100,25 +98,25 @@ def wc_reduced_frugal_resolvent_splitting(L, M, l_values, mu_values, alpha=1, ga (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: SCS); optimal value: 0.694166349969032 - (PEPit) Postprocessing: solver's output is not entirely feasible (smallest eigenvalue of the Gram matrix is: -5.11e-08 < 0). + (PEPit) Solver status: optimal (wrapper:cvxpy, solver: SCS); optimal value: 0.694445372649345 + (PEPit) Postprocessing: solver's output is not entirely feasible (smallest eigenvalue of the Gram matrix is: -6.83e-08 < 0). Small deviation from 0 may simply be due to numerical error. Big ones should be deeply investigated. In any case, from now the provided values of parameters are based on the projection of the Gram matrix onto the cone of symmetric semi-definite matrix. (PEPit) Primal feasibility check: - The solver found a Gram matrix that is positive semi-definite up to an error of 5.105166188821449e-08 - All the primal scalar constraints are verified up to an error of 2.4801379455707817e-07 + The solver found a Gram matrix that is positive semi-definite up to an error of 6.830277482116174e-08 + All the primal scalar constraints are verified up to an error of 5.489757276544438e-07 (PEPit) Dual feasibility check: - The solver found a residual matrix that is positive semi-definite up to an error of 1.44343037744619e-15 - 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.480013558669696e-07 - (PEPit) Final upper bound (dual): 0.694166983156479 and lower bound (primal example): 0.694166349969032 - (PEPit) Duality gap: absolute: 6.331874469189813e-07 and relative: 9.121552016274325e-07 + The solver found a residual matrix that is positive semi-definite up to an error of 8.242470267419779e-16 + 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.710626393869319e-06 + (PEPit) Final upper bound (dual): 0.6944445325135953 and lower bound (primal example): 0.694445372649345 + (PEPit) Duality gap: absolute: -8.401357497467288e-07 and relative: -1.209793862606597e-06 *** Example file: worst-case performance of reduced parameterized frugal resolvent splitting *** - PEPit guarantee: ||w_(t+1)^0 - w_(t+1)^1||^2 <= 0.694167 ||w_(t)^0 - w_(t)^1||^2 + PEPit guarantee: ||w_(t+1)^0 - w_(t+1)^1||^2 <= 0.694445 ||w_(t)^0 - w_(t)^1||^2 >>> comparison() ---------------------------------------------------------------- Contraction factors of different designs with constant step size 0.5 and optimized step size - when n=4 and each function has l=2, mu=1. + with 4 smooth strongly convex functions having l=2, mu=1. ---------------------------------------------------------------- Contraction factor with Contraction factor with Design constant step size 0.5 optimized step size @@ -137,15 +135,12 @@ def wc_reduced_frugal_resolvent_splitting(L, M, l_values, mu_values, alpha=1, ga n = L.shape[0] d = M.shape[0] + # Get problem operators + operators = problem.list_of_functions + # Validate input sizes are consistent assert L.shape == (n,n) - assert M.shape[1] == len(l_values) == len(mu_values) == n - - # Instantiate PEP - problem = PEP() - - # Declare operators (works for LipschitzStronglyMonotoneOperator as well) - operators = [problem.declare_function(SmoothStronglyConvexFunction, L=l, mu=mu) for l, mu in zip(l_values, mu_values)] + assert M.shape[1] == len(operators) == n # Then define the starting points v0 and v1 n = L.shape[0] @@ -227,27 +222,41 @@ def comparison(): [-0.707, -0.707, 0. , 1.414]]) print('----------------------------------------------------------------') - print('Contraction factors of different designs with constant step size 0.5 and optimized step size\n when n=4 and each function has l=2, mu=1.') + print('Contraction factors of different designs with constant step size 0.5 and optimized step size\n with 4 smooth strongly convex functions having l=2, mu=1.') print('----------------------------------------------------------------') print('\t', 'Contraction factor with\t', 'Contraction factor with') print('Design\t', 'constant step size 0.5\t\t', 'optimized step size') print('----------------------------------------------------------------') # Malitsky-Tam [3] - tau_MT = wc_reduced_frugal_resolvent_splitting(L_MT, M_MT, l_values, mu_values, gamma=0.5, verbose=-1) - tau_MT_opt_step = wc_reduced_frugal_resolvent_splitting(L_MT, M_MT, l_values, mu_values, gamma=1.405, verbose=-1) - # string format for the output of the function rounding to 3 decimal places with tab separation - print('MT \t {:.3f} \t\t\t\t {:.3f}'.format(tau_MT, tau_MT_opt_step)) + taus = [] + for gamma in [0.5, 1.405]: + problem = PEP() + for l, mu in zip(l_values, mu_values): + problem.declare_function(SmoothStronglyConvexFunction, L=l, mu=mu) + taus.append(wc_reduced_frugal_resolvent_splitting(L_MT, M_MT, problem, gamma=gamma, verbose=-1)) + + print('MT \t {:.3f} \t\t\t\t {:.3f}'.format(*taus)) # Fully Connected - tau_f = wc_reduced_frugal_resolvent_splitting(L_full, M_full, l_values, mu_values, gamma=0.5, verbose=-1) - tau_f_opt_step = wc_reduced_frugal_resolvent_splitting(L_full, M_full, l_values, mu_values, gamma=1.09, verbose=-1) - print('Full \t {:.3f} \t\t\t\t {:.3f}'.format(tau_f, tau_f_opt_step)) + taus = [] + for gamma in [0.5, 1.09]: + problem = PEP() + for l, mu in zip(l_values, mu_values): + problem.declare_function(SmoothStronglyConvexFunction, L=l, mu=mu) + taus.append(wc_reduced_frugal_resolvent_splitting(L_full, M_full, problem, gamma=gamma, verbose=-1)) + + print('Full \t {:.3f} \t\t\t\t {:.3f}'.format(*taus)) # 2-Block [1] - tau_b = wc_reduced_frugal_resolvent_splitting(L_block, M_block, l_values, mu_values, gamma=0.5, verbose=-1) - tau_b_opt_step = wc_reduced_frugal_resolvent_splitting(L_block, M_block, l_values, mu_values, gamma=1.12, verbose=-1) - print('Block \t {:.3f} \t\t\t\t {:.3f}'.format(tau_b, tau_b_opt_step)) + taus = [] + for gamma in [0.5, 1.12]: + problem = PEP() + for l, mu in zip(l_values, mu_values): + problem.declare_function(SmoothStronglyConvexFunction, L=l, mu=mu) + taus.append(wc_reduced_frugal_resolvent_splitting(L_block, M_block, problem, gamma=gamma, verbose=-1)) + + print('Block \t {:.3f} \t\t\t\t {:.3f}'.format(*taus)) print('----------------------------------------------------------------') print(''' Optimized step sizes found using the dual of the PEP as in [1]. MT is the Malitsky-Tam algorithm from [3]. @@ -258,15 +267,17 @@ def comparison(): if __name__ == "__main__": # Douglas-Rachford [5] print("\n1. Basic test using Douglas-Rachford matrices\n") + + # Instantiate PEP + problem = PEP() + + # Declare operators + problem.declare_function(SmoothStronglyConvexFunction, L=2, mu=1) + problem.declare_function(MonotoneOperator) pepit_tau = wc_reduced_frugal_resolvent_splitting( L=array([[0,0],[2,0]]), M=array([[1,-1]]), - l_values=[2, 1000], - mu_values=[1, 0], - alpha=1, - gamma=0.5, - wrapper="cvxpy", - verbose=1) + problem=problem) # Comparison for 4 operators for Malitsky-Tam, Fully Connected, and 2-Block designs # with and without optimized step sizes from [1] diff --git a/tests/additional_complexified_examples_tests/__init__.py b/tests/additional_complexified_examples_tests/__init__.py index 0c390a57..5af1b507 100644 --- a/tests/additional_complexified_examples_tests/__init__.py +++ b/tests/additional_complexified_examples_tests/__init__.py @@ -11,6 +11,7 @@ from .proximal_point_LMI import wc_proximal_point_complexified3 from .randomized_coordinate_descent_smooth_convex import wc_randomized_coordinate_descent_smooth_convex_complexified from .randomized_coordinate_descent_smooth_strongly_convex import wc_randomized_coordinate_descent_smooth_strongly_convex_complexified +from .frugal_resolvent_splitting import wc_reduced_frugal_resolvent_splitting_dr, wc_frugal_resolvent_splitting_dr __all__ = ['gradient_descent_blocks', 'wc_gradient_descent_blocks', @@ -26,4 +27,6 @@ 'proximal_point_LMI', 'wc_proximal_point_complexified3', 'randomized_coordinate_descent_smooth_convex', 'wc_randomized_coordinate_descent_smooth_convex_complexified', 'randomized_coordinate_descent_smooth_strongly_convex', 'wc_randomized_coordinate_descent_smooth_strongly_convex_complexified', + 'frugal_resolvent_splitting', + 'wc_reduced_frugal_resolvent_splitting_dr', 'wc_frugal_resolvent_splitting_dr' ] diff --git a/tests/additional_complexified_examples_tests/frugal_resolvent_splitting.py b/tests/additional_complexified_examples_tests/frugal_resolvent_splitting.py new file mode 100644 index 00000000..14e16f89 --- /dev/null +++ b/tests/additional_complexified_examples_tests/frugal_resolvent_splitting.py @@ -0,0 +1,80 @@ +from numpy import array +from PEPit import PEP +from PEPit.operators import LipschitzStronglyMonotoneOperator, StronglyMonotoneOperator +from PEPit.examples.monotone_inclusions_variational_inequalities import wc_frugal_resolvent_splitting +from PEPit.examples.monotone_inclusions_variational_inequalities import wc_reduced_frugal_resolvent_splitting + + +def wc_reduced_frugal_resolvent_splitting_dr(l, mu, alpha, gamma, wrapper="cvxpy", solver=None, verbose=1): + """ + See description in `PEPit/examples/monotone_inclusions_variational_inequalities/reduced_frugal_resolvent_splitting.py`. + This example is for testing purposes. It validates that the reduced frugal resolvent splitting with the matrices below correctly gives the worst-case value for Douglas-Rachford splitting with :math:`l`-Lipschitz and maximally monotone and (maximally) :math:`\\mu`-strongly + monotone operators. + + Args: + l (float): the Lipschitz parameter. + mu (float): the strongly monotone parameter. + alpha (float): resolvent scaling parameter. + gamma (float): step size parameter. + 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 + """ + problem = PEP() + problem.declare_function(StronglyMonotoneOperator, mu=mu) + problem.declare_function(LipschitzStronglyMonotoneOperator, L=l, mu=0) + wc = wc_reduced_frugal_resolvent_splitting( + L=array([[0,0],[2,0]]), + M=array([[1,-1]]), + problem=problem, + alpha=alpha, + gamma=gamma, + wrapper=wrapper, + solver=solver, + verbose=verbose) + return wc + +def wc_frugal_resolvent_splitting_dr(l, mu, alpha, gamma, wrapper="cvxpy", solver=None, verbose=1): + """ + See description in `PEPit/examples/monotone_inclusions_variational_inequalities/frugal_resolvent_splitting.py`. + This example is for testing purposes. It validates that the frugal resolvent splitting with the matrices below correctly gives the worst-case value for Douglas-Rachford splitting with :math:`l`-Lipschitz and maximally monotone and (maximally) :math:`\\mu`-strongly + monotone operators. + + Args: + l (float): the Lipschitz parameter. + mu (float): the strongly monotone parameter. + alpha (float): resolvent scaling parameter. + gamma (float): step size parameter. + 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 + """ + problem = PEP() + problem.declare_function(StronglyMonotoneOperator, mu=mu) + problem.declare_function(LipschitzStronglyMonotoneOperator, L=l, mu=0) + wc = wc_frugal_resolvent_splitting( + L=array([[0,0],[2,0]]), + W=array([[1,-1],[-1,1]]), + problem=problem, + alpha=alpha, + gamma=gamma, + wrapper=wrapper, + solver=solver, + verbose=verbose) + return wc diff --git a/tests/test_uselessly_complexified_examples.py b/tests/test_uselessly_complexified_examples.py index be8f0d3f..9c21b1fd 100644 --- a/tests/test_uselessly_complexified_examples.py +++ b/tests/test_uselessly_complexified_examples.py @@ -22,7 +22,7 @@ from tests.additional_complexified_examples_tests import wc_randomized_coordinate_descent_smooth_convex_complexified from tests.additional_complexified_examples_tests import wc_gradient_descent_useless_blocks from tests.additional_complexified_examples_tests import wc_gradient_descent_blocks - +from tests.additional_complexified_examples_tests import wc_reduced_frugal_resolvent_splitting_dr, wc_frugal_resolvent_splitting_dr class TestExamples(unittest.TestCase): @@ -178,3 +178,21 @@ def test_gradient_descent_blocks_one_block(self): wc_modified, theory = wc_gradient_descent_blocks(L=[L], n=n, verbose=self.verbose) self.assertAlmostEqual(wc_modified, theory, delta=self.relative_precision * theory) + + def test_reduced_frugal_resolvent_splitting_dr(self): + l = 1 + mu = .1 + alpha = 1.3 + gamma = 0.9 + ref = 0.928771 # see wc_douglas_rachford_splitting + wc = wc_reduced_frugal_resolvent_splitting_dr(l, mu, alpha, gamma, verbose=self.verbose) + self.assertAlmostEqual(wc, ref, delta=self.relative_precision * ref) + + def test_frugal_resolvent_splitting_dr(self): + l = 1 + mu = .1 + alpha = 1.3 + gamma = 0.9 + ref = 0.928771 # see wc_douglas_rachford_splitting + wc = wc_frugal_resolvent_splitting_dr(l, mu, alpha, gamma, verbose=self.verbose) + self.assertAlmostEqual(wc, ref, delta=self.relative_precision * ref) \ No newline at end of file