c# "Nonlinear Least Squares Problem without Constraints"

[Gianluca_G_1]

Thank you very much for your patience and your helpful suggestions.

We are here again, because we are not sure about how the procedure works.

Following your on line example, we made something similar in C#.

But it seems that the algorithm not always converge, so we would like your opinion.

In the attached test, we have a set of measures, our y values.

We have inputted also initial values very near the solutions (that in this specific case we know), so we expect that the algorithm converges to something close as result.

What we obtain is a very far solution, and in many cases it never exit from the main loop.

We have tried with many Trust Region Size (rs parameter in dtrnlspbc_init) but nothing changed.

Could you help us?

Our value of the functional is:

F = y – f(x) = measured value - calculated value.

Is it correct?

Sometime we noted that if:

F = ((measured value - calculated value) / measured value)^2

Trust region method works well.

Why this differences?

Moreover (and maybe this is the crucial question), we noted that if a minimum = 0 exist and then all F = 0, the algorithm works perfectly.

Vice versa, if minimum > 0 the algorithm in the tested cases does not work, as already said.

See attached draft.

Finally consider that, we know that in our case, Jacobian matrix coefficients cannot be calculate in a very precise way because the calculated value function is an infinite integrals calculated in a numerical way, but if this could be a problem, it will be a problem also when minimum = 0 exist.    

thank you very much

 
Gianluca

[mecej4]

Please provide a mathematical description of the problem, after seeing which someone may be able to suggest things to try.

Invoke my_obj_function() with your "known solution" and see if the results are as expected. 

[Gianluca_G_1]

Here we are,

I attach a document with the matematical and phisical details. 

For us it is very important issue to understand if we are using well MKL API ... so if you prefer that we submit this issue in a more appropriate channel, tell us.

We have bought Intel Parallel Studio XE with MKL and we have a regular license.

Best Regards

Gianluca

[mecej4]

Your least squares problem is complicated by the presence of the semi-infinite range of the integral, and the oscillatory nature of the integrand in that integral (because of the Bessel function J₀).

Furthermore, for certain data values, the residuals become quite insensitive to one or more of the fit coefficients, and this fact can lead to divergence of the NLSQ algorithm(s) used. Likewise, without good starting guesses for the fit coefficients, such algorithms can fail. 

I was able to obtain converged solutions for the three data sets in your Program.cs file that correspond to the 3-layer model, assuming that B(λ) = B₃(λ). I used the Quadpack routine DQAGI for quadrature, instead of using the Simpson rule as you did.

Here are the resulting fit coefficients:

TP3  : 2.0106E+02  3.7798E+02  1.8857E+02  2.0062E+00  7.2084E+00 ||F|| =   6.197E+00
WTGA2: 4.5597E+02  1.1364E+00  7.2194E+01  2.9438E+00  2.3344E-01 ||F|| =   1.382E+02
WTGA3: 1.2529E+02  3.0021E+01  1.4700E+02  7.0832E-01  7.1831E+00 ||F|| =   4.767E+00

I realize that there may be other local minima, and that the fit coefficients that I found may be unsuitable based on physics reasons. One point to note is that the quadrature often failed when the starting guess was not good.

I note that there are some 2-layer data in Program.cs, but since you did not provide the equations for the 2-layer model, i.e., the formula for B₂(λ), I did not process those. 

[Gianluca_G_1]

Hello,

In the test code, I have provided to you, there is only an example for a three layer model.
We know that Trust Region method is quite robust and also reliable, so in the example with very good starting values, the method should converge.

Instead it seems the method find solutions worst than initial values and this is not possible.

Bessel has an oscillatory nature but with the internal exponential this oscillation should be shaded. In the following some typical curves.

(attached figure)

We observed that in theoretical condition, when minimum of error function is zero, the algorithm finds always a solution, also with bad initial conditions (see attached).

This is a very strange aspect.

Our competitors for the same problem use the Levemberg Marquardt method that of course converge. Literature said that Trust Region is powerful than Levemberg Marquardt.

Moreover, we inform you that before MKL Trust Region we used a free code based on Simplex method and we observed that this code convergence.

We think the problem could be in Jacobian calculation because we cannot calculate partial derivate in a precise way starting on a function calculate with a numerical integration. But as we wrote, when minimum of error function is zero, the algorithm finds always a solution and this is not understandable.

If you comment the measures in "CDSGE 3 Layer WTGA2" (row 83) and use the misures in "Test Perfect 3 Layer" (row 111) you can see that the method works well.

Please, let us know if is there an alternative way to proceed?

What do you suggest?

[Gianluca_G_1]

Hello,

In these days we have analyzed the problem more deeply, and we arrived at these conclusions:

1. The method not always converges, but generally yes. Good starting values are suggested.

For this reason it is neccessary, to insert a timeout control inside the main loop

while (!bSuccessful)

{ ...         

2. The mathematical minimums not always have phisical meanings.

For this reason good lower and upper bounds are essentials, in order to filter the results. We didn't use good values in the example we provided to you.

3. It is better to iterate more times the algorith, using the previous step results as input of each new iteration.

Now we are able to compare our method with the results of our competitors.

BR

Gianluca
 

[mecej4]

You have rediscovered well-known properties of nonlinear least-squares algorithms. I hinted at one or two of them in #4. Here is another tip: it is often useful to use more than one algorithm on a single problem, separately or in sequence, because one algorithm may succeed where the other fails, or one algorithm may yield a better minimum when started with the final result of the other.