The eventual route I followed was to adapt a version of the NM-simplex written in C and link to it with MathLink. Once I've tested this out fully I'll post it here.

I've been playing around with Oleksandr's compiled Nelder-Mead Minimization function presented here: Shaving the last 50 ms off NMinimize

First, I'll create some sample data. This creates a noisy version of a Gaussian function.

gf1[x_] := 1/Sqrt[2 Pi 5^2] Exp[-((x - 25)^2/(2 5^2))];
addSomeGaussianNoise[x_, sig_] := Map[# + RandomVariate[
    NormalDistribution[0., sig]] &, x, {1}];

interval = 0.01;
values = Range[0, 50, interval];
cleandata = gf1[#] & /@ values;
data = addSomeGaussianNoise[cleandata, 0.01];

enter image description here

There are a few ways to minimise the model. I ran NelderMeadMinimize from the linked question:

fitter = Block[{data = #}, NelderMeadMinimize[Block[{x = values}, 
      Norm[data - 1/Sqrt[2 Pi c^2] Exp[-((x - b)^2/(2 c^2))]]],
      {b,c}]] &;
fitter2 = Block[{data = #}, NelderMeadMinimize[Block[{x = values}, 
      EuclideanDistance[data, 1/Sqrt[2 Pi c^2] Exp[-((x - b)^2/(2 c^2))]]],
      {b, c}]] &;
fitter3 = Block[{data = #}, NelderMeadMinimize[
        Abs[data[[x/interval + 1]] - 1/Sqrt[2 Pi c^2] Exp[-((x - b)^2/(2 c^2))]]^2, 
      {b, c}]] &;
fitter4 = With[{inter = interval, x = Range@Length[#]},
        Abs[#[[x]] - 1/Sqrt[2 Pi c^2] Exp[-(((x - 1)*inter - b)^2/(2 c^2))]]^2]], 
      {b,c}]] &;

fits = fitter /@ {data} // AbsoluteTiming
(*{0.299809, {{0.712906, {b -> 24.985, c -> 5.0084}}}}*)

fits2 = fitter2 /@ {data} // AbsoluteTiming
(*{0.240542, {{0.712906, {b -> 24.985, c -> 5.0084}}}}*)

fits3 = fitter3 /@ {data} // AbsoluteTiming
(*{6.575716, {{0.712906, {b -> 24.985, c -> 5.0084}}}}*)

fits4 = fitter4 /@ {data} // AbsoluteTiming
(*{0.291925, {{0.712906, {b -> 24.985, c -> 5.0084}}}}*)

They all return the same result, but in varying times. Interestingly EuclideanDistance is slightly faster than Norm.

Option 3 isn't well-written - Option 4 is better. Both 3 and 4 are the most obvious ones for me to choose to extend to a 2D (multivariate) Gaussian function, as I understand them better.

Meanwhile, the first 2 options using Norm and EuclideanDistance are good, but don't allow for a different minimization function such as that in a question I asked recently. In that scenario I was happy using NMinimize, but now I'm looking to apply this alternative Nelder-Mead method to see what performance gains can be had.

My question is therefore:

How can I get any of the options to work with 2D data rather than 1D data? i.e. fitting a Gaussian function in $x$ and $y$.

Option 3/4 can be improved further, such that the fitting only takes 0.04 seconds, namely:

   epsilon = $MachineEpsilon,
   minimizer = NelderMeadMinimize`Dump`CompiledNelderMead[
     Function[{b, c},
      Block[{y = Range@Length[data]},
         Abs[data[[y]] - 1/Sqrt[2 Pi c^2] Exp[-(((y - 1)*interval - b)^2/(2 c^2))]]^2]
        ]]],{b, c},"ReturnValues" -> "AugmentedOptimizedParameters"]
  serialFitter = 
    Compile[{{data, _Real, 1},{interval,_Real}}, 
     minimizer[RandomReal[{0, 1}, {2 + 1, 2}], epsilon, -1], 
     CompilationOptions -> {"InlineCompiledFunctions" -> True}, 
     RuntimeOptions -> {"Speed", "EvaluateSymbolically" -> False}, 
     CompilationTarget -> "C"];

serialFitter[data,interval] // AbsoluteTiming
(* {0.046861, {0.711633, 25.0001, 5.02656}} *)

Update #1

This creates the 2D data

gf[x_, y_] := 1/Sqrt[2 Pi 5^2] Exp[-((x - 25)^2/(2 5^2) + (y - 25)^2/(2 5^2))];
addSomeGaussianNoise[x_, sig_] := 
       Map[# + RandomVariate[NormalDistribution[0., sig]] &, x, {2}];
interval = 1.;
cleandata = Table[gf[x, y], {x, 0, 50, interval}, {y, 0, 50, interval}];
data = addSomeGaussianNoise[cleandata, 0.01];

I then try and fit with this, where I've flattened the 2D data to a 1D vector prior to fitting.

fitter = Compile[{{data, _Real, 1}, {interval, _Real}},
      {b, c, d},
      Block[{x = Range@Length[data], n = Sqrt[Length[data]]},
           data[[x]] - 
            1/Sqrt[2 Pi c^2] Exp[-((Floor[Mod[x - 1, n]]*interval - b)^2/(2 c^2) +  
              (Quotient[x - 1, n]*interval - d)^2/(2 c^2))]]^2]
     {b, c, d}, "ReturnValues" -> "AugmentedOptimizedParameters"]
    [RandomReal[{0, 1}, {3 + 1, 3}], $MachineEpsilon, -1],
   CompilationOptions -> {"InlineCompiledFunctions" -> True}, 
   RuntimeOptions -> {"Speed", "EvaluateSymbolically" -> False},
   CompilationTarget -> "C"];

fitter[Flatten[data], 1.] // AbsoluteTiming

However, I get errors:

CompiledFunction::cfte: Compiled expression {0.874982,0.207452,-0.119287,0.635137} should be a rank 2 tensor of machine-size real numbers.

CompiledFunction::cfexe: Could not complete external evaluation; proceeding with uncompiled evaluation.

It's nearly there as a solution to my problem, but not quite!

Update #2

The first error (cfte) refers to the "ReturnValues" option - if I change it from "AugmentedOptimizedParameters" to "OptimizedParameters" the error now says

CompiledFunction::cfte: "Compiled expression {-0.053765,0.291736,0.432465} should be a rank 2 tensor of machine-size real numbers.

See how it's changed length? The function is expecting the return values to be a rank 2 tensor, but they're not...

  • 1
    $\begingroup$ Perhaps Oleksandr is the best-placed person to answer this (it's his code after all), but any help appreciated! $\endgroup$ Commented Jul 10, 2014 at 20:39
  • 1
    $\begingroup$ Could you please clarify what you mean by two-dimensional? The sum of two of these Gaussian curves as in the linked question?--or is it that rather than a univariate Gaussian, you want to find the best fitting multivariate Gaussian? NM minimizes a function $f(\vec{v}) : \mathbb{R}^n \rightarrow \mathbb{R}^1$ with respect to $\vec{v}$, so as long as you can cast your problem in this form, you can do it. $\endgroup$ Commented Jul 13, 2014 at 13:40
  • 1
    $\begingroup$ Also, NelderMeadMinimize compiles whatever is passed to it, so this should be a self-contained problem (without any references to external symbols) if calls out of the VM are to be avoided. Otherwise you have to use NelderMeadMinimize`Dump`CompiledNelderMead to produce the code for the minimizer and write the surrounding code yourself so that all of the necessary symbols are defined within the VM on entry. Failure to do so will still produce a result, but the many calls out of the VM present a performance problem. $\endgroup$ Commented Jul 13, 2014 at 13:51
  • 1
    $\begingroup$ Your question seems like it should be reasonably straightforward to address, but I have things I need to do until later this evening. Can I get back to you in a few hours? $\endgroup$ Commented Jul 13, 2014 at 14:15
  • 1
    $\begingroup$ Sorry about the delay in replying--I had been playing around with your chi-square goodness of fit statistic from the other post and found it a bit difficult to work with. But that's not a fundamental limitation and I'll write an answer later. Basically the problem can be reduced to changing your residual matrix into a residual vector and treating it as usual; the mathematics is the same either way, so you can just Flatten it if you want... $\endgroup$ Commented Jul 14, 2014 at 11:18


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.