4
$\begingroup$

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.

$\endgroup$
10
  • 2
    $\begingroup$ If someone asked me for an estimate about the quality of your post, I would say it's in the best 5% of first posts here. Therefore, welcome here. Please skim over the FAQ and read how voting and accepting of answers works here. (I'll remove the last sentence of your post. Really not needed) $\endgroup$
    – halirutan
    Feb 13, 2015 at 18:50
  • $\begingroup$ What is u? It is not defined. $\endgroup$
    – halirutan
    Feb 13, 2015 at 19:04
  • $\begingroup$ u cant be a numeric argument to P and the integration variable inside..(?). $\endgroup$
    – george2079
    Feb 13, 2015 at 19:30
  • $\begingroup$ @george2079 Ahh, sorry, haven't seen it in the integrate call. That makes even less sense. $\endgroup$
    – halirutan
    Feb 13, 2015 at 19:48
  • $\begingroup$ MaxIterations->1 makes no sense, it basically says do one eval and return the initial guess. $\endgroup$
    – george2079
    Feb 13, 2015 at 20:05

1 Answer 1

4
$\begingroup$

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.

$\endgroup$
7
  • $\begingroup$ This was extremely useful actually. Thank you. Although, you are right, this does not really make the program work, but it does mean that I need to fundamentally rethink how quickly I either evaluate this function, or how I perform the fitting. Not sure which yet though. $\endgroup$ Feb 14, 2015 at 3:26
  • $\begingroup$ @GregoryFahs Please see my edit. I approximated your integral with a stupid sum and it doesn't seem to work that bad. Maybe you can use this as an initial guess for a more expensive method. $\endgroup$
    – halirutan
    Feb 14, 2015 at 4:01
  • $\begingroup$ I agree with you, reducing the number of sampling points could help. The number of points I am using is actually defined by the detector we use, but that number is fairly arbitrary as far as I am concerned, really I only need to use a number that follow 2^n for the sampling number. Additionally, this fit now only takes 15 minutes even if I include 1000 points, I think the key was to make the integral run faster by compiling it as you did. This gives me a lot to play around with, I really appreciate all of your help. $\endgroup$ Feb 14, 2015 at 6:18
  • 1
    $\begingroup$ Just for completeness I figured I would add in what this program actually does. It is computing theoretical data for small angle x-ray scattering of dilute solutions of ellipsoids. Essentially the function FHali as halirutan has defined it is the form factor for ellipsoids of rotation (This is a Fourier transform of all orientations of an ellipsoid). Here is a link to the equation I am fitting: ncnr.nist.gov/resources/sansmodels/Ellipsoid.html $\endgroup$ Feb 14, 2015 at 6:22
  • $\begingroup$ NMinimize I think will use differential evolution by default, which (being a population-based method) is notoriously profligate with its use of objective function evaluations per iteration. That's why I specified Nelder-Mead in my comment above. $\endgroup$ Feb 14, 2015 at 14:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.