From Introduction to Local Minimization:
With Method -> Automatic, the Wolfram Language uses the quasi-Newton method unless the problem is structurally a sum of squares, in which case the Levenberg-Marquardt variant of the Gauss-Newton method is used.
To confirm, I use the first example in FindFit >> Options >> Method and compare the output for various settings for the Method option:
Possible settings for Method include "ConjugateGradient", "Gradient", "LevenbergMarquardt", "Newton", "NMinimize", and "QuasiNewton", with the default being Automatic.
model = a Exp[-b (x - c)^2] + d Sin[ω x + ϕ];
data = Table[{x, model /. {a -> 2, b -> 1, c -> 0, d -> 2, ω -> 0.67, ϕ -> 0.1}},
{x, -5, 5, .1}] + RandomReal[.25, 101];
methods = {Automatic, "ConjugateGradient", "Gradient", "LevenbergMarquardt",
"Newton", "NMinimize", "QuasiNewton"};
Grid[Partition[Labeled[Quiet@NonlinearModelFit[data, model, {a, b, c, d, ω, ϕ},
x, Method -> #]["ParameterTable"], #, Top] & /@ methods, 3], Spacings -> {5, 5}]
