Gauss-Newton algorithm implementation

I am trying to implement a Gauss-Newton algorithm for non-linear optimisation. Unfortunately despite searching through the library and the internet I can't figure out quite how to use it in my case. As such I am hoping that I have come to the right place to ask for some assistance? Thank you to anyone willing to help me.

My problem consists of minimising a function $f(\vec x)$ where $\mathbb R^n \ni \vec x = (x_1, \dots ,x_n)$ and $f:\mathbb R^n\rightarrow \mathbb R$. Physically $f$ is a chi squared value $f=y-e(\vec x)$ where $y$ is the target value and $e(\vec x)$ is the test function evaluated at guess $\vec x$.

I begin by computing $f(\vec x+\delta \vec x)$ and $f(\vec x)$. Here $\vec x+\delta \vec x=x_1+\delta x_1,\dots , x_n+\delta x_n$. I then evaluate numerical derivatives using the following quotient formula, $$\frac{\partial f(\vec x)}{\partial x_i} = \frac{f(x_1,\dots , x_i+\delta x_i,\dots ,x_n)-f(\vec x)}{\delta x_i}\qquad \forall i =1,\dots ,n$$ I must then combine this into the Jacobian $\vec J$, $$\vec J = \bigg(\frac{\partial f(\vec x)}{\partial x_1}\ \frac{\partial f(\vec x)}{\partial x_2} \cdots \frac{\partial f(\vec x)}{\partial x_n}\bigg)$$ I can take the transpose of $\vec J$ denoted $\vec J^T$.

So far I can code the above (in C++). However The next part becomes confusing for me personally. In particular I do not know how to write the following matrix equation. $$(\vec J \vec J^T)\delta \vec x = \vec Jf$$ In fact I am not certain that this is the correct procedure?

Therefore, with numerical derivatives in place how can I obtain $\delta \vec x$ for the next iteration of the descent? Again many thanks for your time.

• More suitable for scicomp.stackexchange.com ? – Anton Trunov Mar 22 '16 at 21:21
• @AntonTrunov I think so, yes. The question seems to be principally about implementation, which is off-topic, here. – David Richerby Mar 22 '16 at 23:20
• I think this is fine and on-topic here. We don't migrate questions that are on-topic here (even if it might be a slightly better fit on another site, we don't migrate it -- the poster chose where they want to post it, and if it's on-topic here, we don't migrate it against their will -- we only consider migrating if it is off-topic here). So it's fine for this to stay here. – D.W. Mar 23 '16 at 4:54
• Maybe this would be on-topic and well received for Mathematics as well. – Juho Apr 5 '16 at 17:52
• Do you really mean Gauss-Newton? Gauss-Newton is for minimizing sums of squares, but I can't see any evidence that your $f(x)$ has the form of sum of squares. – D.W. Jul 3 '17 at 15:03

Solve the equation multiplying both terms of the equation for the inverse of the hessian matrix $\vec{J} \vec{J}^T$ on the left side:
$$\delta \vec{x} = (\vec{J} \vec{J}^T)^{-1} \cdot \vec{J}f$$