Skip to content

Instantly share code, notes, and snippets.

@agirault
Created June 1, 2014 22:06
Show Gist options
  • Save agirault/9e84b7e1c233951a3dbd to your computer and use it in GitHub Desktop.
Save agirault/9e84b7e1c233951a3dbd to your computer and use it in GitHub Desktop.
function [histout,costdata, x] = linesearch(x0,f,g,tolPos,tolGrad,maxit)
% LINESEARCH finds the minimum of an objective function using
% - Steepest descent Algorithm
% - Armijo as a condition for sufficient decrease (could try Wolf?)
% - Step length selection
%
% INPUT:
% - x0 = initial iterate position
% - f = objective function,
% - g = grad f (COLUMN vector)
% - tolPos = termination criterion norm(dx) < tolPos
% (optional) default = 1.d-6
% - tolGrad = termination criterion norm(grad) < tolGrad
% (optional) default = 1.d-6
% - maxit = maximum iterations
% (optional) default = 1000
%
% OUTPUT:
% - x = solution (minimum found)
% - histout = iteration history. Each row of histout is :
% [norm(grad), f, number of step length reductions, iteration count]
% - costdata = [num f, num grad]
%
% Requires: stepLengthSelection.m
%% INITIALISATION
% linesearch parameters
if nargin < 6
maxit = 1000;
end
if nargin < 5
tolGrad = 1.d-6;
end
if nargin < 4
tolPos = 1.d-6;
end
c1 = 1.d-4;
% current values for x, f(x) and g(x)
dx = 1.d10;
x_cur = x0;
f_cur = feval(f,x_cur); numf = 1;
g_cur = feval(g,x_cur); numg = 1;
% Hist
it_count = 0;
ithist = zeros(maxit,4);
ithist(1,1) = norm(g_cur);
ithist(1,2) = f_cur;
ithist(1,3) = 0;
ithist(1,4) = it_count;
%% ITERATIONS
while((norm(g_cur) > tolGrad) & (dx > tolPos) & (it_count <= maxit))
it_count = it_count+1;
% Compute alpha0 and the values to test Armijo condition
alpha_cur = min(1,100/(1+norm(g_cur))); %fixup for very long steps
x_test = x_cur-alpha_cur*g_cur;
phi_alpha = feval(f,x_test); numf=numf+1;
phi_0 = f_cur;
phiP_0 = -g_cur'*g_cur;
% Test Armijo Condition
it_armijo = 0;
while( phi_alpha > (phi_0 + c1*alpha_cur*phiP_0) )
% test # of iterations
it_armijo = it_armijo+1;
if(it_armijo > 10)
disp(' Armijo error in steepest descent')
histout = ithist(1:it_count,:); costdata=[numf, numg, numh];
return;
end
% compute new alpha with step length selection
if it_armijo == 1 % quadratic model
alpha = stepLengthSelection(phi_0, phiP_0, alpha_cur, phi_alpha);
else % cubic model
alpha = stepLengthSelection(phi_0, phiP_0, alpha_cur, phi_alpha, alpha_prev, phi_prev);
end
% store previous values
phi_prev = phi_alpha;
alpha_prev = alpha_cur;
alpha_cur = alpha;
% calculate new value to test for Armijo
x_test = x_cur-alpha_cur*g_cur;
phi_alpha = feval(f,x_test); numf = numf+1;
end
% Update new values x, f(x) and g(x)
dx = norm( x_test - x_cur );
x_cur = x_test;
f_cur = phi_alpha;
g_cur = feval(g,x_cur); numg = numg+1;
% Update Hist
ithist(it_count,1) = norm(g_cur);
ithist(it_count,2) = f_cur;
ithist(it_count,3) = it_armijo;
ithist(it_count,4) = it_count;
end
%% OUTPUTS
x = x_cur;
histout = ithist(1:it_count,:);
costdata = [numf, numg];
end
function [alpha]=stepLengthSelection(phi_0, phiP_0, alpha_cur, phi_alpha, alpha_prev, phi_prev)
% Cubic/quadratic polynomial step length selection
%
% Finds minimizer alpha of the quadratic (first stepsize reduction)
% or the cubic (higher stepsize reduction) polynomial Phi on the interval
% [blow * alpha_cur, bhigh * alpha_cur]
%
%% INIT
bhigh=.5; blow=.1;
alpha_min = alpha_cur*blow;
alpha_max = alpha_cur*bhigh;
%% STEP LENTH SELECTION
if nargin == 4 % quadratic model
alpha = - phiP_0*alpha_cur^2/(2*(phi_alpha - phi_0 - phiP_0*alpha_cur) );
if alpha < alpha_min alpha = alpha_min; end
if alpha > alpha_max alpha = alpha_max; end
else % cubic model
M1 = [alpha_prev^2, -alpha_cur^2; -alpha_prev^3, alpha_cur^3];
M2 = [phi_alpha - phi_0 - phiP_0*alpha_cur; phi_prev - phi_0 - phiP_0*alpha_prev];
M = M1 * M2 / (alpha_prev^2 * alpha_cur^2 * (alpha_cur - alpha_prev));
a = M(1);
b = M(2);
alpha = (-b + sqrt(b*b - 3*a*phiP_0))/(3*a);
if alpha < alpha_min alpha = alpha_min; end
if alpha > alpha_max alpha = alpha_max; end
end
end
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment