Optimizer (fminsearch) doesn't work as expected

Hi all,

I am working on an optimization problem, where I need to find the global minimum for , where

There is also a constraint or penalty term: .

To test the code, I use the exact solution as the guess, however, the optimizer cannot reproduce the exact solution. And I limit my x to [0, 20] rather than [,∞], which theoretically affects little.

Some thoughts: My variables are , but is involded in the calculation a lot, I'm wondering whether the calculation of introduces a good amount of errors, which accumulates during the calculation, and affects the results significantly.

The codes are attached below. opt_main.m is the main function and functional_elastic_main.m calculates the objective function to optimize. A figure shows the comparison between the exact solution and the solution produced by the optimizer is attached below as well. Looking at the figure, the optimizer works better for the left-size curve but works not properly for the right-size curve.

No matter what kind of ideas you have, please share those with me. I appreciate very much.

The main function:

function opt_main
%% main function
tic
clear all; clc; close all;

% ---------- user input ----------
% given parameters
mu = 1;
nu = 0.3;
vb = [1,0,0];
b = norm(vb);
lambda = 100; % lagrange multiplier 

x = [0:1:20];

% exact solution, used here as input
% it is shifted a bit for better demonstration
del0 = (b/pi*atan(2*(1-nu)*(x-10)/b)+1/2*b); 

n = length(x);

% calculate the middle point of x_i and x_i+1
% and dx between x_i and x_i+1
% and length(x1) = n-1;
[x1,dx1] = discrete_x(x);

%% --------------- optimization ---------------
opts = optimset;
opts.TolX = 1e-10;
opts.TolFun = 1e-10;
opts.MaxIter = 1000000;
opts.MaxFunEvals = 1000000;

% ----- perturbation -----
% del0(2:10) = del0(2:10) - 0.1;

objectiveFun = @(d) functional_elastic_main(n,mu,nu,x,x1,dx1,...
            d,b,lambda);
[nd2] = fminsearch(objectiveFun,del0,opts);

%% --------------- plotting the results ---------------
plot(x,del0,'--b^','linewidth',2,'MarkerSize',10);
hold on
plot(x,nd2,'-ko','linewidth',2,'MarkerSize',10)

ax = gca;
ax.XAxis.FontSize = 20;
ax.YAxis.FontSize = 20;
    
xlabel('X','Interpreter','Latex','FontSize',20,'Color','k');
ylabel('$\delta$','Interpreter','Latex','FontSize',25,'Color','k');
legend({'exact','w/ minimization'},'fontsize',16,...
    'location','southeast')
print('-dpng','-r400',sprintf('delta_elastic.png'));
%%
toc
fprintf('ALL DONE \n');
end

Objective function:

function [func] = functional_elastic_main(n,mu,nu,x,x1,dx1,del0,b,lambda)

% calculate the middle point of delta_i and delta_i+1
% calculate the slope between delta_i and delta_i+1
[del1,ddeldx1] = discrete_del(del0,x,n);

% ------- calculate U_1 -------
U_1 = 0;
u_1 = zeros(n-1,1);
for ii = 1:n-1
    temp1 = -mu/(2*pi*(1-nu))*ddeldx1(ii);
    for jj = 1:n-1
        if ii ~= jj
            temp2 = ddeldx1(jj)*log(abs(x1(ii)-x1(jj)))*dx1(jj);
            u_1(ii,1) = u_1(ii,1) + temp1*temp2; 
        end
    end
    U_1 = U_1 + u_1(ii,1)*dx1(ii);
end

% ------- calculate U_2 -------
U_2 = 0;
u_2 = zeros(n-1,1);
for ii = 1:n-1
    u_2(ii) = mu*b/(4*pi^2)*(1-cos(2*pi*del1(ii)/b));
    U_2 = U_2 + u_2(ii)*dx1(ii);
end

% ------- add the contraint or penalty term-------
R = 0;
for ii = 1:n-1
    R = R + abs(ddeldx1(ii)*dx1(ii));
end
R = lambda*(R - 1)^2;

func = U_1 + U_2 + R;

end

Answer 1

[Matt J]

It appears that you have 21 unknown variables. That is too many. fminsearch is only guaranteed to converge theoretically for one unknown variable, and empirically starts to perform poorly after about 6 variables.

You might try using ga instead, if you have the Global Optimization Toolbox.

[Joe Zhang]

Thansk for your suggestion. I tried ga, but the result is even worse and the result varies every time. Can I get better optimization by optimizing variables one by one?

[Matt J]

You can't optimize one variable at a time unless your objective is a separable sum of 1D functions.

[Joe Zhang]

Since both fminsearch and ga don't work well for this case. Do you know any other optimizers that may work for a multi-variable problem? Or do you think the problem come from any mistakes I made in my code?

[Matt J]

Well, I haven't really seen what you've done with ga(). One thing you could try is to include the fminsearch solution, which seems to be fair, in ga()'s initial population using the InitialPopulationMatrix option,

https://www.mathworks.com/help/gads/ga.html#mw_4a8bfdb9-7c4c-4302-8f47-d260b7a43e26

That way, ga() shoudln't do worse than the fminsearch solution.

[Joe Zhang]

Matt, thanks for your patience.

So, first, I don't think fminsearch() really does its job properly. I shift the middle part of the exact solution a bit and use this as my initial guess, and the reasult from fminsearch() doesn't look good, although existflag is 1. The shifting is like:

del0(18:25) = del0(18:25) - 0.05;

For ga(), I modified is the optimization part, and I also droped the penalty term in functional_elastic_main.m, added a condition. The code is following.

%% --------------- optimization ---------------
opts = optimset;
opts.TolX = 1e-10;
opts.TolFun = 1e-10;
opts.MaxIter = inf;
opts.MaxFunEvals = inf;
% opts.Display = 'iter';

lb = 0*ones(n,1);
ub = 1*ones(n,1);
nonlcon = @condition;
[nd2,fval,exitflag] = ga(objectiveFun,n,[],[],[],[],lb,ub,nonlcon,opts);

condiftion.m

function [c,ceq] = condition(d)
vb = [1,0,0];
b = norm(vb);
x = [0:1:20]/b;
n = length(x);
[~,dx1] = discrete_x(x);
[~,ddeldx1] = discrete_del(d,x,n);
r = 0;
for ii = 1:n-1
    r = r + abs(ddeldx1(ii)*dx1(ii));
end
c = [];
ceq = r-1;
end

I also tried what you suggested that give an initial population to ga(). If I understand correctly, it is the initial guess. And I use the exact solution as the initial populations. If I misunderstood what you suggested, please let me know. The result is not better than what I got from fminsearch().

%% --------------- optimization ---------------
opts = optimset;
opts.TolX = 1e-10;
opts.TolFun = 1e-10;
opts.MaxIter = inf;
opts.MaxFunEvals = inf;

opts.InitialPopulation = del0;
lb = 0*ones(n,1);
ub = 1*ones(n,1);
nonlcon = @condition;
[nd2,fval,exitflag] = ga(objectiveFun,n,[],[],[],[],lb,ub,nonlcon,opts);

Sorry I cannot attach any files due to the community limit.

[Matt J]

To test the code, I use the exact solution as the guess, however, the optimizer cannot reproduce the exact solution.

You said you initially used the exact theoretical solution as the initial guess in fminsearch. Does fminsearch manage to reduce the objective function value relative to this initial point? Does ga() manage to? If so, you know that the theoretical solution is not the global minimizer of the objective that you've provided.

[Joe Zhang]

Thanks for the suggestion. Theoretically, this is not possible, the exact solution comes from a balanced physical system. The most confusing thing is whether optimizers work properly or not. If the optimizer works fine, then I can go back and look into the derivation of the exact solution. However, the results from optimizers don't look correct or not reliable. As you mentioned, fminsearch() is likely uncapable of solving this many-variable problem. And if the setup for ga() was not wrong, then ga() didn't provide better results either. And I also tried fmincon() and fminuc().

[Matt J]

I think you may have missed my point. You need to compare the objective function value of the solution you're getting to the objective function value of a solution that you think would be better. If the solution you're getting has the lower objective value, you know it can't be the optimizers' fault. It's the fault of the objective function. The objective function disagrees with you on what a good solution is.

[Torsten]

@Joe Zhang

The surprising thing for me is that "func" returned from "functional_elastic_main" becomes negative.

This should never happen in fitting problems. The error you try to minimize should always be positive and equal zero in the best case.

I think the solution for delta you think is optimal might be optimal under all infinitly often differentiable functions, but numerically, you will get

delta(x) = 0 , x <= 9
delta(x) linear for 9 <= x <= 11
delta(x) = 1, x >= 11.

[Joe Zhang]

@Torsten

Thanks for the reply and the information. Do you believe the fact that "func" becomes negative is due to any wrong calculation in the script? And yes, theoretically, integrating from -infinity to +infinity can be more ideal, but it is infeasible. Besides, delta(x) is getting closer to 0 with x getting closer to 0, rather than being euqal to 0. Do you agree that increacing the range of integration helps?

[Torsten]

Do you believe the fact that "func" becomes negative is due to any wrong calculation in the script?

At least, I don't see that negative values for U_1 cannot happen.

Do you agree that increacing the range of integration helps?

Helps for what ? My guess is that the numerical solution you'll get (in the limit) for delta is 0 for x<= 10 , 1 for x>10. This gives u1(x) = u2(delta(x)) = 0 and thus U1 = U2 = 0. And if you take eps > 0 and define delta to be linear in between 10 - eps and 10 + eps, the integral of d(delta)/dx equals (2*eps * 1/eps)/2 = 1.

[Joe Zhang]

@Matt J

Sorry for misunderstanding your point. Yes, the total value does decrease after minimization. But it is not desired. Is there anything I can do to improve it? By the way, with a simple test , minimizing f = (x1 - 1)^2 + (x2 - 2)^2 + .... + (xn - n)^2 using fminsearch(), it turns out that fminsearch works fine with n <= 21, and when n > 21, fminserach() doesn't work properly, which proves what you claimed that fminsearch() may not work properly with so many variables.

[Matt J]

But it is not desired. Is there anything I can do to improve it?

It means your objective function design needs to be improved, but that is not a Matlab issue. It's a problem-definition issue.

You mentioned in your initial post, though, that the discretization of the integral in the problem could be creating deviations from the ideal computations. Have you tried finer sampling? Or maybe using the integral() command?

[Joe Zhang]

@Matt J

Thanks for the suggestion. I'd like to try integral(). But I have a problem writing a function handle with ddelta_dx and ln|x-y|,or ln|x_i-x_j|, knowning my variable is delta rather than x, and integrate it using integral().

[Matt J]

integral() doesn't know or care what the variable is called, e.g.,

integral(@(x) x,0,2.5)

ans = 3.1250

integral(@(delta) delta,0,2.5)

ans = 3.1250

[Joe Zhang]

@Matt J

But in my calculation of , is needed and this item dependens on my discretization of x. Also, δis my unknown, but my integrant is integrated over x, which means I need to specify x as unknowns as well, intergal() doesn't work for this case, based on I tested so far. But I think int() can work once I specify x and delta as symbols.

[Matt J]

You should probably be parametrizing delta in terms of unknown spline coeffiicients rather than explicit samples. Given the spline coefficients, the continuous-domain derivatives of delta can be obtained by various means, as discussed here,

https://www.mathworks.com/matlabcentral/answers/95194-how-do-i-find-the-derivative-of-a-spline-curve-in-matlab-7-9-r2009b

Once you have continuous-domain derivatives, you can integrate them continuously (with integral()) as well.