I have a nonlinear least squares problem, in which I am trying to minimize residuals which can be divided into four classes: $$ \min_x ||\epsilon(x)||^2 + ||\xi(x)||^2 + ||\delta(x)||^2 + ||s(x)||^2 $$ The parameter vector $x$ can be naturally divided into four classes of parameters, $x = \left( x_1, x_2, x_3, x_4 \right)$. Some of the residuals are linear in the parameters, and others are not. The full Jacobian matrix is: $$ J = \pmatrix{ \mathcal{E}_1 & \mathcal{E}_2 & 0 & \mathcal{E}_4(x) \\ \Xi_1 & \Xi_2 & \Xi_3 & 0 \\ \Delta_1 & \Delta_2 & 0 & 0 \\ S_1(x) & S_2(x) & S_3(x) & 0 } $$ So you can see some blocks in the Jacobian don't depend on the parameters $x$ (these correspond to the residuals which are linear in $x$). So for example, $\mathcal{E}_i = \partial \epsilon / \partial x_i$.
Now since its a nonlinear least squares problem, I need to iterate, and solve the following linear least squares problem at each iteration: $$ \min_p \left| \left| \pmatrix{ \epsilon(x) \\ \xi(x) \\ \delta(x) \\ s(x) } - J p \right| \right|^2 $$ and then set $x := x + p$.
I want to take advantage of the sparse structure of $J$, and also the parts of $J$ which don't depend on $x$ to make the calculation as fast as possible. Another issue, is that $J$ is too large to fit in memory at once, so I am using the normal equations approach. $J$ has roughly 10 million rows and 20,000 columns.
With the normal equations approach, I can separate $J^T J$ into a part which can be precomputed prior to the iteration, and a part which needs to be computed on each iteration: $$ J^T J = \pmatrix{ \mathcal{E}_1^T \mathcal{E}_1 + \Xi_1^T \Xi_1 + \Delta_1^T \Delta_1 &&& \\ \mathcal{E}_2^T \mathcal{E}_1 + \Xi_2^T \Xi_1 + \Delta_2^T \Delta_1 & \mathcal{E}_2^T \mathcal{E}_2 + \Xi_2^T \Xi_2 + \Delta_2^T \Delta_2 && \\ \Xi_3^T \Xi_1 & \Xi_3^T \Xi_2 & \Xi_3^T \Xi_3 & \\ 0 & 0 & 0 & 0 } + \pmatrix{ S_1(x)^T S_1(x) &&& \\ S_2(x)^T S_1(x) & S_2(x)^T S_2(x) && \\ S_3(x)^T S_1(x) & S_3(x)^T S_2(x) & S_3(x)^T S_3(x) & \\ \mathcal{E}_4(x)^T \mathcal{E}_1 & \mathcal{E}_4(x)^T \mathcal{E}_2 & 0 & \mathcal{E}_4(x)^T \mathcal{E}_4(x) } $$ So the first term doesn't depend on $x$, is small enough to fit in memory, and can be precomputed. The second term does depend on $x$ and must be computed on each iteration. I wrote only the lower triangles since these matrices are symmetric. Once I have the full $J^T J$, I use a Cholesky decomposition to compute the step $p$.
My question is: is the above approach the best way to implement something like this?
In the future, I would like to add more parameters to my model (i.e. $x_5, x_6, \dots$) which would increase the number of columns in my Jacobian. But then the $J^T J$ matrix will get more and more complicated, and it is already becoming difficult to compute all the different terms in the matrix (without making mistakes!). Does anyone have any suggestions for a strategy to compute all of these terms efficiently and simply? Is there a better method I should be using, such as TSQR?
