Features/fista

PEPit/examples/composite_convex_minimization/__init__.py

@@ -1,5 +1,6 @@
 from .accelerated_douglas_rachford_splitting import wc_accelerated_douglas_rachford_splitting
 from .accelerated_proximal_gradient import wc_accelerated_proximal_gradient
+from .accelerated_proximal_gradient_simplified import wc_accelerated_proximal_gradient_simplified
 from .bregman_proximal_point import wc_bregman_proximal_point
 from .douglas_rachford_splitting import wc_douglas_rachford_splitting
 from .douglas_rachford_splitting_contraction import wc_douglas_rachford_splitting_contraction
@@ -13,6 +14,7 @@
 
 __all__ = ['accelerated_douglas_rachford_splitting', 'wc_accelerated_douglas_rachford_splitting',
            'accelerated_proximal_gradient', 'wc_accelerated_proximal_gradient',
+           'accelerated_proximal_gradient_simplified', 'wc_accelerated_proximal_gradient_simplified',
            'bregman_proximal_point', 'wc_bregman_proximal_point',
            'douglas_rachford_splitting', 'wc_douglas_rachford_splitting',
            'douglas_rachford_splitting_contraction', 'wc_douglas_rachford_splitting_contraction',

PEPit/examples/composite_convex_minimization/accelerated_proximal_gradient.py

@@ -1,3 +1,5 @@
+from math import sqrt
+
 from PEPit import PEP
 from PEPit.functions import SmoothStronglyConvexFunction
 from PEPit.functions import ConvexFunction
@@ -14,7 +16,7 @@ def wc_accelerated_proximal_gradient(mu, L, n, wrapper="cvxpy", solver=None, ver
     and where :math:`h` is closed convex and proper.
 
     This code computes a worst-case guarantee for the **accelerated proximal gradient** method,
-    also known as **fast proximal gradient (FPGM)** method.
+    also known as **fast proximal gradient (FPGM)** method or FISTA [1].
     That is, it computes the smallest possible :math:`\\tau(n, L, \\mu)` such that the guarantee
 
     .. math :: F(x_n) - F(x_\\star) \\leqslant \\tau(n, L, \\mu) \\|x_0 - x_\\star\\|^2,
@@ -26,31 +28,26 @@ def wc_accelerated_proximal_gradient(mu, L, n, wrapper="cvxpy", solver=None, ver
     :math:`\\tau(n, L, \\mu)` is computed as the worst-case value of
     :math:`F(x_n) - F(x_\\star)` when :math:`\\|x_0 - x_\\star\\|^2 \\leqslant 1`.
 
-    **Algorithm**: Accelerated proximal gradient is described as follows, for :math:`t \in \\{ 0, \\dots, n-1\\}`,
+    **Algorithm**: Initialize :math:`\\lambda_1=1`, :math:`y_1=x_0`. One iteration of FISTA is described by
 
     .. math::
-        :nowrap:
 
         \\begin{eqnarray}
-            x_{t+1} & = & \\arg\\min_x \\left\\{h(x)+\\frac{L}{2}\|x-\\left(y_{t} - \\frac{1}{L} \\nabla f(y_t)\\right)\\|^2 \\right\\}, \\\\
-            y_{t+1} & = & x_{t+1} + \\frac{i}{i+3} (x_{t+1} - x_{t}),
+            \\text{Set: }\\lambda_{t+1} & = & \\frac{1 + \\sqrt{4\\lambda_t^2 + 1}}{2}\\\\
+            x_t & = & \\arg\\min_x \\left\\{h(x)+\\frac{L}{2}\|x-\\left(y_t - \\frac{1}{L} \\nabla f(y_t)\\right)\\|^2 \\right\\}\\\\
+            y_{t+1} & = & x_t + \\frac{\\lambda_t-1}{\\lambda_{t+1}} (x_t-x_{t-1}).
         \\end{eqnarray}
 
-    where :math:`y_{0} = x_0`.
-
-    **Theoretical guarantee**: A **tight** (empirical) worst-case guarantee for FPGM is obtained in
-    [1, method FPGM1 in Sec. 4.2.1, Table 1 in sec 4.2.2], for :math:`\\mu=0`:
-
-    .. math:: F(x_n) - F_\\star \\leqslant \\frac{2 L}{n^2+5n+2} \\|x_0 - x_\\star\\|^2,
+    **Theoretical guarantee**: The following worst-case guarantee can be found in e.g., [1, Theorem 4.4]:
 
-    which is attained on simple one-dimensional constrained linear optimization problems.
+    .. math:: f(x_n)-f_\\star \\leqslant \\frac{L}{2}\\frac{\\|x_0-x_\\star\\|^2}{\\lambda_n^2}.
 
     **References**:
-
-    `[1] A. Taylor, J. Hendrickx, F. Glineur (2017).
-    Exact worst-case performance of first-order methods for composite convex optimization.
-    SIAM Journal on Optimization, 27(3):1283–1313.
-    <https://arxiv.org/pdf/1512.07516.pdf>`_
+    
+    `[1] A. Beck, M. Teboulle (2009).
+    A Fast Iterative Shrinkage-Thresholding Algorithm for Linear Inverse Problems.
+    SIAM journal on imaging sciences, 2009, vol. 2, no 1, p. 183-202.
+    <https://www.ceremade.dauphine.fr/~carlier/FISTA>`_
 
 
     Args:
@@ -84,19 +81,19 @@ def wc_accelerated_proximal_gradient(mu, L, n, wrapper="cvxpy", solver=None, ver
         (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.05263158422835028
+        (PEPit) Solver status: optimal (wrapper:cvxpy, solver: MOSEK); optimal value: 0.05167329605152958
         (PEPit) Primal feasibility check:
-        		The solver found a Gram matrix that is positive semi-definite up to an error of 5.991982341524508e-09
-        		All the primal scalar constraints are verified up to an error of 1.4780313955381486e-08
+        		The solver found a Gram matrix that is positive semi-definite up to an error of 6.64684463996332e-09
+        		All the primal scalar constraints are verified up to an error of 1.6451693951591295e-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 7.783914601477293e-08
-        (PEPit) Final upper bound (dual): 0.052631589673196755 and lower bound (primal example): 0.05263158422835028 
-        (PEPit) Duality gap: absolute: 5.444846476465592e-09 and relative: 1.034520726726044e-07
+        (PEPit) The worst-case guarantee proof is perfectly reconstituted up to an error of 8.587603813802402e-08
+        (PEPit) Final upper bound (dual): 0.051673302055698395 and lower bound (primal example): 0.05167329605152958 
+        (PEPit) Duality gap: absolute: 6.004168814910393e-09 and relative: 1.1619480996379491e-07
         *** Example file: worst-case performance of the Accelerated Proximal Gradient Method in function values***
-        	PEPit guarantee:	 f(x_n)-f_* <= 0.0526316 ||x0 - xs||^2
-        	Theoretical guarantee:	 f(x_n)-f_* <= 0.0526316 ||x0 - xs||^2
+        	PEPit guarantee:	 f(x_n)-f_* <= 0.0516733 ||x0 - xs||^2
+        	Theoretical guarantee:	 f(x_n)-f_* <= 0.0661257 ||x0 - xs||^2
 
     """
 
@@ -118,26 +115,28 @@ def wc_accelerated_proximal_gradient(mu, L, n, wrapper="cvxpy", solver=None, ver
     # Set the initial constraint that is the distance between x0 and x^*
     problem.set_initial_condition((x0 - xs) ** 2 <= 1)
 
-    # Compute n steps of the accelerated proximal gradient method starting from x0
+    # Compute n steps of the accelerated proximal gradient method starting from x0    	
     x_new = x0
     y = x0
+    lam = 1
     for i in range(n):
+        lam_old = lam
+        lam = (1 + sqrt(4 * lam_old ** 2 + 1)) / 2
         x_old = x_new
         x_new, _, hx_new = proximal_step(y - 1 / L * f.gradient(y), h, 1 / L)
-        y = x_new + i / (i + 3) * (x_new - x_old)
+        y = x_new + (lam_old - 1) / lam * (x_new - x_old)
 
-    # Set the performance metric to the function value accuracy
+        # Set the performance metric to the function value accuracy
     problem.set_performance_metric((f(x_new) + hx_new) - Fs)
 
     # Solve the PEP
     pepit_verbose = max(verbose, 0)
     pepit_tau = problem.solve(wrapper=wrapper, solver=solver, verbose=pepit_verbose)
 
-    # Compute theoretical guarantee (for comparison)
-    if mu == 0:
-        theoretical_tau = 2 * L / (n ** 2 + 5 * n + 2)  # tight, see [2], Table 1 (column 1, line 1)
-    else:
-        theoretical_tau = 2 * L / (n ** 2 + 5 * n + 2)  # not tight (bound for smooth convex functions)
+    # Theoretical guarantee (for comparison)
+    theoretical_tau = L / (2 * lam_old ** 2)
+
+    if mu != 0:
         print('Warning: momentum is tuned for non-strongly convex functions.')
 
     # Print conclusion if required

PEPit/examples/composite_convex_minimization/accelerated_proximal_gradient_simplified.py

@@ -0,0 +1,155 @@
+from PEPit import PEP
+from PEPit.functions import SmoothStronglyConvexFunction
+from PEPit.functions import ConvexFunction
+from PEPit.primitive_steps import proximal_step
+
+
+def wc_accelerated_proximal_gradient_simplified(mu, L, n, wrapper="cvxpy", solver=None, verbose=1):
+    """
+    Consider the composite convex minimization problem
+
+    .. math:: F_\\star \\triangleq \\min_x \\{F(x) \equiv f(x) + h(x)\\},
+
+    where :math:`f` is :math:`L`-smooth and :math:`\\mu`-strongly convex,
+    and where :math:`h` is closed convex and proper.
+
+    This code computes a worst-case guarantee for the **accelerated proximal gradient** method,
+    also known as **fast proximal gradient (FPGM)** method.
+    That is, it computes the smallest possible :math:`\\tau(n, L, \\mu)` such that the guarantee
+
+    .. math :: F(x_n) - F(x_\\star) \\leqslant \\tau(n, L, \\mu) \\|x_0 - x_\\star\\|^2,
+
+    is valid, where :math:`x_n` is the output of the **accelerated proximal gradient** method,
+    and where :math:`x_\\star` is a minimizer of :math:`F`.
+
+    In short, for given values of :math:`n`, :math:`L` and :math:`\\mu`,
+    :math:`\\tau(n, L, \\mu)` is computed as the worst-case value of
+    :math:`F(x_n) - F(x_\\star)` when :math:`\\|x_0 - x_\\star\\|^2 \\leqslant 1`.
+
+    **Algorithm**: Accelerated proximal gradient is described as follows, for :math:`t \in \\{ 0, \\dots, n-1\\}`,
+
+    .. math::
+        :nowrap:
+
+        \\begin{eqnarray}
+            x_{t+1} & = & \\arg\\min_x \\left\\{h(x)+\\frac{L}{2}\|x-\\left(y_{t} - \\frac{1}{L} \\nabla f(y_t)\\right)\\|^2 \\right\\}, \\\\
+            y_{t+1} & = & x_{t+1} + \\frac{i}{i+3} (x_{t+1} - x_{t}),
+        \\end{eqnarray}
+
+    where :math:`y_{0} = x_0`.
+
+    **Theoretical guarantee**: A **tight** (empirical) worst-case guarantee for FPGM is obtained in
+    [1, method FPGM1 in Sec. 4.2.1, Table 1 in sec 4.2.2], for :math:`\\mu=0`:
+
+    .. math:: F(x_n) - F_\\star \\leqslant \\frac{2 L}{n^2+5n+2} \\|x_0 - x_\\star\\|^2,
+
+    which is attained on simple one-dimensional constrained linear optimization problems.
+
+    **References**:
+
+    `[1] A. Taylor, J. Hendrickx, F. Glineur (2017).
+    Exact worst-case performance of first-order methods for composite convex optimization.
+    SIAM Journal on Optimization, 27(3):1283–1313.
+    <https://arxiv.org/pdf/1512.07516.pdf>`_
+
+
+    Args:
+        L (float): the smoothness parameter.
+        mu (float): the strong convexity parameter.
+        n (int): number of iterations.
+        wrapper (str): the name of the wrapper to be used.
+        solver (str): the name of the solver the wrapper should use.
+        verbose (int): level of information details to print.
+
+                        - -1: No verbose at all.
+                        - 0: This example's output.
+                        - 1: This example's output + PEPit information.
+                        - 2: This example's output + PEPit information + solver details.
+
+    Returns:
+        pepit_tau (float): worst-case value.
+        theoretical_tau (float): theoretical value.
+
+    Example:
+        >>> pepit_tau, theoretical_tau = wc_accelerated_proximal_gradient_simplified(L=1, mu=0, n=4, wrapper="cvxpy", solver=None, verbose=1)
+        (PEPit) Setting up the problem: size of the Gram matrix: 12x12
+        (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 30 scalar constraint(s) ...
+        			Function 1 : 30 scalar constraint(s) added
+        			Function 2 : Adding 20 scalar constraint(s) ...
+        			Function 2 : 20 scalar constraint(s) added
+        (PEPit) Setting up the problem: additional constraints for 0 function(s)
+        (PEPit) Compiling SDP
+        (PEPit) Calling SDP solver
+        (PEPit) Solver status: optimal (wrapper:cvxpy, solver: MOSEK); optimal value: 0.05263158422835028
+        (PEPit) Primal feasibility check:
+        		The solver found a Gram matrix that is positive semi-definite up to an error of 5.991982341524508e-09
+        		All the primal scalar constraints are verified up to an error of 1.4780313955381486e-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 7.783914601477293e-08
+        (PEPit) Final upper bound (dual): 0.052631589673196755 and lower bound (primal example): 0.05263158422835028 
+        (PEPit) Duality gap: absolute: 5.444846476465592e-09 and relative: 1.034520726726044e-07
+        *** Example file: worst-case performance of the Accelerated Proximal Gradient Method in function values***
+        	PEPit guarantee:	 f(x_n)-f_* <= 0.0526316 ||x0 - xs||^2
+        	Theoretical guarantee:	 f(x_n)-f_* <= 0.0526316 ||x0 - xs||^2
+
+    """
+
+    # Instantiate PEP
+    problem = PEP()
+
+    # Declare a strongly convex smooth function and a convex function
+    f = problem.declare_function(SmoothStronglyConvexFunction, mu=mu, L=L)
+    h = problem.declare_function(ConvexFunction)
+    F = f + h
+
+    # Start by defining its unique optimal point xs = x_* and its function value Fs = F(x_*)
+    xs = F.stationary_point()
+    Fs = F(xs)
+
+    # Then define the starting point x0
+    x0 = problem.set_initial_point()
+
+    # Set the initial constraint that is the distance between x0 and x^*
+    problem.set_initial_condition((x0 - xs) ** 2 <= 1)
+
+    # Compute n steps of the accelerated proximal gradient method starting from x0
+    x_new = x0
+    y = x0
+    for i in range(n):
+        x_old = x_new
+        x_new, _, hx_new = proximal_step(y - 1 / L * f.gradient(y), h, 1 / L)
+        y = x_new + i / (i + 3) * (x_new - x_old)
+
+    # Set the performance metric to the function value accuracy
+    problem.set_performance_metric((f(x_new) + hx_new) - Fs)
+
+    # Solve the PEP
+    pepit_verbose = max(verbose, 0)
+    pepit_tau = problem.solve(wrapper=wrapper, solver=solver, verbose=pepit_verbose)
+
+    # Compute theoretical guarantee (for comparison)
+    theoretical_tau = 2 * L / (n ** 2 + 5 * n + 2)  # tight if mu == 0, see [1], Table 1 (column 1, line 1)
+    if mu != 0:
+        print('Warning: momentum is tuned for non-strongly convex functions.')
+
+    # Print conclusion if required
+    if verbose != -1:
+        print('*** Example file:'
+              ' worst-case performance of the Accelerated Proximal Gradient Method in function values***')
+        print('\tPEPit guarantee:\t f(x_n)-f_* <= {:.6} ||x0 - xs||^2'.format(pepit_tau))
+        print('\tTheoretical guarantee:\t f(x_n)-f_* <= {:.6} ||x0 - xs||^2'.format(theoretical_tau))
+
+    # Return the worst-case guarantee of the evaluated method ( and the reference theoretical value)
+    return pepit_tau, theoretical_tau
+
+
+if __name__ == "__main__":
+    pepit_tau, theoretical_tau = wc_accelerated_proximal_gradient_simplified(L=1, mu=0, n=4,
+                                                                             wrapper="cvxpy", solver=None,
+                                                                             verbose=1)

PEPit/examples/composite_convex_minimization/three_operator_splitting.py

@@ -97,7 +97,7 @@ def wc_three_operator_splitting(mu1, L1, L3, alpha, theta, n, wrapper="cvxpy", s
         (PEPit) Final upper bound (dual): 0.4754523347677999 and lower bound (primal example): 0.4754523346392658 
         (PEPit) Duality gap: absolute: 1.285341277856844e-10 and relative: 2.703407227628939e-10
         *** Example file: worst-case performance of the Three Operator Splitting in distance ***
-        	PEPit guarantee:	 ||w^2_n - w^1_n||^2 <= 0.475452 ||x0 - ws||^2
+        	PEPit guarantee:	 ||w^1_n - w^0_n||^2 <= 0.475452 ||w^1_0 - w^0_0||^2
 
     """
 

PEPit/examples/nonconvex_optimization/gradient_descent_quadratic_lojasiewicz_expensive.py

@@ -106,13 +106,13 @@ def wc_gradient_descent_quadratic_lojasiewicz_expensive(L, mu, gamma, n, wrapper
         		All the dual matrices to lmi are positive semi-definite
         		All the dual scalar values associated with inequality constraints are nonnegative up to an error of 5.671954340368105e-10
         (PEPit) The worst-case guarantee proof is perfectly reconstituted up to an error of 2.0306640495891322e-08
-        (PEPit) Final upper bound (dual): 0.6832669563172779 and lower bound (primal example): 0.6832669556328734
+        (PEPit) Final upper bound (dual): 0.6832669563172779 and lower bound (primal example): 0.6832669556328734 
         (PEPit) Duality gap: absolute: 6.844044220244427e-10 and relative: 1.0016647466735981e-09
         *** Example file: worst-case performance of gradient descent with fixed step-size ***
-        *** 	 (smooth problem satisfying a Lojasiewicz inequality; expert version) ***
+        *** 	 (smooth problem satisfying a Lojasiewicz inequality; expensive version) ***
         	PEPit guarantee:	 f(x_1) - f(x_*) <= 0.683267 (f(x_0)-f_*)
         	Theoretical guarantee:	 f(x_1) - f(x_*) <= 0.727273 (f(x_0)-f_*)
-
+    
     """
     # Instantiate PEP
     problem = PEP()
@@ -160,7 +160,7 @@ def wc_gradient_descent_quadratic_lojasiewicz_expensive(L, mu, gamma, n, wrapper
     # Print conclusion if required
     if verbose != -1:
         print('*** Example file: worst-case performance of gradient descent with fixed step-size ***')
-        print('*** \t (smooth problem satisfying a Lojasiewicz inequality; expert version) ***')
+        print('*** \t (smooth problem satisfying a Lojasiewicz inequality; expensive version) ***')
         print('\tPEPit guarantee:\t f(x_1) - f(x_*) <= {:.6} (f(x_0)-f_*)'.format(pepit_tau))
         print('\tTheoretical guarantee:\t f(x_1) - f(x_*) <= {:.6} (f(x_0)-f_*)'.format(theoretical_tau))