3
$\begingroup$

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<n$ (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?

$\endgroup$
0

1 Answer 1

2
$\begingroup$

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).$$

$\endgroup$
7
  • $\begingroup$ Thanks a lot for the answer (and sorry for the late reply-had to refresh my linear algebra skills at first). This explanation makes sense for me. As I understand it, we assume that the objective function g(x) goes to zero, since our system is underdetermined. The optimization step is then calculated based on the linearization around x. Therefore, the derivation of the method is different to what I did above in my question. Do you think I could also similarly derive this form from Newton's method? $\endgroup$
    – ubongoUga
    Feb 21, 2022 at 9:15
  • $\begingroup$ What I still don't quite understand is the decomposition into Delta_x=y+z. Could you point me to a textbook or web article where I can read more about this? $\endgroup$
    – ubongoUga
    Feb 21, 2022 at 10:37
  • $\begingroup$ @ubongoUga Yes, this is the critical point The key word is orthogonal complement. Your favorite linear algebra textbook will discuss it. Any vector can be split into two components, one is our vector space of choice and the other is in the orthogonal complement of the vector space. $\endgroup$ Feb 21, 2022 at 11:01
  • $\begingroup$ Cool. Thank you again for your help, this was really helpful. $\endgroup$
    – ubongoUga
    Feb 21, 2022 at 11:06
  • $\begingroup$ @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. $\endgroup$ Feb 21, 2022 at 15:47

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .