MATLAB Program to Find A Function Minimum

Using Marquardt's Method

by Namir Shammas

The following program uses the Marquardt method to find the minimum of a function.

The function marquardt has the following input parameters:

N - number of variables
X - array of initial guesses
alpha - Alpha factor (usually has an initial value of 1E+4)
C - two-element array of coefficients (typical values are [0.1 10])
gradToler - tolerance for the norm of the slopes
FxToler - 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 100 iterations, a function tolerance of 1e-7, and a gradient tolerance of 1e-7: The values for the alpha and C parameters are supplied as the recommended default of 1E+4 and [0.1 10].

>> [X,F,Iters] = marquardt(2, [0 0], 1e+4, [0.1 10], 1e-7, 1e-7, 100, 'fx1')

X =

2.0000 -5.0000

F =

Iters =

Here is the MATLAB listing:

function y=fx1(X, N)
  y = 10 + (X(1) - 2)^2 + (X(2) + 5)^2;
end

function [X,F,Iters] = marquardt(N, X, alpha, C, gradToler, FxToler, MaxIter, myFx)
% Function MARQUARDT performs multivariate optimization using the
% Marquardt's method.
%
% Input
%
% N - number of variables
% X - array of initial guesses
% alpha -  Alpha factor
% C - two-element array of coefficients
% gradToler - tolerance for the norm of the slopes
% FxToler - array of tolerance values for the functions
% MaxIter - maximum number of iterations
% myFx - name of the optimized function
%
% Output
%
% X - array of optimized variables
% F - function value at optimum
% Iters - number of iterations
%

alpha_max = 1e+100;
if C(1) <= 0 || C(1) >= 1
  C(1) = 0.1;
end

if C(2) <= 1
  C(2) = 10;
end

bGoOn = true;
Iters = 0;
f2 = feval(myFx, X, N);

while bGoOn

  Iters = Iters + 1;
  if Iters > MaxIter
    break;
  end

  f1 = f2;

  g = FirstDerivatives(X, N, myFx);
  fnorm = norm(g);
  if fnorm < gradToler
    break;
  end

  X2 = X;

  J = SecondDerivatives(X, N, myFx);

  for ii=1:MaxIter

    Iden = eye(N, N);
    DeltaX = g / (J + alpha * Iden);
    X = X - DeltaX;

    f2 = feval(myFx, X, N);

    if f2 < f1
      alpha = C(1) * alpha;
      break;
    else
      alpha = C(2) * alpha;
      X = X2;
      if alpha >= alpha_max
        break;
      end
    end
  end
end

F = feval(myFx, X, N);

% end

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

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

% end

function FirstDerivX = FirstDerivatives(X, N, myFx)

for iVar=1:N
  xt = X(iVar);
  h = 0.01 * (1 + abs(xt));
  X(iVar) = xt + h;
  fp = feval(myFx, X, N);
  X(iVar) = xt - h;
  fm = feval(myFx, X, N);
  X(iVar) = xt;
  FirstDerivX(iVar) = (fp - fm) / 2 / h;
end

% end

function SecondDerivX = SecondDerivatives(X, N, myFx)

for i=1:N
  for j=1:N
    % calculate second derivative?
    if i == j
      f0 = feval(myFx, X, N);
      xt = X(i);
      hx = 0.01 * (1 + abs(xt));
      X(i) = xt + hx;
      fp = feval(myFx, X, N);
      X(i) = xt - hx;
      fm = feval(myFx, X, N);
      X(i) = xt;
      y = (fp - 2 * f0 + fm) / hx ^ 2;
    else
      xt = X(i);
      yt = X(j);
      hx = 0.01 * (1 + abs(xt));
      hy = 0.01 * (1 + abs(yt));
      % calculate fpp;
      X(i) = xt + hx;
      X(j) = yt + hy;
      fpp = feval(myFx, X, N);
      % calculate fmm;
      X(i) = xt - hx;
      X(j) = yt - hy;
      fmm = feval(myFx, X, N);
      % calculate fpm;
      X(i) = xt + hx;
      X(j) = yt - hy;
      fpm = feval(myFx, X, N);
      % calculate fmp
      X(i) = xt - hx;
      X(j) = yt + hy;
      fmp = feval(myFx, X, N);
      X(i) = xt;
      X(j) = yt;
      y = (fpp - fmp - fpm + fmm) / (4 * hx * hy);
    end
    SecondDerivX(i,j) = y;
  end
end

%end

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

  MaxIt = 100;
  Toler = 0.000001;

  iter = 0;
  bGoOn = true;
  while bGoOn
    iter = iter + 1;
    if iter > MaxIt
      lambda = 0;
      break
    end

    h = 0.01 * (1 + abs(lambda));
    f0 = myFxEx(N, X, D, lambda, myFx);
    fp = myFxEx(N, X, D, lambda+h, myFx);
    fm = myFxEx(N, X, D, lambda-h, myFx);
    deriv1 = (fp - fm) / 2 / h;
    deriv2 = (fp - 2 * f0 + fm) / h ^ 2;
    diff = deriv1 / deriv2;
    lambda = lambda - diff;
    if abs(diff) < Toler
      bGoOn = false;
    end
  end

% end

BACK