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?

  • $\begingroup$ I asked a similar question and didn't have a complete answer: scicomp.stackexchange.com/questions/30307/… $\endgroup$ – R zu Oct 25 '18 at 15:15
  • 1
    $\begingroup$ Should the Jacobian be a square matrix? Do you mean something like $\begin{bmatrix} J_{a}(\vec{a})& 0\\ 0& J_{b}(\vec{b})\\ \end{bmatrix}$, or $\begin{bmatrix} J_{aa}(\vec{a}, \vec{b})& J_{ab}(\vec{a}, \vec{b})\\ J_{ba}(\vec{a}, \vec{b}) & J_{bb}(\vec{a}, \vec{b})\\ \end{bmatrix}$? If it is the former, then the linear sub-system is separable from the other parts of the system. That means we can solve the two parts separately. If it is the later, then I guess it is situation in my question. $\endgroup$ – R zu Oct 25 '18 at 15:25

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.

| cite | improve this answer | |
  • $\begingroup$ Thanks, this is an interesting approach. In my application, $A$ and $B(y)$ are indeed long and skinny. However, $A$ is dense and $B(y)$ is sparse...I wonder how this would affect the solution for $\Delta y$, since I need to form $A v$ twice and $A^T v$ once. $A$ is so big in fact that I cannot store it in memory, so when forming $A v$ I need to build $A$ each time from scratch. But its probably worth trying it to see. $\endgroup$ – vibe Oct 27 '18 at 7:54
  • $\begingroup$ How big is $A^T A$ compared to $B^T(y) B(y)$? A fast matrix-vector product with $A$ and $B(y)$ can be (approximately) implemented using finite differences if you have a fast way to evaluate $f(x,y)$. See some of the literature on "matrix-free" Krylov-Newton methods $\endgroup$ – Richard Zhang Oct 29 '18 at 16:19

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.