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..f4fb0b44 --- /dev/null +++ b/PEPit/examples/monotone_inclusions_variational_inequalities/frugal_resolvent_splitting.py @@ -0,0 +1,300 @@ +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, problem, alpha=1, gamma=0.5, wrapper="cvxpy", solver=None, verbose=1): + """ + Consider the the problem + + .. math:: \\mathrm{Find}\\, x:\\, 0 \\in \\sum_{i=1}^{n} A_i(x), + + 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. + + 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 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, \\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, \\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:`\\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: + + .. 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} + + :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. + + **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. + 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. + 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: + >>> 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]]), + 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 ... + (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: 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 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 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 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 factors of different designs with fixed step size, optimal step size, and optimal W matrix + 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 + --------------------------------------------------------------------------- + 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] + + # Get problem operators + operators = problem.list_of_functions + + # Validate input sizes are consistent + assert L.shape == W.shape == (n,n) + 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)] + 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) + + # 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)) + + # 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))) + + # 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 + l_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('---------------------------------------------------------------------------') + 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] + 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(*taus)) + + # Fully Connected + 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] + 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]. + 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") + + # 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]]), + 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] + 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 new file mode 100644 index 00000000..e2a7f8a4 --- /dev/null +++ b/PEPit/examples/monotone_inclusions_variational_inequalities/reduced_frugal_resolvent_splitting.py @@ -0,0 +1,287 @@ +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, problem, alpha=1, gamma=0.5, wrapper="cvxpy", solver=None, verbose=1): + """ + Consider the the problem + + .. math:: \\mathrm{Find}\\, x:\\, 0 \\in \\sum_{i=1}^{n} A_i(x), + + 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. + + 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], 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, \\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, \\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:`\\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: + + .. 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} + + :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`. + + **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. + 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. + 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: + >>> 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]]), + 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 ... + (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: 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 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 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.694445 ||w_(t)^0 - w_(t)^1||^2 + >>> comparison() + ---------------------------------------------------------------- + Contraction factors of different designs with constant step size 0.5 and optimized step size + 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 + ---------------------------------------------------------------- + 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] + + # Get problem operators + operators = problem.list_of_functions + + # Validate input sizes are consistent + assert L.shape == (n,n) + assert M.shape[1] == len(operators) == n + + # 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)) + + # 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))) + + # 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 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) + 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 + l_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('----------------------------------------------------------------') + 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] + 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 + 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] + 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]. + 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") + + # 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]]), + problem=problem) + + # 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() + + 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 2c323991..7c5f5f7d 100644 --- a/docs/source/whatsnew/0.3.3.rst +++ b/docs/source/whatsnew/0.3.3.rst @@ -1,7 +1,9 @@ What's new in PEPit 0.3.3 ========================= -- Silver step-sizes have been added in the PEPit.examples.unconstrained_convex_minimization module. +- Silver step-sizes have been added in the `PEPit.examples.unconstrained_convex_minimization `_ module. + +- An monotone inclusion example has been added in the `PEPit.examples.monotone_inclusions_variational_inequalities `_ module. It covers parameterized frugal resolvent splittings, including the Malitsky-Tam algorithm, Douglas-Rachford, the Ryu Three Operator Splitting, and a set of novel block splitting methods. This example compares a fixed step size for this class with one optimized using the dual of the PEP. - A fix has been added in the class :class:`SmoothStronglyConvexQuadraticFunction`. Prior to that fix, using this class without stationary point in a PEP solved with the direct interface to MOSEK was problematic due to the late creation of a stationary point. After this fix, a stationary point is automatically created when instantiating this class of functions. 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