Problems extending NelderMeadMinimize into 2D

Question

Update

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];

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[
      Sqrt[Sum[
        Abs[data[[x/interval + 1]] - 1/Sqrt[2 Pi c^2] Exp[-((x - b)^2/(2 c^2))]]^2, 
         {x,values}
      ]],
      {b, c}]] &;
fitter4 = With[{inter = interval, x = Range@Length[#]},
     NelderMeadMinimize[
      Sqrt[Total[
        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:

With[{
   epsilon = $MachineEpsilon,
   minimizer = NelderMeadMinimize`Dump`CompiledNelderMead[
     Function[{b, c},
      Block[{y = Range@Length[data]},
       Sqrt[Total[
         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}},
   NelderMeadMinimize`Dump`CompiledNelderMead[
     Function[
      {b, c, d},
      Block[{x = Range@Length[data], n = Sqrt[Length[data]]},
       Sqrt[
        Total[
         Abs[
           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...