How to compute the Hessian matrix in the robustified Gauss-Newton method for optimization

Question

When robustifying the Gauss-Newton method used to find a local minima of a squared error function expressed by: $$f(\mathbf{x})=\frac{1}{2}\Delta\mathbf{z}(\mathbf{x})^T\Delta\mathbf{z}(\mathbf{x})$$, as explained in works like Bundle adjustment: A Modern Synthesis (page 17), we have to introduce a strong nonlinearity into this cost function, which is denoted by $\rho(\cdot)$ (it can be e.g. the Huber loss).

Thereby to reduce the influence of outliers, we are no longer interested in minimizing $f(\mathbf{x})$, but the robustified version of it: $\rho(f(\mathbf{x}))$. Hence, the gradient vector ($\mathbf{g}$) and the approximate Hessian matrix ($\mathbf{H}$) are given by:

$$ \mathbf{g} = \rho'\mathbf{J}^T\Delta\mathbf{z},\quad\quad \mathbf{H} \approx \mathbf{J}^T(\rho'\,\mathbf{I} + 2\rho''\Delta\mathbf{z} (\Delta\mathbf{z})^T)\mathbf{J} $$ Where $\rho'$ and $\rho''$ are the first and second derivatives of $\rho(f(\mathbf{x}))$ w.r.t $f(\mathbf{x})$. And $\mathbf{J}$ is the jacobian matrix given by: $\mathbf{J} = \partial \Delta\mathbf{z}/\partial \mathbf{x}$.

---

Given this, what I don't understand is how this approximation of the Hessian matrix is computed.

My attempt to reach it, is the following: $$\begin{align} \mathbf{H} &= \frac{\partial^2\rho(f(\mathbf{x}))}{\partial \mathbf{x}^2} = \frac{\partial\rho(f(\mathbf{x}))}{\partial f(\mathbf{x})}\frac{\partial^2 f(\mathbf{x})}{\partial \mathbf{x}^2} + \frac{\partial^2\rho(f(\mathbf{x}))}{\partial f(\mathbf{x})^2}\left(\frac{\partial f(\mathbf{x})}{\partial \mathbf{x}}\right)^2 && \text{By applying chain rule} \end{align} $$ Where: $$\begin{align} &\frac{\partial\rho(f(\mathbf{x}))}{\partial f(\mathbf{x})}\frac{\partial^2 f(\mathbf{x})}{\partial \mathbf{x}^2} \approx \rho'\mathbf{J}^T\mathbf{J}\quad\leftrightarrow\quad\text{Gauss-Newton approximation}\\ \\ &\frac{\partial^2\rho(f(\mathbf{x}))}{\partial f(\mathbf{x})^2}\left(\frac{\partial f(\mathbf{x})}{\partial \mathbf{x}}\right)^2 = \rho''\,\mathbf{J}^T\Delta\mathbf{z} (\Delta\mathbf{z})^T\mathbf{J} \end{align} $$ Thereby I am getting an approximated Hessian where the second term is not multplied by 2: $$ \mathbf{H}_{\text{attempt}} \approx \mathbf{J}^T(\rho'\,\mathbf{I} + \rho''\Delta\mathbf{z} (\Delta\mathbf{z})^T)\mathbf{J} \neq \mathbf{J}^T(\rho'\,\mathbf{I} + \color{red}{2}\rho''\Delta\mathbf{z} (\Delta\mathbf{z})^T)\mathbf{J} $$ Can please someone help me to see where I am doing the wrong step?

Thanks in advance!

Answer 1

I was reading the same paper today. If you look at equation (3), the objective function is actually $$\frac{1}{2}\rho(\mathbf{z}(\mathbf{x})^T \mathbf{z}(\mathbf{x})),$$ not what you stated at first $$\rho(\frac{1}{2} \mathbf{z}(\mathbf{x})^T \mathbf{z}(\mathbf{x})).$$ Then $\left(\frac{\partial \mathbf{z}(\mathbf{x})^T \mathbf{z}(\mathbf{x})}{\partial \mathbf{x}}\right)^2$ will generate a coefficient $2^2=4$ at your second term. So nothing in your derivation is wrong.

Answer 2

Have a look at A.129 and A.158 of my thesis. This was my conclusion after reading multiple articles on the approximation method but I can't find the source articles.