Update to include more examples and operator classes

Author: AdrienTaylor

Examples/ConditionalGradient.m

@@ -54,5 +54,4 @@
 %
 % See Jaggi, Martin. "Revisiting Frank-Wolfe: Projection-free sparse 
 %     convex optimization." In: Proceedings of the 30th International
-%     Conference on Machine Learning (ICML-13), pp. 427–435 (2013)
-
+%     Conference on Machine Learning (ICML-13), pp. 427–435 (2013)

Examples/ConjugateGradientMethod.m

@@ -0,0 +1,50 @@
+clear all; clc;
+% In this example, we use a greedy first-order method (GFOM), or conjugate
+% gradient, for solving the L-smooth (possibly mu-strongly) 
+% convex minimization problem
+%   min_x F(x); for notational convenience we denote xs=argmin_x F(x).
+%
+% We show how to compute the worst-case value of F(xN)-F(xs) when xN is
+% obtained by doing N steps of the gradient method starting with an initial
+% iterate satisfying ||x0-xs||<=1.
+
+
+% (0) Initialize an empty PEP
+P=pep();
+
+% (1) Set up the objective function
+param.L=1;      % Smoothness parameter
+param.mu=0.0;   % Strong convexity parameter
+
+% F is the objective function
+F=P.DeclareFunction('SmoothStronglyConvex',param); 
+
+% (2) Set up the starting point and initial condition
+x0=P.StartingPoint();		 % x0 is some starting point
+[xs,fs]=F.OptimalPoint(); 		 % xs is an optimal point, and fs=F(xs)
+P.InitialCondition((x0-xs)^2<=1); % Add an initial condition ||x0-xs||^2<= 1
+
+% (3) Algorithm
+N=5;		% number of iterations
+
+x=cell(N+1,1);%we store all the x's in a cell (for convenience)
+g=cell(N+1,1);%we store all the g's in a cell (for convenience)
+x{1}=x0;
+g{1}=F.gradient(x{1});
+dirs{1}=g{1};
+for i=1:N
+    [x{i+1}, g{i+1}] = exactlinesearch_step(x{i},F,dirs);
+    dirs{2+(i-1)*2}  = x{i+1} - x{1};
+    dirs{3+(i-1)*2}  = g{i+1};
+end
+
+% (4) Set up the performance measure
+fN=F.value(x{N+1});                % g=grad F(x), f=F(x)
+P.PerformanceMetric(fN-fs); % Worst-case evaluated as F(x)-F(xs)
+
+% (5) Solve the PEP
+P.solve()
+
+% (6) Evaluate the output
+double(fN-fs)   % worst-case objective function accuracy
+% The results are the same as those for the optimized gradient method.

Examples/DouglasRachfordSplitting.m

@@ -24,11 +24,11 @@
 P=pep();
 
 % (1) Set up the objective function
-paramf1.mu=.1;	% Strong convexity parameter
-paramf1.L=1;      % Smoothness parameter
-f1=P.DeclareFunction('SmoothStronglyConvex',paramf1);
-f2=P.DeclareFunction('Convex');
-F=f1+f2; % F is the objective function
+paramf1.mu = .1;        % Strong convexity parameter
+paramf1.L  = 1;         % Smoothness parameter
+f1 = P.DeclareFunction('SmoothStronglyConvex',paramf1);
+f2 = P.DeclareFunction('Convex');
+F  = f1+f2; % F is the objective function
 
 % (2) Set up the starting point and initial condition
 w0=P.StartingPoint(); % x0 is some starting point
@@ -41,22 +41,22 @@
 % the next step evaluates the oracle at the tagged point 'opt' (xs) for
 % recovering the values of g1s and g2s; this allows to guarantee that
 % g1s+g2s=0;
-[g1s,~]=f1.oracle('opt');
-[g2s,~]=f2.oracle('opt');
-lambda=2; ws=xs+lambda*g2s;
+[g1s,~] = f1.oracle('opt');
+[g2s,~] = f2.oracle('opt');
+lambda = 2; ws = xs+lambda*g2s;
 
 % Add an initial condition ||w0-ws||^2<= 1
 P.InitialCondition((w0-ws)^2-1<=0); 
 
 % (3) Algorithm
-N=5;            % number of iterations
-gam=lambda;		% step size
+N   = 5;            % number of iterations
+gam = lambda;		% step size
 
-w=w0;
+w = w0;
 for i=1:N
-    x=proximal_step(w,f2,gam);
-    y=proximal_step(2*x-w,f1,gam);
-    w=y-x+w;
+    x = proximal_step(w,f2,gam);
+    y = proximal_step(2*x-w,f1,gam);
+    w = y-x+w;
 end
 
 % (4) Set up the performance measure

Examples/DouglasRachfordSplitting_MonotoneInclusion.m

@@ -0,0 +1,75 @@
+clear all; clc;
+% In this example, we use a Douglas-Rachford splitting (DRS) 
+% method for solving a monotone inclusion problem
+%   find x st   0 \in Ax + Bx 
+% where A is L-Lipschitz and monotone and B is (maximally) mu-strongly
+% monotone. We denote by JA and JB the respective resolvents of A and B.
+%
+% One iteration of the algorithm starting from point w is as follows:
+%       x = JB( w )
+%       y = JA( 2 * x - w )
+%       z = w - theta * ( x - y )
+% and then we choose as the next iterate the value of z.
+%
+% Given two initial points w0 and w1, we show how to compute the worst-case
+% contraction factor ||z0 - z1||/||w0 - w1|| obtained after doing one
+% iteration of DRS from respectively w0 and w1.
+% Note that we allow the user to choose a stepsize alpha in the resolvent.
+%
+% This setting is studied in
+% (**) Ernest K. Ryu, Adrien B. Taylor, C. Bergeling, and P. Giselsson.
+%      "Operator Splitting Performance Estimation: Tight contraction
+%       factors and optimal parameter selection." arXiv:1812.00146, 2018.
+%
+% (0) Initialize an empty PEP
+P=pep();
+
+% (1) Set up the class of monotone inclusions
+paramA.L  =  1; paramA.mu = 0; % A is 1-Lipschitz and 0-strongly monotone
+paramB.mu = .1;                % B is .1-strongly monotone
+
+A = P.DeclareFunction('LipschitzStronglyMonotone',paramA);
+B = P.DeclareFunction('StronglyMonotone',paramB);
+
+% (2) Set up the starting points
+w0=P.StartingPoint();
+w1=P.StartingPoint();
+P.InitialCondition((w0-w1)^2<=1);  % Normalize the initial distance ||w0-ws||^2 <= 1
+
+% (3) Algorithm
+alpha = 1.3;		% step size (in the resolvents)
+theta = .9;         % overrelaxation
+
+x0 = proximal_step(w0,B,alpha);
+y0 = proximal_step(2*x0-w0,A,alpha);
+z0 = w0-theta*(x0-y0);
+
+x1 = proximal_step(w1,B,alpha);
+y1 = proximal_step(2*x1-w1,A,alpha);
+z1 = w1-theta*(x1-y1);
+
+% (4) Set up the performance measure: ||z0-z1||^2
+P.PerformanceMetric((z0-z1)^2);
+
+% (5) Solve the PEP
+P.solve()
+
+% (6) Evaluate the output
+double((z0-z1)^2)   % worst-case contraction factor
+
+% Results to be compared with WC below (from (**)):
+L    = alpha*paramA.L; mu = alpha*paramB.mu;
+C    = sqrt(((2*(theta-1)*mu+theta-2)^2+L^2*(theta-2*(mu+1))^2)/(L^2+1));
+if theta*(theta+C)/(mu+1)^2/C * (C+mu*((2*(theta-1)*mu+theta-2)-L^2*...
+        (theta-2*(mu+1)))/(L^2+1)) >=0
+    WC   = ((theta+C)/2/(mu+1))^2;
+elseif L<=1 && mu >= (L^2+1)/(L-1)^2 && theta<=-(2*(mu+1)*(L+1)*(mu+...
+        (mu-1)*L^2-2*mu*L-1))/(mu+L*(L^2+L+1)+2*mu^2*(L-1)+mu*L*(1-(L-3)*L)+1)
+    WC   = (1-theta*(L+mu)/(L+1)/(mu+1))^2;
+else
+    WC   = (2-theta)/4/mu/(L^2+1) * ...
+        (theta*(1-2*mu+L^2)-2*mu*(L^2-1))*...
+        (theta*(1+2*mu+L^2)-2*(mu+1)*(L^2+1))/...
+        (theta*(1+2*mu-L^2)-2*(mu+1)*(1-L^2));
+end
+WC

Examples/GradientExactLineSearch.m

@@ -36,11 +36,11 @@
 P.InitialCondition(f0-fs<=1); % Add an initial condition f0-fs<= 1
 
 % (3) Algorithm
-N=2;
-x=x0;
-for i=1:N
-    [g,~]=F.oracle(x);
-    x=exactlinesearch_step(x,F,g);
+N = 2;
+x = x0;
+for i = 1:N
+    g = F.gradient(x);
+    x = exactlinesearch_step(x,F,g);
 end
 
 % (4) Set up the performance measure

Examples/HalpernIteration.m

@@ -0,0 +1,40 @@
+clear all; clc;
+% In this example, we use the Halpern iteration for finding a fixed point
+% to the non-expansive operator A :
+%   find x such that  x = Ax
+%
+% (**) Lieder, Felix. "On the Convergence Rate of the Halpern-Iteration." 
+%      (2017)
+
+
+% (0) Initialize an empty PEP
+P=pep();
+
+% (1) Set up the objective function
+paramA.L=1;      % A is 1-Lipschitz (non-expansive)
+A = P.DeclareFunction('Lipschitz',paramA);
+
+% (2) Set up the starting point and initial condition
+x0 = P.StartingPoint();		 % x0 is some starting point
+xs = fixedpoint(A);
+P.InitialCondition((x0-xs)^2<=1); % Add an initial condition ||x0-xs||^2<= 1
+
+% (3) Algorithm
+N = 10;
+lambda = @(k)(1/(k+2));
+x=x0;
+for i=1:N
+    x = lambda(i-1) * x0 + (1-lambda(i-1)) * A.evaluate(x);
+end
+xN  = x;
+AxN = A.evaluate(xN);
+% (4) Set up the performance measure
+P.PerformanceMetric((xN-AxN)^2); % Worst-case squared residual
+
+% (5) Solve the PEP
+P.solve()
+
+% (6) Evaluate the output
+double((xN-AxN)^2)   % worst-case squared residual
+
+% The result should be (2/(N+1))^2 as in (**)

Examples/KrasnoselskiiMann.m

@@ -0,0 +1,41 @@
+clear all; clc;
+% In this example, we use Krasnoselskii-Mann iterations for finding a 
+% fixed point to the non-expansive operator A :
+%   find x such that  x = Ax
+%
+% This scheme was first studied using PEPs in:
+% (**) Felix Lieder. "Projection Based Methods for Conic Linear Programming 
+%       Optimal First Order Complexities and Norm Constrained Quasi Newton 
+%       Methods."  PhD thesis (2018)
+
+
+% (0) Initialize an empty PEP
+P=pep();
+
+% (1) Set up the objective function
+paramA.L=1;      % A is 1-Lipschitz (non-expansive)
+A = P.DeclareFunction('Lipschitz',paramA);
+
+% (2) Set up the starting point and initial condition
+x0 = P.StartingPoint();		 % x0 is some starting point
+xs = fixedpoint(A);
+P.InitialCondition((x0-xs)^2<=1); % Add an initial condition ||x0-xs||^2<= 1
+
+% (3) Algorithm
+N = 10;
+lambda = @(k)(1/(k+2));
+x=x0;
+for i=1:N
+    x = lambda(i-1) * x + (1-lambda(i-1)) * A.evaluate(x);
+end
+xN  = x;
+AxN = A.evaluate(xN);
+% (4) Set up the performance measure
+P.PerformanceMetric((xN-AxN)^2); % Worst-case squared residual
+
+% (5) Solve the PEP
+P.solve()
+
+% (6) Evaluate the output
+double((xN-AxN)^2)   % worst-case squared residual
+

Examples/OptimizedGradientMethod.m

@@ -4,17 +4,17 @@
 %   min_x F(x); for notational convenience we denote xs=argmin_x F(x).
 %
 % We show how to compute the worst-case value of F(xN)-F(xs) when xN is
-% obtained by doing N steps of the gradient method starting with an initial
-% iterate satisfying ||x0-xs||<=1.
+% obtained by doing N steps of OGM starting with an initial iterate
+% satisfying ||x0-xs||<=1.
 %
 % Note that OGM is developped in the following two works:
-%(1)Drori, Yoel, and Marc Teboulle. 
-%   "Performance of first-order methods for smooth convex minimization:
-%   a novel approach." Mathematical Programming 145.1-2 (2014): 451-482.
+%(1) Drori, Yoel, and Marc Teboulle. 
+%    "Performance of first-order methods for smooth convex minimization:
+%    a novel approach." Mathematical Programming 145.1-2 (2014): 451-482.
 %
-%(2)Kim, Donghwan, and Jeffrey A. Fessler. 
-%   "Optimized first-order methods for smooth convex minimization."
-%   Mathematical programming 159.1-2 (2016): 81-107.
+%(2) Kim, Donghwan, and Jeffrey A. Fessler. 
+%    "Optimized first-order methods for smooth convex minimization."
+%    Mathematical programming 159.1-2 (2016): 81-107.
 
 
 % (0) Initialize an empty PEP

Examples/OptimizedGradientMethodForGradient.m

@@ -0,0 +1,70 @@
+clear all; clc;
+% In this example, we use the optimized gradient method for gradient norm
+% (OGM-G) for solving the L-smooth convex minimization problem
+%   min_x F(x); for notational convenience we denote xs=argmin_x F(x).
+%
+% We show how to compute the worst-case value of ||F'(xN)||^2 when xN is
+% obtained by doing N steps of OGM-G starting with an initial
+% iterate satisfying F(x0)-F(x*)<=1.
+%
+% Note that OGMG is developped in the following work:
+% (1) Kim, D., & Fessler, J. A. (2018). "Optimizing the Efficiency of 
+%     First-order Methods for Decreasing the Gradient of Smooth 
+%     Convex Functions." preprint arXiv:1803.06600.
+%
+
+
+% (0) Initialize an empty PEP
+P=pep();
+L = 1;
+% (1) Set up the objective function
+param.L=L;      % Smoothness parameter
+
+F=P.DeclareFunction('SmoothStronglyConvex',param); % F is the objective function
+
+% (2) Set up the starting point and initial condition
+x0=P.StartingPoint();           % x0 is some starting point
+[xs,fs]=F.OptimalPoint(); 		% xs is an optimal point, and fs=F(xs)
+[g0,f0]=F.oracle(x0);
+P.InitialCondition(f0-fs<=1);   % Add an initial condition F(x0)-F(x*)<=1.
+
+% (3) Algorithm
+gam=1/param.L;		% step size
+N=5;		% number of iterations
+
+x=cell(N+1,1);%we store all the x's in a cell (for convenience)
+x{1}=x0;
+g{1}=g0;
+y{1}=x0;
+theta(1)=1;
+for i=1:N
+    if i<N
+        theta(i+1)=(1+sqrt(4*theta(i)^2+1))/2;
+    else
+        theta(i+1)=(1+sqrt(4*theta(i)^2+1))/2;
+    end
+end
+th_ti = @(i)theta(N+1-i);
+
+for i = 1:N
+    y{i+1} = x{i} - 1/L * g{i};
+    cc       = (2*th_ti(i)-1)/(2*th_ti(i-1)-1);
+    x{i+1} = y{i+1} + (th_ti(i-1)-1)/th_ti(i-1)*cc*(y{i+1}-y{i}) + cc * (y{i+1}-x{i});
+    g{i+1}   = F.gradient(x{i+1});
+end
+
+% (4) Set up the performance measure
+% gN=F.gradient(x{N+1});                % g=grad F(x), f=F(x)
+obj = (g{N+1})^2;
+P.PerformanceMetric(obj); % Worst-case evaluated as F(x)-F(xs)
+
+% (5) Solve the PEP
+P.solve();
+
+% (6) Evaluate the output
+double(obj)   % worst-case objective function accuracy
+
+% The result should be 2/theta(end)^2
+% see: Kim, D., & Fessler, J. A. (2018). "Optimizing the Efficiency of 
+%      First-order Methods for Decreasing the Gradient of Smooth 
+%      Convex Functions." preprint arXiv:1803.06600.