Gauss-Newton Method

For a nonlinear least-squares objective $\Phi (\sigma )=\frac{1}{2}\Vert r(\sigma ){\Vert }^2$ with residual $r(\sigma )=\mathcal{F}(\sigma )-y$ and Jacobian $J=r^{\prime }(\sigma )$:

• gradient: $\nabla \Phi =J^{\top }r$;
• Hessian: ${\nabla }^2\Phi =J^{\top }J+{\sum }_k{r}_k{\nabla }^2{r}_k$.

Gauss–Newton drops the second-order term, which is small near a good fit, and solves at each step the linearised least-squares problem

$$\underset{\delta }{\min }\Vert r+J\delta {\Vert }^2⟺J^{\top }J\delta =-J^{\top }r,$$

then updates $\sigma \leftarrow \sigma +\tau \delta$ with a step length $\tau$ from a Line Search.

Jacobian in EIT. Row $(i,j)$ of $J$ is the derivative of the $j$-th boundary measurement of pattern $i$. By the linearisation identity its entries are integrals of $-\nabla u_i\cdot \nabla v_j$, where $v_j$ is the solution driven by the $j$-th measurement functional. Building $J$ explicitly costs one extra solve per measurement. Matrix-free variants only need products $J\delta$ (linearised forward) and $J^{\top }w$ (adjoint solves), and solve the normal equations with CG or LSQR. With a factorisation of the forward problem, each product costs one back-substitution per pattern. A dense Gauss–Newton step costs $O(n^3)$ and $O(mn)$ memory; an inexact CG step costs a few dozen products and no dense matrix at all. On an EIT problem with $16384$ pixels, $255$ electrodes and $16384$ residuals, seven Levenberg–Marquardt iterations took 3.2 hours with dense steps and 104 seconds with CG steps, reaching the same reconstruction. The diagonal of $J^{\top }J$ (for Jacobi preconditioning and sensitivity damping) can be accumulated from row blocks of $J$ without storing it. So can the full Gram matrix, $J^{\top }J={\sum }_b{B}_b^{\top }{B}_b$ over row blocks $B_b$, when $n^2$ fits into memory but $mn$ does not. Each block is a symmetric rank-$b$ update. It should go into one triangle only (BLAS syrk), with the triangle mirrored once at the end. Mirroring after every update is a transpose copy of all $n^2$ entries, which is limited by memory bandwidth rather than arithmetic. For $n=16384$ and blocks of 64 rows, the copy took 5 to 10 seconds per block and the update 0.07 seconds.

If only per-pattern gradients $\nabla J_i$ are available (one adjoint solve per pattern, see Adjoint State Method), a cheaper reduced Gauss–Newton model treats each pattern’s misfit as one scalar residual ${\rho }_i=\sqrt{2J_i}$, with gradient $\nabla {\rho }_i=\nabla J_i/{\rho }_i$. Since ${\sum }_i{J}_i=\frac{1}{2}{\sum }_i{\rho }_i^2$, Gauss–Newton applies with the $N\times n$ Jacobian whose rows are $\nabla {\rho }_i^{\top }$. It combines the $N$ pattern gradients into one search direction, at the price of a coarser curvature model than the full measurement Jacobian.

Properties. Locally quadratic convergence for zero-residual problems and linear convergence for small residuals. The matrix $J^{\top }J$ is extremely ill-conditioned in EIT, so the step must be damped, which gives the Levenberg-Marquardt Method. With a regulariser $\beta \mathcal{R}$ one solves $(J^{\top }J+\beta {\nabla }^2\mathcal{R})\delta =-(J^{\top }r+\beta \nabla \mathcal{R})$.

In ModularEIT.jl: GaussNewton, residual_and_jacobian!, jacobian_operator, jacobian_column_norms.

References

1. J. Nocedal, S. J. Wright (2006). Numerical Optimization, 2nd ed., Ch. 10. Springer. doi:10.1007/978-0-387-40065-5
2. W. R. B. Lionheart (2004). EIT reconstruction algorithms: pitfalls, challenges and recent developments. Physiol. Meas. 25(1), 125–142. doi:10.1088/0967-3334/25/1/021