Return to Revisions

You mentioned in a related question that you would like to do large numbers of these fits quickly.

If each data set has the same dimensions, you can write a fairly fast implementation of Rahul Narain's method by precomputing arrays of x and y coordinates for the data grid, and flattening the data and using Dot to calculate the mean and the elements of the covariance matrix:

x = Table[i, {i, 9}, {j, 9}]//N;
y = Transpose[x];
x = Flatten[x]; y = Flatten[y];

semiaxes[data_] := Module[{min, p, mx, my},
  min = Min[data];
  p = Flatten[data] - min;
  p /= Total[p];
  mx = x.p;
  my = y.p;
  With[{a = (x - mx)^2.p, b = ((x - mx) (y - my)).p, c = (y - my)^2.p},
   Sqrt @ Eigenvalues[{{a, b}, {b, c}}]]]

semiaxes[data]
(* {1.86325, 1.50567} *)

This runs in about 340 µs on my PC

Compiling can give you even more speed, but you need to replace Eigenvalues with the explicit symbolic expressions:

semiaxesc = 
  With[{x = x, y = y}, 
   Compile[{{data, _Real, 2}}, Block[{min, p, mx, my, a, b, c},
     min = Min[data];
     p = Flatten[data] - min;
     p /= Total[p];
     mx = x.p;
     my = y.p;
     a = (x - mx)^2.p;
     b = ((x - mx) (y - my)).p;
     c = (y - my)^2.p;
      {Sqrt[1/2 (a + c - Sqrt[a^2 + 4 b^2 - 2 a c + c^2])],
       Sqrt[1/2 (a + c + Sqrt[a^2 + 4 b^2 - 2 a c + c^2])]}], 
    CompilationTarget -> "C", RuntimeOptions -> "Speed"]];

semiaxesc[data]
(* {1.50567, 1.86325} *)

This runs in about 5.7 µs on my PC.