Why am I getting the same result as the Initial guess when using lsqnonlin?

Hello,

I am using the nonlinear least square function (lsqnonlin) to determine parameter values of power spectral density fitting curve, but it is always giving the same result with my initial guess. The number of parameters to be determine is 7 as shown in the program below;

clear
load LWAL
t=1:length(y);
z=fft(y);
pyy=z.*conj(z)/length(y);
f1=4*(0:fix(length(pyy)/2)-1)/length(pyy);
f1=1./f1;

pyy2=(1:fix(length(pyy)/2));
loglog(f1,pyy(1:fix(length(pyy)/2)),'linewidth',1.5),title('power spectral density'); grid on
pyy2=pyy2(2:32767);
f1=f1(2:32767);

fun=@(x) (x(1)*((f1.^2+x(2)*f1.^3+x(3))/(f1.^4+x(4)*f1.^3+x(5)*f1.^2+x(6)*f1+x(7)))-pyy2);
x0=[0.2476 -3.4562 5.4871 3.5951 5.2421 -1.0198 0.098];

options = optimoptions(@lsqnonlin,'Algorithm','levenberg-marquardt','checkgradients',true);
for i=1:50
x=lsqnonlin(fun,x0);
end
pyy=pyy(4:65535);
for nn=1:32766
    y1(nn)=x(1)*(f1(nn).^2+x(2)*f1(nn).^3+x(3))/(f1(nn).^4+x(4)*f1(nn).^3+x(5)*f1(nn).^2+x(6)*f1(nn)+x(7));
end

f1=1./f1;
loglog(f1,pyy(1:fix(length(pyy)/2)),'linewidth',1.5),title('power spectral density');
hold on
loglog(f1,y1(1:fix(length(pyy)/2)),'linewidth',1.5),title('power spectral density')
xlabel('Spatial wavelength m^{-1}'),ylabel('PSD/(mm^2.m)');
grid on

[Ameer Hamza]

Why are you running lsqnonlin 50 times using the same initial guess? When I run you code, lsqnonlin does not give the same result as your initial guess.

[Umar Faruk Aminu]

Hi Ameer,

Thank you for taking your time to check the program.

The lsqnonlin only change the first and second values but the remaining 5 are the same as the initial guess. Please check it again

Best regards

[Adam Danz]

You're not storing the outputs on each iteration.

x = nan(50,7); 
for i=1:50
    x(i,:) =lsqnonlin(fun,x0);
end

Of course you'll have 50 of the same result since you're using the same initial guesses. Maybe you're intention is to jitter the initial guesses to ensure you've reached a global minima. That involves adding/subtracting a tiny bit from each initial guess on each iteration of the loop.

The first two parameters are different than you initial guesses but the others are not. That suggests that initial guesses 3:7 are already at a local (or global?) minima.

Answer 1

[Ameer Hamza]

Ah!!! What an unfortunately trivial mistake

fun=@(x) (x(1)*((f1.^2+x(2)*f1.^3+x(3))./(f1.^4+x(4)*f1.^3+x(5)*f1.^2+x(6)*f1+x(7)))-pyy2);
                                      %^ use a dot here

The equation in your current code is not doing what you might expect it to do. I guess you want element-wise division, so you need to use (./) operator.

Although it still does not converge to an optimal solution with this initial guess, you will need to see how you can make a good initial guess, or consider using ga() from global optimization toolbox, that does not require an initial guess.

Edit:

I spend some time to understand your code, and I found that lsqnolin was working fine, the problem was the way you was plotting the results. You took the reciprocal of f vector at the end, which was not needed. I refactored your code. Try this now

clear
load LWAL
t=1:length(y);
z=fft(y);
pyy=z.*conj(z)/length(y);
f1=4*(0:fix(length(pyy)/2)-1)/length(pyy);
f1=1./f1;
pyy = pyy(1:fix(length(pyy)/2));

f1 = f1(2:end);
pyy = pyy(2:end);
loglog(f1,pyy,'linewidth',1.5),title('power spectral density'); grid on

x0=[0.2476 -3.4562 5.4871 3.5951 5.2421 -1.0198 0.098];
fun = @(x) x(1)*((f1.^2+x(2)*f1.^3+x(3))./(f1.^4+x(4)*f1.^3+x(5)*f1.^2+x(6)*f1+x(7)));
obj_fun = @(x) (fun(x)-pyy.');
options = optimoptions(@lsqnonlin,'Algorithm','levenberg-marquardt','checkgradients',true);
x = lsqnonlin(obj_fun,x0, [], [], options);
y1 = fun(x);

% f1=1./f1; % no need here
loglog(f1,pyy,'linewidth',1.5),title('power spectral density');
hold on
loglog(f1,y1,'linewidth',1.5),title('power spectral density')
xlabel('Spatial wavelength m^{-1}'),ylabel('PSD/(mm^2.m)');
grid on

[Umar Faruk Aminu]

Hi Ameer,

Thank you for your time,

I tried the element-wise division, but it didn't give me a good curve fitting. Can you please explain more about the ga() from the optimization toolbox?

Thank you

[Ameer Hamza]

This curve-fitting problem is proving to be quite difficult to solve. It seems to have several local minima, which make it difficult for optimizer to find a good optimal solution. Initial guess also has an effect on the final solution.

I have modified the code to work with ga(), a global optimizer (Genetic algorithm) provided by the global optimization toolbox. It does not require an initial guess and can have a good performance. Try the following code. It was taking quite a long time to run, so I terminated it. You can let this code run and wait for it to complete. Hopefully, it will give a good optimal solution. If there is still an issue, then you can let me know in the comment. Maybe, I can suggest something else if this failed too.

clear
load LWAL
t=1:length(y);
z=fft(y);
pyy=z.*conj(z)/length(y);
f1=4*(0:fix(length(pyy)/2)-1)/length(pyy);
f1=1./f1;

pyy2=(1:fix(length(pyy)/2));
loglog(f1,pyy(1:fix(length(pyy)/2)),'linewidth',1.5),title('power spectral density'); grid on
pyy2=pyy2(2:32767);
f1=f1(2:32767);

fun = @(x) (x(1)*((f1.^2+x(2)*f1.^3+x(3))./(f1.^4+x(4)*f1.^3+x(5)*f1.^2+x(6)*f1+x(7)))-pyy2);
obj_fun = @(x) norm(fun(x));

x = ga(obj_fun, 7);

pyy=pyy(4:65535);
for nn=1:32766
    y1(nn)=x(1)*(f1(nn).^2+x(2)*f1(nn).^3+x(3))/(f1(nn).^4+x(4)*f1(nn).^3+x(5)*f1(nn).^2+x(6)*f1(nn)+x(7));
end
% y1 = fun(x)

f1=1./f1;
loglog(f1,pyy(1:fix(length(pyy)/2)),'linewidth',1.5),title('power spectral density');
hold on
loglog(f1,y1(1:fix(length(pyy)/2)),'linewidth',1.5),title('power spectral density')
xlabel('Spatial wavelength m^{-1}'),ylabel('PSD/(mm^2.m)');
grid on

[Umar Faruk Aminu]

Hi Ameer,

Thank you for the suggestion.

I run the modified code but it didn't give me an optimal solution and the curve obtained did not fit the data. I observed when i run my program which is attached above, it is showing me

'solver stopped prematurely

lsqnonlin stopped because it exceeds function evaluation limit

options.MaxFunctionEvaluations=700'

is this the reason why i am not getting a good curve fitting solution?

How can I solve this problem?

Thank you

Best regards

[Ameer Hamza]

Check the code in my updated answer. I have pointed out the issue with you current approach.

[Umar Faruk Aminu]

Hi Ameer,

Thank you very much for taking your time to check the code again. I am now getting a meaningful solution

Best wishes

[Ameer Hamza]

I am glad to be of help.