2. PEPit Examples

2.1. Frugal Resolvent Splitting 

oars.pepit.frugal_resolvent_splitting.wc_frugal_resolvent_splitting(L, W, lipschitz_values, mu_values, operator=PEPit.functions.SmoothStronglyConvexFunction, alpha=1, gamma=0.5, wrapper='cvxpy', solver=None, verbose=1)[source]

Consider the the monotone inclusion problem

\[\mathrm{Find}\, x:\, 0 \in \sum_{i=1}^{n} A_i(x),\]

where \(A_i\) is the subdifferential of an \(l_i\)-Lipschitz smooth and \(\mu_i\)-strongly convex function for all \(i \leq n\). We denote by \(J_{\alpha A_i}\) the resolvent of \(\alpha A_i\). We denote the lifted vector operator \(\mathbf{A}\) as \(\mathbf{A} = [A_1, \dots, A_n]\), and use lifted \(\mathbf{x} = [x_1, \dots, x_n]\) and \(\mathbf{v} = [v_1, \dots, v_n]\). We denote by \(L, W \in \mathbb{R}^{n \times n}\) the algorithm design matrices, and by \(l\) and \(\mu\) the vectors of Lipschitz and strong convexity constants of the lifted operator \(\mathbf{A}\). \(L\) is assumed to be strictly lower diagonal.

This code computes a worst-case guarantee for any frugal resolvent splitting with design matrices \(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 \(\mathbf{v}^{(0)}_t\) and \(\mathbf{v}^{(1)}_t\) (each of which sums to 0), this code computes the smallest possible \(\tau(L, W, l, \mu, \alpha, \gamma)\) (a.k.a. “contraction factor”) such that the guarantee

\[\|\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 \(\mathbf{v}^{(0)}_{t+1}\) and \(\mathbf{v}^{(1)}_{t+1}\) are obtained after one iteration of the frugal resolvent splitting from respectively \(\mathbf{v}^{(0)}_{t}\) and \(\mathbf{v}^{(1)}_{t}\).

In short, for given values of \(L\), \(W\), \(l\), \(\mu\), \(\alpha\) and \(\gamma\), the contraction factor \(\tau(L, W, \mu, \alpha, \theta)\) is computed as the worst-case value of \(\|\mathbf{v}^{(0)}_{t+1} - \mathbf{v}^{(1)}_{t+1}\|^2\) when \(\|\mathbf{v}^{(0)}_{t} - \mathbf{v}^{(1)}_{t}\|^2 \leqslant 1\).

Algorithm: One iteration of the parameterized frugal resolvent splitting is described as follows:

\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.

Parameters:

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
``

2.2. Reduced Frugal Resolvent Splitting 

oars.pepit.reduced_frs.wc_reduced_frugal_resolvent_splitting(L, M, lipschitz_values, mu_values, operator=PEPit.functions.SmoothStronglyConvexFunction, alpha=1, gamma=0.5, wrapper='cvxpy', solver=None, verbose=1)[source]

Consider the the monotone inclusion problem

\[\mathrm{Find}\, x:\, 0 \in \sum_{i=1}^{n} A_i(x),\]

where \(A_i\) is the subdifferential of an \(l_i\)-Lipschitz smooth and \(\mu_i\)-strongly convex function for all \(i \leq n\). We denote by \(J_{\alpha A_i}\) the resolvent of \(\alpha A_i\). We denote the lifted vector operator \(\mathbf{A}\) as \(\mathbf{A} = [A_1, \dots, A_n]\), and use lifted \(\mathbf{x} = [x_1, \dots, x_n]\) and \(\mathbf{w} = [w_1, \dots, w_d]\). We denote by \(L \in \mathbb{R}^{n \times n}\) and \(M \in \mathbb{R}^{n-1 \times n}\) the reduced algorithm design matrices, and by \(l\) and \(\mu\) the vectors of Lipschitz and strong convexity constants of the lifted operator \(\mathbf{A}\). \(L\) is assumed to be strictly lower diagonal.

This code computes a worst-case guarantee for any reduced frugal resolvent splitting with design matrices \(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 \(\mathbf{w}^{(0)}_t\) and \(\mathbf{w}^{(1)}_t\) this code computes the smallest possible \(\tau(L, M, l, \mu, \alpha, \gamma)\) (a.k.a. “contraction factor”) such that the guarantee

\[\|\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 \(\mathbf{w}^{(0)}_{t+1}\) and \(\mathbf{w}^{(1)}_{t+1}\) are obtained after one iteration of the frugal resolvent splitting from respectively \(\mathbf{w}^{(0)}_{t}\) and \(\mathbf{w}^{(1)}_{t}\).

In short, for given values of \(L\), \(M\), \(l\), \(\mu\), \(\alpha\) and \(\gamma\), the contraction factor \(\tau(L, M, \mu, \alpha, \theta)\) is computed as the worst-case value of \(\|\mathbf{w}^{(0)}_{t+1} - \mathbf{w}^{(1)}_{t+1}\|^2\) when \(\|\mathbf{w}^{(0)}_{t} - \mathbf{w}^{(1)}_{t}\|^2 \leqslant 1\).