Avoiding local minimum with fmincon

Hi, i'm doing a non linear optimization using the fmincon (interior point) function with matlab, but i'm stucking in local minimum. I also did this optimization using the gradient descent and together with it, I implemented a heuristic to avoid local minimum, which was:

if (actual objective function value > anterior objective function value)

step of gradient descent reduces in 10%

actual objective function increases in a little step.

This implementation returned me good results.

I've already done this for trying to use this heuristic with fmincon, but I didn't get out of the local minimum. So my question is: Is there any way to avoid local minimum using fmincon or fminunc?

PS: I've already done the GlobalSearch and MultiStart implementation and the results were worse thatn the gradient descent.

[Matt J]

Could you explain your terminology: "actual objective function value" and "anterior objective function value" ? Also, does "gradient descent" mean "steepest descent" here?

[Felippe Trigueiro Angelo]

Hi Matt, the actual value is the value of the objective function in the actual iteration. For exemple:

In the first iteration I get the value 1.23 in my objective function and in the second I get the value 1.5, so in this case I will have to do the heuristic above.

I don't know the steepest descent algorithm, but the gradient descent algorithm just actualizes the parameters with the negative of the gradient of the objective function multiplied with a step.

[Matt J]

So whenever the algorithm takes a non-descending step, you dampen the step-size? I don't see why that would help to avoid local minima. If the algorithm tries to climb uphill to escape a local minimum, your heuristic will pull it back down.

[Felippe Trigueiro Angelo]

It is not just descreasing the step, but also changing the objective function value in a little value to force the algorithm finding the minimum.

[Matt J]

No, I'm afraid I still don't really have a clear picture of the method. The iterations of steepest/gradient descent are entirely guided by the objective's gradient and the step-size. The value of the objective at each iteration doesn't figure into it, unless maybe you are talking about some sort of line search that you are doing to help choose the step-size. But even then, its not clear how perturbing the objective values during a line search would be helpful for climbing out of a local minimum.

It is also still not clear how you would have applied your heuristic to fmincon, as you claim to have done. fmincon doesn't offer the user any options for controlling/altering its step-size selections.

[Felippe Trigueiro Angelo]

Here is the code I'm using. It is part of the Medical Image Registration Toolkit by Adriy Myronenko. The while part is the heuristic he implemented and that I'm using.

Just for a clarification, gamma is the step size of the gradient descent and the anneal is the percentage that it is reduced, in this case I'm testing with 0.9.

% MIRT3D_REGISTERATION  non-rigid registration 
% of a single pair of 3D images at a given hierarchical level

% Copyright (C) 2007-2010 Andriy Myronenko ([e-mail])
% also see http://www.bme.ogi.edu/~myron/matlab/MIRT/
%
% This file is part of the Medical Image Registration Toolbox (MIRT).


function [X, im]=mirt3D_registration(X, main, optim)

%% normalize the initial optimization step size
% compute the gradient of the similarity measure
[Xx,Xy,Xz]=mirt3D_nodes2grid(X, main.F, main.okno);
[f, ddx, ddy,ddz]=mirt3D_similarity(main, Xx, Xy,Xz);
% divide the initial optimization step size by the std of the gradient
% this somewhat normalizes the initial step size for different possible
% similarity measures used
optim.gamma=optim.gamma/std([ddx(:); ddy(:); ddz(:)],0);
clear ddx ddy ddz Xx Xy Xz;

%% Start the main registration
% compute the objective function and its gradient
[fe, T, im]=mirt3D_grad(X,  main);              % compute the similarity measure and its gradient
[Xp, fr]=mirt3D_regsolve(X,T,main, optim, 1);   % compute the regularization term and update the transformation
f=fe+fr;                                        % compute the value of the total objective function (similarity + regularization)

fchange=optim.fundif+1; % some big number
iter=0;

% do while the relative function difference is below the threshold and
% the meximum number of iterations has not been reached
while (abs(fchange)>optim.fundif) && (iter<optim.maxsteps)
                
        % find the new positions of B-spline control points,
        % given their currect positions (X) and gradient in (T)
        [Xp, fr]=mirt3D_regsolve(X,T,main, optim);
               
        % compute new function value and new gradient
        [fe, Tp, imb]=mirt3D_grad(Xp,  main);
        fp=fe+fr;
        
        % compute the relative objective function change
        fchange=(fp-f)/f;
        % check if the step size is appropriate
        if ((fp-f)>0),
            % if the new objective function value does not decrease,
            % then reduce the optimization step size and
            % slightly increase the value of the objective function 
            % (this is an optimization heuristic to avoid some local minima)
            optim.gamma=optim.gamma*optim.anneal;
            f=f+0.001*abs(f);
        else
            % if the new objective function value decreased
            % then accept all the variable as true and show the progress
            X=Xp; f=fp; T=Tp; im=imb;
            
           % mesh_epsplot(X(:,:,round(end/2),1),X(:,:,round(end/2),2)); drawnow;
            % show progress
           disp([upper(main.similarity) ' ' num2str(f) ' dif ' num2str(fchange) ' sub ' num2str(main.level) ' cyc ' num2str(main.cycle) ' vol ' num2str(main.volume) ' iter ' num2str(iter) ' gamma = ' num2str(optim.gamma)]);
        end;
      iter=iter+1;
end

Answer 1

[Matt J]

So my question is: Is there any way to avoid local minimum using fmincon or fminunc?

There is no systematic, problem-independent way. You must, through problem-specific ingenuity, choose an initial guess close enough to the global minimum so that your algorithm falls into its capture basin.

The fact that GlobalSearch and MultiStart are having such a hard time suggests to me that your cost function is extremely dense with local minima. You might want to try ga() or one of the other global solvers.