MATLAB Program to Find A Function Minimum

Using the Conjugate Gradient Method

by Namir Shammas

The following program calculates the minimum point of a multi-variable function using the Fletcher-Reeves conjugate gradient method.

The function conjgrad has the following input parameters:

N - number of variables
X - array of initial guesses
Eps_Fx - tolerance for function
MaxIter - maximum number of iterations
myFx - name of the optimized function

The function generates the following output:

X - array of optimized variables
F - function value at optimum
Iters - number of iterations

Here is a sample session to find the optimum for the following function:

y = 10 + (X(1) - 2)^2 + (X(2) + 5)^2

The above function resides in file fx1.m. The search for the optimum 2 variables has the initial guess of [0 0]. The search employs a maximum of 1000 iterations and a function tolerance of 1e-7 :

>> [X,LastF,Iters]=conjgrad(2,[0 0],1e-7,1000,'fx1')

X =

2.0000 -5.0000

LastF =

Iters =

Here is the MATLAB listing:

function [X,LastF,Iters]=conjgrad(N,X,Eps_Fx,MaxIter,myFx)
% Function performs multivariate optimization using the
% Hooke-Jeeves search method.
%
% Input
%
% N - number of variables
% X - array of initial guesses 
% Eps_Fx - tolerance for diffence in successive function values
% MaxIter - maximum number of iterations
% myFx - name of the optimized function
%
% Output
%
% X - array of optimized variables
% LastF - function value at optimum
% Iters - number of iterations
%

initStep = 0.1;
minStep = 0.000001;
LastF = feval(myFx, X, N);
[dfnorm,Deriv] = getgradients(X, N, myFx);

lambda = 0;
lambda = linsearch(X, N, lambda, initStep, minStep, Deriv, myFx);
X = X + lambda * Deriv;

bGoOn = true;
Iters = 0;

while bGoOn

  Iters = Iters + 1;
  if Iters > MaxIter
    break;
  end
  
  dfnormOld = dfnorm;
  DerivOld = Deriv;
  [dfnorm,Deriv] = getgradients(X, N, myFx);
  Deriv = (dfnorm / dfnormOld)^2 * DerivOld - Deriv;
  if dfnorm < Eps_Fx
    break;
  end
  lambda = 0;
  lambda = linsearch(X, N, lambda, initStep, minStep, Deriv, myFx);
  X = X + lambda * Deriv;
  F = feval(myFx, X, N);
  if abs(F - LastF) < Eps_Fx
    bGoOn = false;
  else
    LastF = F;
  end  
  
end

% end

function y = myFxEx(N, X, DeltaX, lambda, myFx)

  X = X + lambda * DeltaX;
  y = feval(myFx, X, N);

% end

function [fnorm,Deriv] = getgradients(X, N, myFx)

  for i=1:N
    xx = X(i);
    h = 0.01 * (1 + abs(xx));
    X(i) = xx + h;
    Fp = feval(myFx, X, N);
    X(i) = xx - h;
    Fm = feval(myFx, X, N);
    X(i) = xx;
    Deriv(i) = (Fp - Fm) / 2 / h;
  end
  fnorm = norm(Deriv);
% end  

function lambda = linsearch(X, N, lambda, initStep, minStep, D, myFx)

  f1 = myFxEx(N, X, D, lambda, myFx);
  while initStep > minStep
    f2 = myFxEx(N, X, D, lambda + initStep, myFx)  ;
    if f2 < f1
      f1 = f2;
      lambda = lambda + initStep;
    else  
      f2 = myFxEx(N, X, D, lambda - initStep, myFx);
      if f2 < f1
        f1 = f2;
        lambda = lambda - initStep;
      else
        % reduce search step size
        initStep = initStep / 10;
      end
    end
  end 
% end

BACK