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

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?

• nlfb is Nash variant of Marquardt nonlinear least squares solution via qr linear solver. The nls.lm function provides an R interface to lmder and lmdif from the MINPACK library, for solving nonlinear least-squares problems by a modification of the Levenberg-Marquardt algorithm, with support for lower and upper parameter bounds. The implementation can be used via nls-like calls using the nlsLM function. Well, yuk. Maybe think better methods for a start. – Carl Mar 30 '17 at 9:13
• @Carl they seem quite good methods to me. It's not the NLS algorithm which fails to converge, it's the code which I'm trying to calibrate. Anyway, which better methods would you suggest for a start? – DeltaIV Mar 30 '17 at 12:47
• Depends on the data, no one method suits all. Nelder-Mead is my default for well behaved problems (fast). Some people prefer conjugate gradient. Random search for poorly chosen models (slow). Also, I do not use least squares without testing if appropriate, and often is not. Last, explain please what is not converging. – Carl Mar 30 '17 at 16:54
• @Carl ok, tomorrow I will add some convergence histories to my post. Did you read my other question I link to in this post, and also its comments? I think it helps clarifying the problem a bit. Anyway, the concept is this: I'm trying to calibrate not a "usual" statistical model, such as a linear regression model, or a spline, etc. but a computer code. I don't know which field you are from, but if you're familiar with Computational Fluid Dynamics codes, or Finite Element Analysis codes, etc., the idea is that I have a computer code (an executable), actually, which I can run for [1/2] – DeltaIV Mar 30 '17 at 17:43
• [2/2] arbitrary values of its inputs, by changing an input file. The code writes an output file, which I read. Based on years of domain knowledge at my company, with the support of the design engineers I divided the code inputs in inputs, which I take from the tests, and calibration parameters, which I want to calibrate with NLS, as an introduction to a future calibration using MCMC. So the code is the "function" $s=f(\mathbf{x},\boldsymbol{\theta})$ I want to calibrate. I don't have an explicit, analytical expression for the function, but I can evaluate it for arbitrary values ctd... – DeltaIV Mar 30 '17 at 17:47

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.

• Thanks for the answer, but what's not converging here is not the NLS algorithm, it's the computer code (executable) I'm trying to calibrate with NLS. I still would like to know what to do with those points for which I cannot compute the residuals, because my executable doesn't run: set the residual to a very high value? to Not a Value? Do something completely different? – DeltaIV Mar 30 '17 at 7:49
• @DeltaIV nls is nonlinear least squares for those who do not code R-language. The answer provided is not R-language specific. A lack of regression convergence is a failure on the part of the user that means that a regression problem was either ineptly, or more often-inappropriately, posited. One should select a regression algorithm by first answering "What is the exact purpose of applying regression, i.e., the goal or information to be extracted?" – Carl Jan 3 '19 at 19:30