FindFit runs extremely slow for my Function,how do I optimize this?

Question

My function is the first order Bessel Function, and I am basically just trying to get this to fit in a way that does not take forever.

I define the function as follows, (it seems to work faster if I explicity define the function instead of using the special function SphericalBesselJ):

F[q_?NumericQ, ra_?NumericQ, rb_?NumericQ, px_?NumericQ] :=px NIntegrate[((3 (4 Pi)/3 ra rb^2 1/(q rb (1 + x^2 ((ra/rb)^2 - 1)))^3 (Sin[q rb (1 + x^2 ((ra/rb)^2 - 1))] - (q rb (1 + x^2 ((ra/rb)^2 - 1)) Cos[q rb (1 + x^2 ((ra/rb)^2 - 1))])))^2),{x,0,1}]

In order to test the function I have simulated the function with random noise component that adds a 10% relative error in the data

TestData = Array[#1/1024*2.5 + 0*#2 &, {1024, 2}];

Do[
TestData[[n, 2]] = 
F[TestData[[n, 1]], 4.5, 6.8, 1]*RandomReal[{1 - .1, 1 + .1}]
, {n, 1024}]

In order to compute the fit I use a NormFunction that computes that uses the relative residuals instead of the absolute residuals as follows. The norm function I use is called weightedNorm[residuals].

weightedNorm[residuals_] := Norm[residuals/TestData[[All, 2]]]
IntensityWeightedParams = FindFit[TestData, {F[q, ra, rb,p],
                                 {0 < ra < 10, 5 < rb < 15}}, {{ra,4.5}, {rb, 6.8},{p,1}}, q,
                                 NormFunction -> weightedNorm, MaxIterations -> 10, Method -> NMinimize]

This fitting works particularly well when I am not solving a numerical integral, i.e. if I use a function that is purely algebraic, however, when I run this particular fit with even a single iteration as shown here the fit takes approximately 250s.

Answer 1

Not an answer at all but some insight. Let's first define two global variables we will use for looking at how fast your integral is evaluated:

calledIntegrate = 0;
evalStep = 0;

Now, let me redefine your target function by compiling the integrand to make it faster. Additionally, we will increase calledIntegrate on each call:

With[{
  cf = Compile[{{q, _Real, 0}, {ra, _Real, 0}, {rb, _Real, 0}, {x, _Real, 0}},
    ((3 (4 Pi)/3 ra rb^2 1/(q rb (1 + x^2 ((ra/rb)^2 - 1)))^3 (Sin[
           q rb (1 + x^2 ((ra/rb)^2 - 1))] - (q rb (1 + 
              x^2 ((ra/rb)^2 - 1)) Cos[q rb (1 + x^2 ((ra/rb)^2 - 1))])))^2)]
 },
 Module[{int},
  (* int is only a wrapper to force x to be numeric *)
  int[q_?NumericQ, ra_?NumericQ, rb_?NumericQ, x_?NumericQ] := cf[q, ra, rb, x];
  FHali[q_?NumericQ, ra_?NumericQ, rb_?NumericQ, px_?NumericQ] := 
    (calledIntegrate++; 
     px NIntegrate[int[q, ra, rb, x], {x, 0, 1}])
  ]
 ]

Now initialize you test data

TestData = Array[#1/1024*2.5 + 0*#2 &, {1024, 2}];
Do[
  TestData[[n, 2]] = FHali[TestData[[n, 1]], 4.5, 6.8, 1]*RandomReal[{1 - .1, 1 + .1}], 
{n, 1024}];

and now run the fit

Dynamic[{calledIntegrate, evalStep}]

weightedNorm[residuals_] := Norm[residuals/TestData[[All, 2]]]
IntensityWeightedParams = 
 FindFit[TestData, {FHali[q, ra, rb, p], {0 < ra < 10, 
    5 < rb < 15}}, {{ra, 4.5}, {rb, 6.8}, {p, 1}}, q, 
  NormFunction -> weightedNorm, MaxIterations -> 10, 
  Method -> "NMinimize", EvaluationMonitor :> evalStep++]

After some seconds, I have over 30000 calls to fHali but only 30 evaluation steps. I guess no matter how fast you can get your target function, it will never be fast enough so that you can use it with this approach.

Edit

One very crude idea is to skip NIntegrate. You will loose all the fancy stepping and adaption algorithms and I'm sure if the integrand is evil, the world will explode but maybe, it will give you an initial guess for your parameters. What I do is taking the same compiled code only that I'm simply dividing the interval [0,1] equally and replace the integration by a simple sum:

With[{cf = 
   Compile[{{q, _Real, 0}, {ra, _Real, 0}, {rb, _Real, 0}, {x, _Real, 
      0}}, ((3 (4 Pi)/
         3 ra rb^2 1/(q rb (1 + x^2 ((ra/rb)^2 - 1)))^3 (Sin[
           q rb (1 + x^2 ((ra/rb)^2 - 1))] - (q rb (1 + 
              x^2 ((ra/rb)^2 - 1)) Cos[
             q rb (1 + x^2 ((ra/rb)^2 - 1))])))^2), 
    Parallelization -> True, CompilationTarget -> "C", 
    RuntimeAttributes -> {Listable}],
  points = Table[x, {x, 0, 1, 1/100.}]},
 FHali[q_?NumericQ, ra_?NumericQ, rb_?NumericQ, px_?NumericQ] :=   
  px Plus @@ cf[q, ra, rb, points]/Length[points]
 ]

When I compare the execution time of 1000 calls with your original F (you need to define it again!)

AbsoluteTiming[Do[#[1.1, 4.5, 6.8, 1], {1000}];] & /@ {FHali, F}

it is 0.04s compared to 5.3s which is a factor of about 130. I'll take the TestData that was created with your proper function, but I use my crude approximation for the fit. Note that I removed the initial guesses for the parameters (and additionally changed the Method setting as suggested by Olek!)

steps = 0;
Dynamic[steps]

weightedNorm[residuals_] := Norm[residuals/TestData[[All, 2]]]
IntensityWeightedParams = 
 FindFit[TestData, {FHali[q, ra, rb, p], {0 < ra < 10, 
    5 < rb < 15}}, {ra, rb, p}, q, NormFunction -> weightedNorm, 
  Method -> {"NMinimize", Method -> "NelderMead"},
  EvaluationMonitor :> steps++]

After only 16s (!!) I got as answer

{ra -> 4.51598, rb -> 6.77539, p -> 0.993664}

which is not bad at all if we think about that we used just 100 sampling points. Maybe this helps you think about an alternative approach for your problem.