BFGS in comparison with Gaussian newton I'm reading about quasi-newton method, and I get the key idea is to find an approximated Hessian Matrix. The most popular one is for sure L-BFGS. But I've got see Gaussian newton method simply estimate the Hessian with $\nabla J \cdot \nabla J^T$. Isn't it much more simpler than BFGS update? why it is much less used?
 A: Consider the residual: $r_n(\beta) = \|y_n - f(x_n;\beta)\|_2$. We generally want to minimise the sum of square residuals: 
\begin{equation}
S(\beta) = \sum_{i = 1}^N  \|y_n - f(x_i;\beta)\|_2^2 
\end{equation}
If we differentiate this then notice,
\begin{align}
\frac{\partial S(\beta)}{\partial \beta_i} &= 2 \sum_{n = 1}^N  \|y_n - f(x_n;\beta)\|_2\cdot \frac{\partial \|y_n - f(x_n;\beta)\|_2}{\partial \beta_i} \\
&= 2 \sum_{n = 1}^N  r_n \cdot \frac{\partial r_n}{\partial \beta_i}
\end{align}
Let's differentiate once more:
\begin{align}
\frac{\partial^2 S(\beta)}{\partial \beta_i\partial \beta_j} &= 2 \sum_{n = 1}^N \left( \frac{\partial r_n}{\partial \beta_i} \cdot \frac{\partial r_n}{\partial \beta_j} + r_n \frac{\partial^2 r_n}{\partial \beta_i \partial \beta_j}\right)
\end{align}
Now if we set $\frac{\partial^2 r_n}{\partial \beta_i \partial \beta_j} = 0$, notice that we are left with just the multiplication of two Jacobian terms $ \frac{\partial r_n}{\partial \beta_i} \cdot \frac{\partial r_n}{\partial \beta_j} = J^T \cdot J$, as your question asks. 
Therefore, it is an approximation which ignores any cross-coupling terms in the Hessian. This is good for systems when we expect a large degree of independence between terms (very little cross coupling).
This is Gauss-Newton's method with an approximation on the Hessian, which naturally arises from first principles, by differentiating the cost function. 
Now, methods like BFGS, are quasi-Newton methods. Quasi-Newton methods also try to avoid using the Hessian directly, but instead they work to approx. it by (in the case of BFGS), progressively updating an approx. per iteration. This is shown below, where B is approximate Hessian (taken from wiki)
\begin{equation}
B_{k+1}=B_{k}+{\frac {\mathbf {y} _{k}\mathbf {y} _{k}^{\mathrm {T} }}{\mathbf {y} _{k}^{\mathrm {T} }\mathbf {s} _{k}}}-{\frac {B_{k}\mathbf {s} _{k}\mathbf {s} _{k}^{\mathrm {T} }B_{k}^{\mathrm {T} }}{\mathbf {s} _{k}^{\mathrm {T} }B_{k}\mathbf {s} _{k}}}.
\end{equation}
As for when to use which when, it is hard to say. It is problem specific. I would look into literature your topic and see what they say. Since these are all approximation methods, the underlying objective function is not guaranteed to be traversed accurately in either case. For example, in general LGBS type solutions are popular for non-linear least squares problems, because they work well historically, and have nice and optimized libraries (as far as I understand). 
Why is that Gauss-Newton update much less used? Well as we discussed before it places a large assumption on the objective space (negligible cross terms). It can be detrimental for certain applications. That's why I recommend looking at the literature for your topic in particular. 
