lsqnonlin initial conditions (transcendental equation)

4 次查看（过去 30 天）
bil
bil 2023-3-20
编辑： Matt J 2023-3-20
Ran in:
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)
First-order Norm of Iteration Func-count Resnorm optimality Lambda step 0 5 30.8735 4.5e+05 0.01 1 13 25.0777 2.5e+05 10 0.0800468 2 18 22.909 1.49e+09 1 0.0593238 3 25 22.8683 2.07e+04 100 0.0010927 4 30 22.531 1.21e+05 10 0.00652635 5 35 22.5103 3.04e+08 1 0.0514446 6 45 22.3993 7.89e+04 100000 0.000507316 7 50 22.3525 3.14e+04 10000 0.000157084 8 55 22.3327 1.71e+04 1000 0.0005333 9 60 22.326 7.43e+03 100 0.00057971 10 65 22.3236 3.02e+08 10 0.00486204 11 75 22.3235 1.49e+03 1e+06 5.98625e-05 12 81 22.3235 583 1e+07 3.31727e-07 13 86 22.3235 224 1e+06 3.21627e-07 14 91 22.3235 117 100000 3.1936e-06 15 96 22.3233 327 10000 2.66618e-05 16 101 22.3227 1.21e+03 1000 9.88313e-05 17 106 22.3224 6.84e+08 100 0.000484278 18 115 22.3223 3.02e+08 1e+06 6.04123e-05 19 126 22.3223 79.5 1e+12 1.54083e-08 20 129 7.94303e-13 Local minimum possible. lsqnonlin stopped because the relative size of the current step is less than the value of the step size tolerance.
k =
0.0010 + 0.0000i 1.0263 + 0.0021i 4.9869 - 0.0016i 2.1974 + 0.0242i
resnorm = 22.3223
residual =
Columns 1 through 10 0.0866 + 0.0926i -0.2591 + 0.0924i 0.2009 - 2.1180i -0.1175 + 0.1314i -0.1175 + 0.1314i -0.9598 - 1.5939i 0.4355 + 0.0215i 0.4355 + 0.0215i 0.4276 - 0.0045i 0.4276 - 0.0045i Columns 11 through 20 0.1274 + 0.0405i 0.1274 + 0.0405i 0.0000 + 0.0000i -0.1159 - 0.0532i -0.3642 - 0.1519i -0.3642 - 0.1519i -0.3532 - 0.1716i -0.3532 - 0.1716i 0.7534 - 1.4662i 0.7534 - 1.4662i Columns 21 through 27 0.0908 - 0.1256i -0.1557 - 1.8665i -0.1557 - 1.8665i 0.2056 - 0.1334i 0.2056 - 0.1334i -0.0721 - 0.0905i -0.3719 - 0.1255i
exitflag = 4
output = struct with fields:
iterations: 20 funcCount: 129 stepsize: 7.9430e-13 cgiterations: [] firstorderopt: 79.5088 algorithm: 'levenberg-marquardt' message: 'Local minimum possible.↵lsqnonlin stopped because the relative size of the current step is less than↵the value of the step size tolerance.↵↵<stopping criteria details>↵↵Optimization stopped because the relative norm of the current step, 1.432370e-13,↵is less than options.StepTolerance = 1.000000e-12.' bestfeasible: [] constrviolation: []
Thanks.
  2 个评论
Matt J
Matt J 2023-3-20
编辑：Matt J 2023-3-20
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:
https://www.mathworks.com/help/optim/ug/complex-numbers-in-optimization-toolbox-solvers.html
bil
bil 2023-3-20
Hi, Could it be the case that my function is still analytical? Regardless, do you have any suggestion regarding how else I could fit the parameters of the objective function to the data set?

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Matt J
Matt J 2023-3-20
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 个评论
显示 2更早的评论
bil
bil 2023-3-20
编辑：bil 2023-3-20
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.
Matt J
Matt J 2023-3-20
编辑：Matt J 2023-3-20
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

请先登录，再进行评论。

更多回答（1 个）

John D'Errico
John D'Errico 2023-3-20
Ran in:
This does not make sense.
x_data = importdata("x_data.mat");
x_data
x_data = 1×27
1.0e-42 * 0.7175 0.7175 0.7175 0.7175 0.7175 0.7175 0.7175 0.7175 0.7175 0.7175 0.7175 0.7175 0.7175 0.7175 0.7175 0.7175 0.7175 0.7175 0.7175 0.7175 0.7175 0.7175 0.7175 0.7175 0.7175 0.7175 0.7175
Your x_data vector is identically constant values.
unique(x_data)
ans = 7.1746e-43
Why do you expect any model to be fit here?
  1 个评论
bil
bil 2023-3-20
编辑：bil 2023-3-20
Hi,
The transcendental equation should have multiple solutions for the same value of x. The simplest example is how can have multiple (actually infinitely many) solutions x for a single value of .

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Dialog Boxes 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by