Nonlinear least squares when some parameters are linear

Question

Consider the least squares problem, $$ \min_{\mathbf{a},\mathbf{b}} || \mathbf{f}(\mathbf{a},\mathbf{b})||^2 $$ where $\mathbf{a},\mathbf{b}$ represent the unknown parameters to be found. In my problem, the parameters $\mathbf{a}$ appear linearly in the residual vector $\mathbf{f}$ while the parameters $\mathbf{b}$ are nonlinear. So my Jacobian matrix will look like: $$ J = \pmatrix{J_a & J_b(\mathbf{b})} $$ In other words, part of the Jacobian matrix ($J_a$) is fixed and will not change from one least squares iteration to the next. Can this somehow be exploited to speed up the iterative algorithm? On each iteration the following linear least squares problem needs to be solved for the step $(\delta \mathbf{a}, \delta \mathbf{b})$: $$ \pmatrix{J_a & J_b(\mathbf{b})} \pmatrix{ \delta \mathbf{a} \\ \delta \mathbf{b}} = -\mathbf{f}(\mathbf{a},\mathbf{b}) $$ The only solution I've come up with so far is to use the normal equations approach, and the normal equations matrix is, $$ J^T J = \pmatrix{ J_a^T J_a & J_a^T J_b(\mathbf{b}) \\ J_b^T(\mathbf{b}) J_a & J_b^T(\mathbf{b}) J_b(\mathbf{b})} $$ The upper left term $J_a^T J_a$ can be precomputed once while the other terms need to be computed for each iteration. Then I need to solve $$ (J^T J) \pmatrix{ \delta \mathbf{a} \\ \delta \mathbf{b}} = -J^T \mathbf{f} $$ and I can save a little bit of computation when constructing $J^T J$ using the precomputed $J_a^T J_a$.

Just wondering if anyone knows of any other clever methods to solve such a problem. Ideally I'm wondering if its possible to somehow remove the $\mathbf{a}$ parameters and solve a smaller sub-problem involving only $\mathbf{b}$ at each iteration?

Answer 1

One simple way is to pre-factorize the block associated with $J_{a}$. For readability, I will slightly change your notations. Consider solving $$ \min_{x,y}\|f(x,y)\|^{2}\text{ with Jacobian }\nabla f(x,y)=\begin{bmatrix}A & B(y)\end{bmatrix}. $$ We can solve each Gauss-Newton step of the form $$ \begin{bmatrix}A^{T}A & A^{T}B(y)\\ B^{T}(y)A & B^{T}(y)B(y) \end{bmatrix}\begin{bmatrix}\Delta x\\ \Delta y \end{bmatrix}=\begin{bmatrix}c\\ d \end{bmatrix}, $$ by block-LU factorization $$ \begin{bmatrix}A^{T}A & A^{T}B(y)\\ B^{T}(y)A & B^{T}(y)B(y) \end{bmatrix}=\begin{bmatrix}A^{T}A & 0\\ B^{T}(y)A & S(y) \end{bmatrix}\begin{bmatrix}I & (A^{T}A)^{-1}B(y)\\ 0 & I \end{bmatrix},\tag{1} $$ where the Schur complement reads $$ S(y)=B^{T}(y)B(y)-B^{T}(y)A(A^{T}A)^{-1}A^{T}B(y). $$

Specifically, let us precompute the permuted Cholesky of $A^{T}A$, as in $$ A^{T}A=\Pi(LL^{T})\Pi^{T} $$ with permutation matrix $\Pi$ and sparse lower-triangular $L$ with positive diagonals. Then we can invert $A^{T}A$ via two sparse triangular solves: $$ (A^{T}A)\xi=g\quad\iff L\eta=(\Pi^{T}g),\quad L^{T}(\Pi^{T}\xi)=\eta.\tag{2} $$

Combined, we solve $\Delta x$ and $\Delta y$ by inverting the block in (1) analytically, while using (2) to compute matrix-vector products with $(A^{T}A)^{-1}$, as in $$\begin{align*} z & =(A^{T}A)^{-1}c\\ \text{solve }S(y)\Delta y & =d-B^{T}(y)Az\text{ for }\Delta y\\ \Delta x & =z-(A^{T}A)^{-1}B(y)\Delta y \end{align*}$$ Note that this is only efficient if $B(y)$ is long-and-skinny, i.e. has many more rows than columns.

We can also apply the same idea to a QR-based technique. The resulting complexity is the same, though in practice it is slower but much more accurate.