Code Listings
Contents
Code Listings#
newton_simple.m#
function [x, fval, exitflag] = newton_simple (fun, x0, options)
% Find minimum of unconstrained multivariable function.
%
% Input:
% fun - function (returning function value, gradient
% vector value and Hessian matrix value)
% x0 - initial guess
% options - structure with fields
% .MaxIterations - max. number of iterations
% .OptimalityTolerance - tolerance for the norm of gradient
% .Display - display details for 'iter'
%
% Output:
% x - solution
% fval - objective value at solution
% exitflag - 0 = no solution found
% 1 = solution found
%
narginchk(2,3);
opt = struct ();
if ((nargin == 3) && isstruct (options))
opt = options;
end
if (~isfield (opt, 'MaxIterations') || isempty (opt.MaxIterations))
opt.MaxIterations = 30;
end
if (~isfield (opt, 'Display'))
opt.Display = 'off';
end
if (~isfield (opt, 'OptimalityTolerance'))
opt.OptimalityTolerance = 1e-9;
end
if (strcmp(opt.Display, 'iter'))
fprintf('Iteration\tStep-size\tf(x)\t\t||df(x)||\n')
end
exitflag = 0;
x = x0(:);
i = 0;
while (i < opt.MaxIterations)
[fval, fgrad, fhess] = fun(x);
if (norm (fgrad, 'inf') < opt.OptimalityTolerance)
exitflag = 1; % success!
return;
end
d = fhess \ -fgrad;
if (strcmp (opt.Display, 'iter'))
fprintf('%5d\t\t%e\t%e\t%e\n', i, norm (d, 'inf'), ...
fval, norm (fgrad, 'inf'));
end
x = x + d;
i = i + 1;
end
if (i >= opt.MaxIterations)
warning ('MaxIterations = %d reached.\n', opt.MaxIterations);
end
end
nelder_mead.m#
function [x, fval, exitflag] = nelder_mead (fun, x0, options)
% Find minimum of unconstrained multivariable function.
%
% Input:
% fun - function
% x0 - initial guess
% options - structure with fields
% .MaxIterations - max. number of iterations
% .TolFun - tolerance for function values
% .Display - display details for 'iter',
% display animation for 'full'.
%
% Output:
% x - solution
% fval - objective value at solution
% exitflag - 0 = no solution found
% 1 = solution found
%
narginchk(2,3);
opt = struct ();
if ((nargin == 3) && isstruct (options))
opt = options;
end
if (~isfield (opt, 'MaxIterations') || isempty (opt.MaxIterations))
opt.MaxIterations = 30;
end
if (~isfield (opt, 'Display'))
opt.Display = 'off';
end
if (~isfield (opt, 'TolFun'))
opt.TolFun = 1e-9;
end
if (strcmp(opt.Display, 'iter') || strcmp (opt.Display, 'full'))
disp (' Iteration Func-count min f(x) Procedure');
h1 = [];
h2 = [];
h3 = [];
end
exitflag = 0;
% Parameter, here constant.
alpha = 1;
beta = 2;
gamma = 0.5;
len = 1;
% Create start simplex.
n = length (x0);
p = (len / sqrt (2)) * (sqrt (n + 1) - 1) / n;
q = [zeros(n,1), ones(n) * p];
x = repmat (x0(:), 1, n + 1) + q ...
+ [zeros(n,1), (len / sqrt (2)) * eye(n)];
% Compute all function values.
fvals = zeros (1, n + 1);
for i = 1:(n + 1)
fvals(i) = fun(x(:,i));
end
fcount = n + 1;
% Sort simplex (left is best).
[fvals, idx] = sort (fvals);
x = x(:,idx);
i = 0;
while (i < opt.MaxIterations)
% Compute centroid.
m = sum(x(:,1:n), 2) / n;
% Compute reflect.
r = m + alpha * (m - x(:,end));
fr = fun(r);
fcount = fcount + 1;
if (strcmp (opt.Display, 'full'))
delete (h1);
delete (h2);
delete (h3);
h1 = plot3 ([x(1,:), x(1,1)], [x(2,:), x(2,1)], [fvals, fvals(1)]);
h2 = plot3 (x(1,1), x(2,1), fvals(1), 'ro');
h3 = plot3 ([x(1,end), r(1)], [x(2,end), r(2)], [fvals(end), fr], 'ro');
pause (0.5);
end
if (fr <= fvals(2)) % Case 1: reflect
if (fr < fvals(1)) % Case 2: expand
e = m + beta * (m - x(:,end));
fe = fun(e);
fcount = fcount + 1;
if (fe < fr) % Expansion successful =)
r = e;
fr = fe;
end
% Replace worst x(:,end) with r (or s).
x = [r, x(:,1:n)];
fvals = [fr, fvals(1:n)];
step_method = 'expand';
else
% Replace worst x(:,end) with r.
x = [x(:,1), r, x(:,2:n)];
fvals = [fvals(1), fr, fvals(2:n)];
step_method = 'reflect';
end
else % Case 3: Contract / Shrink
if (fr > fvals(end))
step_method = 'contract inside';
c = m - gamma * (m - x(:,end));
else
step_method = 'contract outside';
c = m + gamma * (m - x(:,end));
end
fc = fun(c);
fcount = fcount + 1;
if (fc < fvals(end))
% Replace worst x(:,end) with c.
x = [x(:,1:n), c];
fvals = [fvals(1:n), fc];
else
step_method = 'shrink';
x = (x + repmat (x(:,1), 1, n + 1)) / 2;
% Compute new function values again.
for i = 1:(n + 1)
fvals(i) = fun(x(:,i));
end
fcount = fcount + n + 1;
end
% Sort again.
[fvals, idx] = sort (fvals);
x = x(:,idx);
end
fval = fvals(1);
if (strcmp (opt.Display, 'iter') || strcmp (opt.Display, 'full'))
fprintf (' %2d %2d %+8f %s\n', ...
i, fcount, fval, step_method);
end
% Check stopping criteria.
if (var (fvals) < opt.TolFun)
x = x(:,1)';
exitflag = 1; % success!
return;
end
i = i + 1;
end
x = x(:,1)';
if (i >= opt.MaxIterations)
warning ('MaxIterations = %d reached.\n', opt.MaxIterations);
end
end
lp_solver.m#
function [x, fval, exitflag, output] = lp_solver (c, A, b, x0, options)
% LP_SOLVER Solver for linear programs (LP).
%
% Input:
%
% min c'*x
% s.t. A *x = b
% x >= 0
%
% x0 - initial guess
%
% Output:
% x - solution
% fval - objective value at solution
% exitflag - 0 = no solution found
% 1 = solution found
% output - struct with fields:
% .xpath - history of x
%
narginchk (3, 5);
opt = struct ();
if ((nargin == 5) && isstruct (options))
opt = options;
end
if (~isfield (opt, 'MaxIterations') || isempty (opt.MaxIterations))
opt.MaxIterations = 10;
end
if (~isfield (opt, 'Display'))
opt.Display = 'off';
end
if (~isfield (opt, 'OptimalityTolerance'))
opt.OptimalityTolerance = 1e-8;
end
c = c(:);
b = b(:);
n = length (c);
m = length (b);
% Hessian matrix for Newton step.
H = @(x,s) [ zeros(n), -eye(n), A'; ...
diag(s), diag(x), zeros(n, m); ...
A, zeros(m,n + m) ];
% Gradient vector for Newton step.
G = @(x,s,y,mu) [ c - s + A' * y; ...
x .* s - mu * ones(n, 1); ...
A * x - b ];
% Initial values. x should be an inner LP point.
if ((nargin > 3) && ~isempty (x0))
x = x0(:);
else
x = ones (n, 1);
end
s = ones (n, 1);
y = ones (m, 1);
mu = 0.1;
% Track central path.
output.xpath = x;
if (strcmp (opt.Display, 'iter'))
fprintf('Iteration fval ||dx|| ||dy|| ||ds||\n')
end
fval = c'*x;
exitflag = 0;
% Newton iteration.
for i = 1:opt.MaxIterations
% Newton step on optimality conditions.
d = H(x,s) \ -G(x,s,y,mu);
% Stop iteration if x does not change any more.
if (norm (d(1:n), 'inf') < opt.OptimalityTolerance)
if (strcmp (opt.Display, 'iter'))
fprintf ('\nIPM converged in step %d.\n\n', i);
end
exitflag = 1;
break;
end
% Update variables.
x = x + d(1:n);
s = s + d((1:n) + n);
y = y + d((1:m) + (n * 2));
fval = c'*x;
mu = mu^2; % Shrink mu fast.
% Output step statistics.
output.xpath(:,end+1) = x;
if (strcmp (opt.Display, 'iter'))
fprintf('%5d %.2f %.2e %.2e %.2e\n', ...
i, fval, ...
norm (d(1:n), 'inf'), ...
norm (d((1:m) + n), 'inf'), ...
norm (d((1:m) + (n * 2)), 'inf'));
end
end
if (i == opt.MaxIterations)
warning ('MaxIterations = %d reached.\n', opt.MaxIterations);
end
end
um01_experiment.m#
Requires newton_simple and nelder_mead function.
function um01_experiment ()
disp ('Step 1: plotting the problem')
f = @(x,y) x.^4 + 2.*x.*y + (1 + y).^2;
[x, y] = meshgrid (linspace (-3, 3, 20));
mesh (x, y, f(x,y));
grid on;
xlabel ('x');
ylabel ('y');
zlabel ('f(x,y)');
title ('f(x,y) = x^4 + 2xy + (1 + y)^2');
hold on;
plot3 (1, -2, -2, 'ro');
view (-2, 70);
disp ('Step 2: Try fminsearch')
format long;
xpath = zeros(0, 3); % Memorize path
x0 = [1, 1];
options.Display = 'iter';
[xopt, fval, exitflag] = fminsearch (@fun, x0, options);
disp ('Solution');
disp (xopt);
disp ('Objective value at solution');
disp (fval);
fprintf ('exitflag = %d\n', exitflag);
animate_path (xpath, {'bo-'});
disp ('Step 3: Try fminunc')
xpath = zeros(0, 3); % Memorize path
%x0 = [1, 1];
options.Display = 'iter';
options.GradObj = 'on';
[xopt, fval, exitflag] = fminunc (@fun, x0, options);
disp ('Solution');
disp (xopt);
disp ('Objective value at solution');
disp (fval);
fprintf ('exitflag = %d\n', exitflag);
animate_path (xpath, {'go-'});
disp ('Step 4: Try newton_simple')
xpath = zeros(0, 3); % Memorize path
%x0 = [1, 1];
options.Display = 'iter';
[xopt, fval, exitflag] = newton_simple (@fun, x0, options);
disp ('Solution');
disp (xopt);
disp ('Objective value at solution');
disp (fval);
fprintf ('exitflag = %d\n', exitflag);
animate_path (xpath, {'ro-'});
disp ('Step 5: Try nelder_mead')
xpath = zeros(0, 3); % Memorize path
%x0 = [1, 1];
options.Display = 'iter';
options.MaxIterations = 30;
[xopt, fval, exitflag] = nelder_mead (@fun, x0, options);
disp ('Solution');
disp (xopt);
disp ('Objective value at solution');
disp (fval);
fprintf ('exitflag = %d\n', exitflag);
animate_path (xpath, {'mo-'});
% Nested function to pass to fminsearch and fminunc.
function [fx, gx, hx] = fun (x)
fx = x(1).^4 + 2.*x(1).*x(2) + (1 + x(2)).^2;
gx = [ 4.*x(1).^3 + 2.*x(2); ...
2.*x(1) + 2.*(1 + x(2)) ];
hx = [ 12.*x(1).^2, 2; ...
2, 2 ];
xpath = [xpath; x(:)', fx]; % Memorize path
end
end
function animate_path (xpath, LineSpec)
for i = 2:size (xpath, 1)
plot3 (xpath(i-1:i,1), xpath(i-1:i,2), xpath(i-1:i,3), LineSpec{:});
pause (0.1);
drawnow ();
end
end
um02_experiment.m#
Requires newton_simple and nelder_mead function.
function um02_experiment ()
disp ('Step 1: plotting the problem')
f = @(x,y) x.^3 - 6.*x.*y + 8.*y.^3;
[x, y] = meshgrid (linspace (-2, 2, 20));
mesh (x, y, f(x,y));
grid on;
xlabel ('x');
ylabel ('y');
zlabel ('f(x,y)');
title ('f(x,y) = x^3 - 6xy + 8y^3');
hold on;
plot3 (1, 0.5, -1, 'ro');
view (-60, 50);
disp ('Step 2: Try fminsearch')
format long;
xpath = zeros(0, 3); % Memorize path
x0 = [1, 1];
options.Display = 'iter';
[xopt, fval, exitflag] = fminsearch (@fun, x0, options);
disp ('Solution');
disp (xopt);
disp ('Objective value at solution');
disp (fval);
fprintf ('exitflag = %d\n', exitflag);
animate_path (xpath, {'bo-'});
disp ('Step 3: Try fminunc')
xpath = zeros(0, 3); % Memorize path
%x0 = [1, 1];
options.Display = 'iter';
options.GradObj = 'on';
[xopt, fval, exitflag] = fminunc (@fun, x0, options);
disp ('Solution');
disp (xopt);
disp ('Objective value at solution');
disp (fval);
fprintf ('exitflag = %d\n', exitflag);
animate_path (xpath, {'go-'});
disp ('Step 4: Try newton_simple')
xpath = zeros(0, 3); % Memorize path
%x0 = [1, 1];
options.Display = 'iter';
[xopt, fval, exitflag] = newton_simple (@fun, x0, options);
disp ('Solution');
disp (xopt);
disp ('Objective value at solution');
disp (fval);
fprintf ('exitflag = %d\n', exitflag);
animate_path (xpath, {'ro-'});
disp ('Step 5: Try nelder_mead')
xpath = zeros(0, 3); % Memorize path
%x0 = [1, 1];
options.Display = 'iter';
options.MaxIterations = 4;
[xopt, fval, exitflag] = nelder_mead (@fun, x0, options);
disp ('Solution');
disp (xopt);
disp ('Objective value at solution');
disp (fval);
fprintf ('exitflag = %d\n', exitflag);
animate_path (xpath, {'mo-'});
% Nested function to pass to fminsearch and fminunc.
function [fx, gx, hx] = fun (x)
fx = x(1).^3 - 6.*x(1).*x(2) + 8.*x(2).^3;
gx = [ 3.*x(1).^2 - 6.*x(2); ...
24.*x(2).^2 - 6.*x(1) ];
hx = [ 6.*x(1), -6; ...
-6, 48.*x(2) ];
xpath = [xpath; x(:)', fx]; % Memorize path
end
end
function animate_path (xpath, LineSpec)
for i = 2:size (xpath, 1)
plot3 (xpath(i-1:i,1), xpath(i-1:i,2), xpath(i-1:i,3), LineSpec{:});
pause (0.1);
drawnow ();
end
end
um03_experiment.m#
Requires newton_simple and nelder_mead function.
function um03_experiment ()
disp ('Step 1: plotting the problem')
f = @(x,y) 100 * (y - x.^2).^2 + (1 - x).^2;
[x, y] = meshgrid (linspace (-2, 2, 30), linspace (-1, 4, 30));
mesh (x, y, f(x,y));
grid on;
xlabel ('x');
ylabel ('y');
zlabel ('f(x,y)');
title ('Rosenbrock function');
hold on;
plot3 (1, 1, 0, 'ro');
view (-30, 50);
disp ('Step 2: Try fminsearch')
format long;
xpath = zeros(0, 3); % Memorize path
x0 = [-0.5, -0.5];
options.Display = 'iter';
[xopt, fval, exitflag] = fminsearch (@fun, x0, options);
disp ('Solution');
disp (xopt);
disp ('Objective value at solution');
disp (fval);
fprintf ('exitflag = %d\n', exitflag);
animate_path (xpath, {'bo-'});
disp ('Step 3: Try fminunc')
xpath = zeros(0, 3); % Memorize path
%x0 = [1, 1];
options.Display = 'iter';
options.GradObj = 'on';
[xopt, fval, exitflag] = fminunc (@fun, x0, options);
disp ('Solution');
disp (xopt);
disp ('Objective value at solution');
disp (fval);
fprintf ('exitflag = %d\n', exitflag);
animate_path (xpath, {'go-'});
disp ('Step 4: Try newton_simple')
xpath = zeros(0, 3); % Memorize path
%x0 = [1, 1];
options.Display = 'iter';
[xopt, fval, exitflag] = newton_simple (@fun, x0, options);
disp ('Solution');
disp (xopt);
disp ('Objective value at solution');
disp (fval);
fprintf ('exitflag = %d\n', exitflag);
animate_path (xpath, {'ro-'});
disp ('Step 5: Try nelder_mead')
xpath = zeros(0, 3); % Memorize path
%x0 = [1, 1];
options.Display = 'iter';
options.MaxIterations = 4;
[xopt, fval, exitflag] = nelder_mead (@fun, x0, options);
disp ('Solution');
disp (xopt);
disp ('Objective value at solution');
disp (fval);
fprintf ('exitflag = %d\n', exitflag);
animate_path (xpath, {'mo-'});
% Nested function to pass to fminsearch and fminunc.
function [fx, gx, hx] = fun (x)
fx = 100 * (x(2) - x(1).^2).^2 + (1 - x(1)).^2;
gx = [-400 * x(1) * (x(2) - x(1).^2) - 2 * (1 - x(1)); ...
200 * (x(2) - x(1).^2)];
hx = [(1200 * x(1).^2 - 400 * x(2) + 2), -400 * x(1); ...
-400 * x(1), 200];
xpath = [xpath; x(:)', fx]; % Memorize path
end
end
function animate_path (xpath, LineSpec)
for i = 2:size (xpath, 1)
plot3 (xpath(i-1:i,1), xpath(i-1:i,2), xpath(i-1:i,3), LineSpec{:});
pause (0.1);
drawnow ();
end
end