lsqnonlin initial conditions (transcendental equation)
Show older comments
Hey all,
I am trying to fit a transcendental equation to some data points but I am coming across the issue where lsqnonlin ends up giving solutions ("possible minima") that heavily depend upon the choice of initial condition. I am wondering if this is because my choice of initial condition simply is incorrect or if there is something else going on and how I could better probe what exactly is happening. I have attached the data_points below and my code is:
a1 = 20;
a2 = 10^-4;
x_data = importdata("x_data.mat");
y_data = importdata("y_data.mat");
options = optimoptions(@lsqnonlin,'Algorithm','levenberg-marquardt','StepTolerance',1e-12,'FunctionTolerance',1e-12,'Display','iter');
a0 = [10^-3,1,5,2]; %initial guess%
m = @(a)sqrt(a2^2 + (x_data.*a(2)).^2 + a(3).*(1./a1).^a(4));
factor1 = @(a)sqrt((x_data+1i.*(m(a)./a(2)))./(x_data-1i.*(m(a)./a(2))));
b = @(a)(1/2) - 1i.*(y_data./a(2));
factor2 = @(a)double(gamma(sym(complex(1-b(a))))./gamma(sym(complex(b(a)))));
factor3 = @(a)(2./(a(1).*sqrt(x_data.^2+(m(a).^2)./a(2).^2))).^(2.*1i.*y_data./a(2));
func = @(a)factor1(a).*factor2(a).*factor3(a) - 1i;
[k,resnorm,residual,exitflag,output] = lsqnonlin(func,a0,[],[],options)
Thanks.
2 Comments
It's pretty amazing that the optimization runs to completion at all, given that your objective is complex-valued. Generally speaking, that isn't allowed, but for exceptions discussed here:
bil
on 20 Mar 2023
Accepted Answer
Matt J
on 20 Mar 2023
You can split the model function into real and imaginary parts,
func=@(a) [real(func(a)); imag(func(a))];
[k,resnorm,residual,exitflag,output] = lsqnonlin(func,a0,[],[],options)
4 Comments
bil
on 20 Mar 2023
Torsten
on 20 Mar 2023
Do you have a well-arranged graphics of your model function ?
My imagination blocks when I see these nested functions that finally lead to "func".
Could you elaborate on what you mean by well-arranged graphic. Do you mean like the explicit function itself?
Edit: Here is the function anyway
where
and
where 
are defined as originally in my code and the fitted parameters are 
and 
correspond to the x_data and y_data respectively.
are defined as originally in my code and the fitted parameters are and correspond to the x_data and y_data respectively.When you say you are getting highly different solutions for different initial points, are you getting very different resnorms in each case as well? If not, the problem is simply ill-posed: you have non-unique solutions.
If you are getting very different resnorms, then you are probably getting stuck at local minima, in which case you could contemplate doing a 4D parameter search over the different k(i) combinations for a more accurate initial guess k0. This shouldn't be very computationally demanding, since the operations in func() are very well-vectorized. Can you write down lb,ub bounds on the parameters? If so, you can download ndgridVecs and create grid search data quite easily, e.g. as below,
ranges=arrayfun(@(l,u)linspace(l,u,30), lb,ub,'uni',0); %30-point discretization
[~,~,K1,K2,K3,K4]=ndgridVecs(1,1,ranges{:}); %grid search coordinates
values = vecnorm(func(K1,K2,K3,K4)); %rewrite func to support 4-argument input syntax
[~,i]=min(values,[],'all','linear');
k0=[K1(i), K2(i), K3(i), K4(i)]; %initial point
More Answers (1)
This does not make sense.
x_data = importdata("x_data.mat");
x_data
Your x_data vector is identically constant values.
unique(x_data)
Why do you expect any model to be fit here?
1 Comment
Categories
Find more on Creating, Deleting, and Querying Graphics Objects in Help Center and File Exchange
