I am working on a non-linear least squares problem with standard form, in which I need to calibrate a parameter vector $\Theta$ to a set of inputs $\mathbf{x}$ and outputs $\mathbf{y}$:

$$\begin{align} & \min_\Theta ||R(\mathbf{x},\mathbf{y};\Theta)||_2^2 \\[4pt] &\quad \text{st: } R(\mathbf{x},\mathbf{y};\Theta):=\mathbf{y}-f(\mathbf{x};\Theta) \end{align}$$

More specifically, I want to fit the 2-dimensional vector $\Theta=(\theta_1,\theta_2)$ to $\mathbf{y}=\{y_1,y_2\}$ based on inputs $\mathbf{x}=\{x_1,x_2\}$. To do so, I have implemented a Gauss-Newton algorithm. The function $f$ is non-linear and requires Monte Carlo methods to evaluate (I am fixing a given simulation seed to avoid noise during the fitting procedure).

I have started by writing a 1-dimensional Newton algorithm for benchmarking purposes. I have then tried to fit $\theta\in\{\theta_1,\theta_2\}$ to one point $y\in\{y_1,y_2\}$ $-$ so yielding four fitting programs in total. All four optimizations above have worked, meaning satisfactory convergence in few iterations.

Next, I have written a Gauss-Newton routine to fit $\Theta$ to $\mathbf{y}$, and this is where my problems start. My Gauss-Newton update at the first iteration is way too large compared to what I was obtaining when fitting individual parameters. For example:

  • When fitting $\theta_1$ to $y_1$ using Newton, on the first iteration I was obtaining a Newton update approximately equal to $0.002$. However, under the Gauss-Newton algorithm, I get an update 3 orders of magnitude larger for $\theta_1$, i.e. approximately $-4.268$.
  • Similarly, when fitting $\theta_2$ to $y_2$ using Newton, on the first iteration of the Newton algorithm I get an update of approximately $2.6$, whereas under Gauss-Newton I get again 3 orders of magnitude more, $1616.5$.

I've done the following checks:

  • My Gauss-Newton algorithm yields the same results than my Newton one when applied to individual parameter/data-point fitting.
  • The four entries on the Jacobian matrix from the first iteration of Gauss-Newton, coincide with the derivatives I was obtaining on the first iteration of each one of the four single parameter-target problems above.
  • Using the solutions from the individual Newton fits as initializers for $\Theta^{(0)}$.

Below is the Jacobian I am getting on the first iteration to 7 decimals, as well as the initial state of the squared residuals $R^2_0$: $$ J_{\text{iter}_1}=\begin{pmatrix} 4.998\times 10^{-4} & 1.320\times 10^{-6} \\ 1.520\times 10^{-3} & 4.021\times 10^{-5} \end{pmatrix} \qquad R_0^2 =\begin{pmatrix} 1.189\times 10^{-6}\\ 1.048\times 10^{-6} \end{pmatrix} $$ such that the Gauss-Newton update is: $$\Theta^{(1)} = \Theta^{(0)}-J_{\text{iter}_1}^{-1}R_0^2$$

I am unsure what the problem is here. Starts to look like there is a fundamental issue with applying the Gauss-Newton method, but for the problem I am tackling this algorithm should be applicable.

Any pointers as to what the problem might be or what to look for to continue debugging are appreciated.

  • $\begingroup$ Why is the question being downvoted? $\endgroup$ – Daneel Olivaw Apr 1 at 16:31
  • 1
    $\begingroup$ The Gauss-Newton update would normally be $\Theta^{1}=\Theta^{0}-({J^{0}}^{T}J^{0})^{-1}{J^{0}}^{T}R^{0}$. $\endgroup$ – Brian Borchers Apr 1 at 16:54

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.