When trying to solve for a (over)determined non-linear least square method:

$$\underset{x}{\min}||f(x)||^2_2, f: \mathbb{R}^n \rightarrow \mathbb{R}^m, x\in \mathbb{R}^n, m\geq n$$

we use the Gauss-Newton algorithm to solve it iteratively.

At each iteration $k\in\mathbb{N}$, the system is linearized at the current iteration of the independent variables $x_k$, using the Taylor expansion and discarding the higher order terms:

$$f(x)\approx f(x_k)+J(x_k)(x-x_k)$$

where $J(x_k)\in \mathbb{R}^{m \times n}$ is the jacobian matrix of $f(x)$ with respect to $x$, evaluated at $x_k$. It becomes a linear least-square problem:

$$\underset{x}{\min}||f(x_k)+J(x_k)(x-x_k)||^2_2 = \underset{x}{\min} (f-Jx_k)^T(f-Jx_k) + 2x^TJ^T(f-Jx_k) + x^TJ^TJx$$

The right-hand side is after expansion and the calculation of the 2-norm. Sorry for the abuse of notation here by dropping the $(x_k)$. This minimization has a closed-form solution by taking the derivative with respect to $x$ and setting the result to be $0$:

$$J^TJx+J^T(f-Jx_k)\stackrel{!}{=} 0$$

In the case where $J$ is full column rank and $J^TJ$ is invertible, the result would be:

$$x_{k+1} = -(J^TJ)^{-1}J^T(f-Jx_k)=x_k -(J^TJ)^{-1}J^Tf$$

This is the Gauss-Newton algorithm. In the case where $J$ is not full column rank and $J^TJ$ becomes non-invertible, the Levenberg–Marquardt algorithm proposes adding a regularization term to the objective:

$$\underset{x}{\min}||f(x_k)+J(x_k)(x-x_k)||^2_2 + \lambda||x-x_k||^2_2$$

The solution to this minimization becomes:

$x_{k+1} = x_k -(J^TJ+\lambda I)^{-1}J^Tf$, where $(J^TJ+\lambda I)$ becomes invertible.

My question is, after such a long bedding: at step $k$, after the linearization:

$$f(x)\approx f(x_k)+J(x_k)(x-x_k)=J(x_k)x + f(x_k) - J(x_k)x_k$$

Can we directly equate this function to $0$ and solve for $x$ in the least-square sense by taking the pseudo-inverse of $J(x_k)$?

$$J(x_k)x + f(x_k) - J(x_k)x_k \stackrel{!}{=} 0$$

$$x_\mathrm{Least square} = J^{\dagger}(x_k)J(x_k)x_k-J^{\dagger}(x_k) f(x_k)$$

Doing this does not require $J(x_k)$ to be full column rank. It does not even need the system to be over(determined) because the pseudo-inverse can be found anyway using the SVD decomposition. It becomes the Gauss-Newton algorithm when $J$ is full column rank. Would it be a more general solution to tackle rank deficiency in the Gauss-Newton algorithm?

  • $\begingroup$ You aren't trying to find a zero of $f(x)$, but rather a minimum of $f(x)$ where $\nabla f(x)=0$. $\endgroup$ Commented Jun 4 at 23:12
  • $\begingroup$ Isn't $f(x)$ the error vector in the least square problems, e.g. in the line fitting problem, the error between the estimated line and data points? Do we not want to equate that vector to the zero vector and solve for $x$? $\endgroup$ Commented Jun 5 at 20:44
  • 1
    $\begingroup$ OK- you've used $f(x)$ for the entire vector of residuals. You want to minimize $\| f(x) \|^{2}$ by setting the gradient of $\| f(x) \|^{2}=0$, not by setting $f(x)=0$. $\endgroup$ Commented Jun 5 at 23:37
  • $\begingroup$ The gradient of $\| f(x) \|^{2}$ is $J(x)^{T}f(x)$. $\endgroup$ Commented Jun 5 at 23:49
  • $\begingroup$ Doesn't setting $f(x)=0$ and solving for $x$ using pseudo-inverse lead to the least square solution? In the case that $J(x_k)$ is rank-deficient, setting the gradient of the norm to $0$ would not lead to a closed-form solution, right? $\endgroup$ Commented Jun 6 at 14:01


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.