Linearizing a nonlinear least squares model

[vibe3]

I have a nonlinear least squares problem with a set of parameters $\bf{g}$, where I need to minimize the function:
$$\chi^2 = \sum_i \left( y_i - M(t_i ; {\bf g}) \right)^2$$
The $t_i$ are some independent parameters associated with the observations $y_i$ and the model function has the form
$$M(t_i ; {\bf g}) = \sqrt{X(t_i; {\bf g})^2 + Y(t_i;{\bf g})^2}$$
The functions $X(t_i;\bf{g})$ and $Y(t_i;\bf{g})$ are linear in the model parameters $\bf{g}$, ie:
$$X(t_i;{\bf g}) = \sum_k g_k X_k(t_i)$$
and
$$Y(t_i;{\bf g}) = \sum_k g_k Y_k(t_i)$$
In order to solve the nonlinear least squares problem, I need to construct the matrix $J^T J$ at each iteration, where $J$ is the Jacobian:
$$J_{ij} = - \frac{\partial}{\partial g_j} M(t_i;{\bf g})$$
My question is can anyone see a clever way to optimize the computation of the matrix $J^T J$, since all of the $X_k(t_i)$ and $Y_k(t_i)$ can be precomputed? I have millions of observations (rows of the Jacobian) and so its extremely slow to compute each row of the Jacobian and add it into the $J^T J$ matrix. I'm hoping that the linearity of the functions $X$ and $Y$ might allow some way to linearize or precompute large portions of the Jacobian matrix, but I don't see an easy way to do this.

 

[jvicens]

I understand your frustration with the slow computation of the Jacobian matrix. However, there are a few ways that you can optimize the computation process.

One approach is to use a sparse matrix representation for J^T J. This means that instead of storing the entire matrix, you only store the non-zero elements. This can greatly reduce the memory and computational requirements, especially if the majority of the elements in the matrix are zero.

Another approach is to use a forward-mode automatic differentiation algorithm, which can efficiently compute the Jacobian matrix in a single pass. This method is particularly useful if your problem has a large number of parameters and a relatively small number of observations. However, it may not be the best option for your specific problem since you have millions of observations.

Additionally, you can try to parallelize the computation of the Jacobian matrix. This means using multiple processors or threads to compute different rows of the matrix simultaneously, which can significantly speed up the process.

Finally, you can also consider using a different optimization algorithm that does not require the computation of the Jacobian matrix, such as a genetic algorithm or a particle swarm optimization algorithm. These methods may be more computationally intensive, but they can often handle large amounts of data more efficiently.

In conclusion, there are several ways to optimize the computation of the Jacobian matrix for your nonlinear least squares problem. Depending on the specific characteristics of your problem, some of these approaches may be more effective than others. I recommend experimenting with different methods to find the most efficient solution for your particular problem.