Calibration of a computer code: how to deal with parameter vectors such that for some data points the code doesn't converge

Question

I have an experimental data set $D=\{\mathbf{x}_i,y_i\}_{i=1}^N$, and a computer code with inputs $\mathbf{x}$ and calibration parameters $\boldsymbol{\theta}$, returning a value $s=f(\mathbf{x},\boldsymbol{\theta})$. Assuming Gaussian iid errors, I want to calibrate my code using the data set $D$ using nls or similar, more advanced functions in R. If you need more details, you can find them in this question. Note that since my model is implicitly defined by my computer code and not by an explicit formula, any nonlinear least squares function must use numerical, not analytical, Jacobians.

Now, the problem is that for some vector $\boldsymbol{\theta}^{(j)}$, the computer code doesn't converge at all points of $D$. It converges for some data points, but it doesn't converge for others. From a theoretical point of view, I believe the likelihood $\mathcal{L}(D;\boldsymbol{\theta}^{(j)})$ is not even defined.

Because of various issues, I finally switched to using nonlinear least squares functions such asnlfb from package nlmrt and nls.lm from package minpack.lm, which do not require a formula argument, but take a function resfn which takes in input the parameter vector $\boldsymbol{\theta}$ and return the vector of residuals $\mathbf{r}=(y_1-f(\mathbf{x}_1,\boldsymbol{\theta}),\dots,y_N-f(\mathbf{x}_N,\boldsymbol{\theta}))$. Suppose that for a certain value $\boldsymbol{\theta}^{(j)}$ proposed by the NLS algorithm, my computer code is not converging at points $\mathbf{x}_{i1},\dots,\mathbf{x}_{im}$. Now, what should I return as the corresponding components $r_{i1},\dots,r_{im}$ of the residual vector $\mathbf{r}$?

• Should I return NA? I don't think that would be handled by nlfb or by nls.lm.
• Should I return an arbitrary, very large number $z$? I think this would "encourage" the NLS algorithm to move away from similar parameter values, and move towards other regions of the parameter space, where hopefully my code would converge on all of $D$. Sounds very ad-hoc.

I tried dealing with this issue empirically: I removed those few points in $D$, which correspond to conditions which are not very interesting for calibration and which are also known to pose convergence issues to my computer code. Since I have $N$=319 experimental points, and only 3 parameters to calibrate (for now!), deleting a few points wasn't an issue. However, the NLS regression is still not converging, and I need to solve this issue. Can you give me some tips?

Answer 1

Convergence should always occur for problems that meet certain (rather minimal) requirements (e.g., the existence of a solution). Although it is always risky to generalize, it is common for some data sets to converge when others do not when (in no particular order)

1) the model is not especially well suited to the data 
2) the numerical precision supplied to the regression algorithm is inadequate 
3) the regression method type is not global enough 
4) not enough trial initial parameter values are used or enough combinations of different starting values are attempted 
5) the starting values of parameters are too far from the solution values to converge to the global solution 
6) the regression norm minimized is not matched to the data type 
7) the regression method is not optimized for the data/model combination.

Ad 1) the model is not especially well suited to the data If, for example, we have a biexponential mixture model and are trying to fit blood plasma drug concentration data, the correct convergent answer may be complex number valued. This is because in the extreme case when we have 4 data entries and 4 parameters to fit, the biexponential model cannot always solve for those points. However, if we then allow for complex field numerical processing a convergent answer will always be found. Moreover, the are $n$-tuple such solutions, and they (obviously) lack uniqueness. Thus, if we have more than 4 data entries and the same problem persists, a convergent answer (the likes of which we are duty bound to find by applying enough elbow grease) is the least of our problems. In that case, the major problem we have is a mismatch between data and fit model.

Ad 2) the numerical precision supplied to the regression algorithm is inadequate It is quite common to need many more places of machine precision (generally twice as many) than the precision of our converged answer. Suppose our data entries are barely accurate to 3 places. We still need to zero pad them, sometimes to 65 places, just to get our regression to converge. Often, if an algorithm does not converge with 20 place machine precision, it will at 30.

Ad 3) the regression method type is not global enough For well behaved problems Nelder-Mead is generally fast and convergent. For not so well behaved problems, we may need to cover the entire parameter space with starting values in a random search routine, which will then always converge globally (but take a long time to do so.)

Ad 4) not enough trial initial parameter values are used or enough combinations of different starting values are attempted Not much to say here beyond that if one did not try enough, try again.

Ad 5) the starting values of parameters are too far from the solution values to converge to the global solution Getting stuck in a non-global solution region is common and the complexity of searching for $n$ parameters requires searching in an $n$-space topology, messy that is, at the best of times. Thus, being close enough to a solution to begin with, avoids getting lost in that topology.

Ad 6) the regression norm minimized is not matched to the data type Alternatively, one can sometimes transform the data to match the regression type. But more generally, changing the norm is more flexible. This is a more advanced topic than many users would consider. However, for example, if the data is proportional type one should be minimizing $||1-\frac{y}{f(x)}||$ not $||f(x)-y||$.

Ad 7) the regression method is not optimized for the data/model combination There is a lot more to regression than minimizing goodness-of-fit. Sometimes, to obtain convergent solutions, it is necessary to abandon goodness-of-fit entirely and use a regression target that is specific to the question being asked. For example, if we want to know what kidney function a patient has, goodness-of-fit to plasma concentrations of an exogenous kidney function marker drug will not help us, but a minimum error of relative kidney function regression target could.