# Derivation of Gauss-Newton method for underdetermined inverse problems

I would like to derive the Gauss-Newton method for solving underdetermined inverse problems, i.e., where the parameter space is larger than the number of constraints. Although my final solution works, I do not understand how the method can be derived for such underdetermined problems:

Let $$f$$ be a scalar objective function $$f:\mathbb{R}^n \rightarrow \mathbb{R}$$ that we wish to minimize using the Gauss-Newton algorithm, i.e., $$f=\sum_{i=1}^m r_i^2(\underline{\alpha}),$$ where $$\underline{\alpha} \in \mathbb{R}^n$$ is a vector of $$n$$ parameters, $$r_i$$ describes the i-th constraint (residual) and $$m$$ is the number of constraints. To the best of my knowledge and according to Wikipedia , the Gauss-newton algorithm can be derived directly from Newton's method, i.e., $$\underline{\alpha}^{\nu+1} = \underline{\alpha}^{\nu} - \underline{\underline{H}}^{-1}\underline{g},\tag{1}\label{eq1}$$ where $$\underline{g}$$ and $$\underline{\underline{H}}$$ are the gradient vector and Hessian matrix of $$f$$. The gradient vector in index notation is given as $$g_j=2\sum_{i=1}^m r_i \frac{\partial r_i}{\partial \alpha_j}$$ and the Hessian is $$H_{jk}=2\sum_{i=1}^m\left(\frac{\partial r_i}{\partial \alpha_j} \frac{\partial r_i}{\partial \alpha_k} +r_i \frac{\partial^2 r_i}{\partial \alpha_j \partial \alpha_k}\right).$$ In order to obtain the Gauss-Newton method, we ignore the second-order derivatives of the Hessian. Thus, $$\underline{g} = 2 \underline{\underline{J}}_r^T \underline{r}$$ and $$\underline{\underline{H}} \approx 2 \underline{\underline{J}}_r^T \underline{\underline{J}}_r,$$ where $$\underline{\underline{J}}_r$$ is the Jacobian of the residual vector $$\underline{r}\in\mathbb{R}^m$$. With this, we can now write the Gauss-Newton algorithm, i.e., $$\underline{\alpha}^{\nu+1} = \underline{\alpha}^{\nu} - \left (\underline{\underline{J}}_r^T \underline{\underline{J}}_r \right)^{-1} \underline{\underline{J}}_r^T \underline{r}.\tag{2}\label{eq2}$$ For the usual case with $$m>n$$ (overdetermined problem), this algorithm works perfectly fine for me. However, if $$m (underdetermined inverse problem), $$\left (\underline{\underline{J}}_r^T \underline{\underline{J}}_r \right)$$ becomes singular and cannot be inverted. I realized that I can overcome this problem if I replace the left pseudoinverse of Equation \ref{eq2} by the right pseudoinverse, i.e., $$\underline{\alpha}^{\nu+1} = \underline{\alpha}^{\nu} - \underline{\underline{J}}_r^T\left (\underline{\underline{J}}_r \underline{\underline{J}}^T_r \right)^{-1} \underline{r}\tag{3}\label{eq3}.$$ With this modification, my algorithm works perfectly fine and converges to a solution. However, I do not understand which assumptions are necessary to derive Equation \ref{eq3}, e.g. similarly to what was done above for Equation \ref{eq2} based on Newton's method (\ref{eq1}). I assume that in my derivation, I implicitly assume $$m>n$$, but I do not know where.

Can anyone help me with this?

You have a system of equations $$g(x) = 0$$ where $$g : \mathbb{R}^m \rightarrow \mathbb{R}^n$$ and $$n \leq m$$ so that there are at least as many unknowns as there are constraints. By Taylor's formula we have $$g(x+ \Delta x) \approx g(x) + Dg(x) \Delta x$$ where $$Dg(x) \in \mathbb{R}^{n \times m}$$ is the Jacobian of $$g$$ at the point $$x \in \mathbb{R}^m$$. We observe that $$Dg(x)$$ is a wide matrix with at least as many columns as rows. In order to derive Newton's method we seek $$\Delta x \in \mathbb{R}^m$$ such that $$g(x) + Dg(x) \Delta x = 0.$$ Suppose that there exists such a $$\Delta x$$, then there are unique $$y \in \text{Ker}(Dg(x)), \quad z \in \text{Ran}(Dg(x)^T)$$ such that $$\Delta x = y + z.$$ Here $$\text{Ker}(A)$$ denotes the kernel or null space of the linear transformation $$A$$ and $$\text{Ran}(A)$$ denotes the range or column space of the linear transformation $$A$$. It is clear that $$y$$ is irrelevant since $$g(x) + Dg(x)\Delta x = g(x) + Dg(x) z.$$ Now if $$z \in \text{Ran}(Dg(x)^T)$$, then $$z = Dg(x)^T p$$ for some $$p \in \mathbb{R}^n$$. It follows that we should consider the equation $$g(x) + Dg(x) Dg(x)^T p = 0$$ If we assume that $$Dg(x)$$ has full row rank, then $$Dg(x) Dg(x)^T \in \mathbb{R}^{n \times n}$$ is nonsingular and our equation has a unique solution $$p$$. It follows Newton's method takes the form $$x \gets x - Dg(x)^T(Dg(x)Dg(x)^T)^{-1} g(x).$$
• @ubongoUga The two cases of $m < n$ and $m > n$ are really quite similar. I think that it would be worth you time to have a second look at $g(x) + Dg(x) \Delta x = 0$. Suppose that $Dg(x)$ is a tall matrix of full rank. Then this equation implies $Dg(x)^T g(x) + Dg(x)^T Dg(x) \Delta x = 0$ and we can solve for $\Delta x$. What we get is the least squares solution of the original overdetermined linear system. What I did above was to find the solution of the smallest 2-norm of an underdetermined system. Feb 21, 2022 at 15:47