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?

  • 3
    $\begingroup$ With 10 million observations, your problem is amazingly over constrained. You would probably get perfectly good results with a random sampling from those 10 million observations- this could speed up your computations significantly. $\endgroup$ – Brian Borchers Dec 15 '19 at 3:34
  • 1
    $\begingroup$ Automatic Differentiation (AD) tools could prove helpful in this situation. However, these tend to be very language specific. What language are you implementing this in? $\endgroup$ – Brian Borchers Dec 15 '19 at 3:34
  • $\begingroup$ I'm using C and I've never used AD tools $\endgroup$ – vibe Dec 15 '19 at 3:54

Your Answer

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

Browse other questions tagged or ask your own question.