From [Introduction to Local Minimization](http://reference.wolfram.com/language/tutorial/UnconstrainedOptimizationIntroductionLocalMinimization.html#509267359): 
> **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](http://reference.wolfram.com/mathematica/ref/FindFit.html) 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}]

![enter image description here][1]


  [1]: https://i.sstatic.net/sOVW5.png