Optimizing systems of non-linear equation using optimvar and lsqcurvefit

Hi there,

I am currently trying to optimise a set of non-linear equations iteratively - I am trying to solve equation (4) in this paper (https://aip.scitation.org/doi/pdf/10.1063/1.3416910?casa_token=i26EckyYkb8AAAAA:UQ1axvJNde0OswFL9jgt24izc_sVX6XhurCYzC0M1BfNObq-K6OOOjuykKQuwt37mR6zALj-QrhD) where the author uses a non-linear least square optimisation approach to solve for 4 variables with 3 non-linear equations. Simply put, the 3 equations can be represented by the following:

1: y1 = sqrt(((1+4.*xdata.^2.*x(3)^2).^(-2) + x(4))./((1+4.*xdata.^2.*x(2)^2).^(-2) + x(4)));

2: y2 = sqrt(((1+4.*xdata.^2.*x(3)^2).^(-2) + x(4))./((1+4.*xdata.^2.*x(1)^2).^(-2) + x(4)));

3: y3 = sqrt(((1+4.*xdata.^2.*x(1)^2).^(-2) + x(4))./((1+4.*xdata.^2.*x(2)^2).^(-2) + x(4)));

where y1, y2, y3, are the respective y values and xdata the x values that I know. Essentially, I have 3 sets of data points that I want to fit with the 4 unknown variables (x(1) to x(4)) through optimisation. I have tried the following:

x = optimvar('x',4);
eq1 = y1 == sqrt(((1+4.*xdata.^2.*x(3)^2).^(-2) + x(4))./((1+4.*xdata.^2.*x(2)^2).^(-2) + x(4)));
eq2 = y2 == sqrt(((1+4.*xdata.^2.*x(3)^2).^(-2) + x(4))./((1+4.*xdata.^2.*x(1)^2).^(-2) + x(4)));
eq3 = y3 == sqrt(((1+4.*xdata.^2.*x(1)^2).^(-2) + x(4))./((1+4.*xdata.^2.*x(2)^2).^(-2) + x(4)));
prob = eqnproblem;
prob.Equations.eq1 = eq1;
prob.Equations.eq2 = eq2;
prob.Equations.eq3 = eq3;
show(prob)
x1.x = [0.5e-3 0.5e-3 0.5e-3 0];
[sol,fval,exitflag] = solve(prob,x1);
disp(sol.x)

But these are the comments generated:

Solving problem using fsolve.
Warning: Trust-region-dogleg algorithm of FSOLVE cannot handle non-square systems; using Levenberg-Marquardt algorithm
instead. 
> In fsolve (line 323)
  In optim.problemdef.EquationProblem/callSolver
  In optim.internal.problemdef.ProblemImpl/solveImpl
  In optim.problemdef.EquationProblem/solve 

No solution found.

fsolve stopped because the relative size of the current step is less than the
value of the step size tolerance, but the vector of function values
is not near zero as measured by the value of the function tolerance.

I was able to fit a single equation using the lsqcurvefit function by defining the following:

fun = @(x,xdata) sqrt(((1+4.*xdata.^2.*x(1)^2).^(-2) + x(3))./((1+4.*xdata.^2.*x(2)^2).^(-2) + x(3)));
x0 = [0.1e-3 0.1e-3 0];
options = optimoptions('lsqcurvefit','Algorithm','levenberg-marquardt');
x = lsqcurvefit(fun,x0,xdata,y1,[],[],options)

However, when I tried defining the the 3 equations with a single function handle with the following codes:

fun = @(x,xdata) [sqrt(((1+4.*xdata.^2.*x(1)^2).^(-2) + x(4))./((1+4.*xdata.^2.*x(2)^2).^(-2) + x(4)));
    sqrt(((1+4.*xdata.^2.*x(1)^2).^(-2) + x(4))./((1+4.*xdata.^2.*x(3)^2).^(-2) + x(4)));
    sqrt(((1+4.*xdata.^2.*x(2)^2).^(-2) + x(4))./((1+4.*xdata.^2.*x(3)^2).^(-2) + x(4)))];

x0 = [0.5e-3 0.5e-3 0.5e-3 0.1e-5];
options = optimoptions('lsqcurvefit','Algorithm','levenberg-marquardt');
x = lsqcurvefit(fun,x0,xdata,y1,y2,y3,[],[],options)

It shows the following:

Warning: Length of lower bounds is > length(x); ignoring extra bounds. 
> In checkbounds (line 27)
  In lsqnsetup (line 77)
  In lsqcurvefit (line 210) 
Warning: Length of upper bounds is > length(x); ignoring extra bounds. 
> In checkbounds (line 41)
  In lsqnsetup (line 77)
  In lsqcurvefit (line 210) 
Exiting due to infeasibility:  4 lower bounds exceed the corresponding upper bounds.

Unfortunately, I am not too familiar with the optimisation procedures with MATLAB, so am wondering if anyone has any insights or suggestions on what else I should try. Thank you for your time.

Best Regards,

Wayne

[Alex Sha]

Usually, if the No. of variables is bigger than the No. of non-linear equations, there will not be accurate results。Moreover, give out the values of y and xdata you known, so others may have a try.

[Wayne]

Thanks Alex for your comments. I was able to do the fit but the results are not very accurate, as you have said. I will be asking this in a separate post and will include the relevant data there.

Answer 1

[Star Strider]

Your ‘fun’ function will work, providing the dependent variable matrix is (3xN), and ‘y1’. ‘y2’ and ‘y3’ are row vectors.

Assuming that, the dependent variables (as row vectors) need to be in a matrix, so that there is one variable for ‘y’, not three.

Try this:

ymtx = [y1; y2; y3];
x = lsqcurvefit(fun,x0,xdata,ymtx,[],[],options)

The error arises because in the code youi posted, ‘y2’ and ‘y3’ appear in the positions in the lsqcurvefit argument list reserved for the bounds vectors.

.

[Wayne]

Thank you - that works.

[Star Strider]

As always, my pleasure!

Answer 2

[Matt J]

Avoid sqrt so as to ensure fun is differentiable

fun = @(x,xdata) [((1+4.*xdata.^2.*x(1)^2).^(-2) + x(4))./((1+4.*xdata.^2.*x(2)^2).^(-2) + x(4));
       ((1+4.*xdata.^2.*x(1)^2).^(-2) + x(4))./((1+4.*xdata.^2.*x(3)^2).^(-2) + x(4));
       ((1+4.*xdata.^2.*x(2)^2).^(-2) + x(4))./((1+4.*xdata.^2.*x(3)^2).^(-2) + x(4))];

x0 = [0.5e-3 0.5e-3 0.5e-3 0.1e-5];
options = optimoptions('lsqcurvefit','Algorithm','levenberg-marquardt');

xdata=xdata(:).';
ydata=( [y1(:),y2(:),y3(:)].' ).^2;

x = lsqcurvefit(fun,x0,xdata,ydata,[],[],options)

[Wayne]

Hi Matt, thanks for your help! Is there a reason why using sqrt may cause the function to be not differentiable?

[Matt J]

As a simpler example, consider the 1D least squares cost function,

It's quite clear that the first two terms on the far right hand side are non-differentiable at x=0 is it not? More importantly, the lsqcurvefit algorithms would try to compute the Jacobian of the residuals which is also non-differentiable at x=0.