lsqnonlin initial conditions (transcendental equation)

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.

[Matt J]

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]

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?

Answer 1

[Matt J]

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)

[bil]

Thank you for the suggestion, the problem still remains where the final fitted parameters heavily depend upon my choice of initial parameters. This behavior confuses me and I am not sure how I should probe it (changing tolerance parameters does not help alleviate the "local minima possible" message either).

[Torsten]

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".

[bil]

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]

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

Answer 2

[John D'Errico]

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?

[bil]

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 .