Inspired by A Recipe for Training Neural Networks, here are some tips and tricks for curve fitting. Assume we are fitting a non-linear, parametric function \(f(x; \boldsymbol{\theta})\) to data \(y\). The problem could be anything. Many highly technical problems are illuminated by curve fitting, but it helps to think of everyday examples such as: personal savings vs. time, apartment price vs. square meters, salary vs. experience, strength vs. bodyweight, world population vs. year, etc.

Modeling

• Determine the purpose. Do you want to extrapolate, or do you want to interpolate? Is prediction the overall goal, or is it inference of latent variables? In other words: raw predictive power or understanding underlying phenomena? Is uncertainty needed? If yes, then what kind? A posterior conditional distribution \(p(x | \hat{\boldsymbol{\theta}})\) or a posterior predictive distribution \(p(x | \boldsymbol{\theta}) p(\boldsymbol{\theta} | X)\) that accounts for parameter uncertainty? Is curve fitting a one-time job, or will it every night or on demand? Will it run on a single data set, or perhaps on thousands of datasets?
• Read the literature. Before deciding on the function \(f\), search the web and review papers. What functions are used on similar problems? Avoid ad-hoc models like unregularized high-degree polynomials or very intricate functions. Apply Occam’s razor.
• Know basic functions and how to combine them. Know the behavior of basic building blocks like linear functions, powers, exponentials, logarithms, sigmoids, hyperbolic functions, etc. Sometimes you can transition between them in more natural ways than the convex combination \(\alpha f(x) + (1- \alpha) g(x)\). For instance, the power transform \((x^\alpha - 1) / \alpha\) lies between the affine function \(x-1\) when \(\alpha=1\) and the logarithmic function \(\ln x\) when \(\alpha \to 0\). Similarly, the generalized pareto distribution lies “between” the exponential \(e^{-x}\) and the hyperbola \((1 + x)^{-1}\). Probability density functions and their antiderivatives are a rich source of modeling functions. Know how to combine functions to create new models.
• Use an identifiable model. Avoid using non-identifiable models like \(f(x; \boldsymbol{\theta}) = \theta_1 \exp(x - \theta_2)\), since they do not have a unique solution. This seems obvious, but I’ve seen papers where non-identifiable models are used without priors.
• Choose an appropriate loss. The least squares loss \(\mathcal{L}(\boldsymbol{\theta}) = \sum_i (y_i - f(x_i; \boldsymbol{\theta}))^2\) has nice properties, but might not be what you want. Know when to use \(p\)-norms, huber loss, pseudohuber, geometric loss, harmonic loss and statistical losses. Remember that fitting log-data to log-model with a loss like \((\log y_i - \log f(x_i; \boldsymbol{\theta}))^2\) is not equivalent to fitting data to a model without logs. Determining the appropriate loss function is a big part of the job. Always align the primary metric with the loss. If you care about MAE as the metric, it makes little sense to use RMSE as the loss.
• Use a prior, or if the model is not statistical, a regularization term. A prior encodes prior knowledge and helps the solver numerically. It also lets you get a solution when you have little data or no data at all. The 2-norm \(\lVert \boldsymbol{\theta} - \boldsymbol{\mu}_{\boldsymbol{\theta}} \rVert_2\) is a good start and \((\boldsymbol{\theta} - \boldsymbol{\mu}_{\boldsymbol{\theta}})^T A (\boldsymbol{\theta} - \boldsymbol{\mu}_{\boldsymbol{\theta}})\) can be even better. Often \(A\) is diagonal.
• Deal with parameter constraints. There are several options: (1) use an optimizer that handles constraints natively, (2) re-parametrize the function as \(f(x; t(\boldsymbol{\theta}))\) with \(t\) mapping from unbounded real numbers to the desired bounds or (3) add a barrier function to the objective. For instance, if \(\theta_1 \in [0, 1]\), then you can (1) use an algorithm like BFGS-B (the B stands for “box constraints”), (2) transform with e.g. a sigmoid like \(\theta_1 = t(z) = (1 + e^{-z})^{-1}\) or (3) add barriers such as \(-\log(\theta_1) - \log (1 - \theta_1 )\) to the objective. The log-barrier is equivalent to a beta-prior.
• Weight the data. Consider a weighted loss like \(\sum_i w_i (y_i - f(x_i; \boldsymbol{\theta}))^2\) if (1) some data points are more important than others or (2) you are working on a temporal forecasting problem. Beware of the difference between (a) inverse variance weights and (b) weighting the log-likelihood. They are the same when estimating the mean of a normal, but not in general (see also this post).
• Consider partial pooling. Use hierarchical models to automatically regularize with hyperpriors. This is relevant if you have many curve fitting instances that you can conceptualize as being drawn from some population, and you want to learn the distribution of that population. It’s also extremely useful if some instances have little data.

Validation

• Cross validate. Split randomly if you’re interpolating. Use a time based split if you’re forecasting. Always match cross validation to how the model will be used. Compute metrics and study them.
• Use your eyes. Curve fitting can be visualized! Sample prior curves and plot them, plot the initial guess, curve fit “by hand” a few times, fit on every instance you have in your database. Visually look for errors, edge cases, etc.
• Implement the forward model. Implement data generation based on your model if it’s statistical. This is called a prior predictive check. Can you tell a simulation apart from real-world observed data? How? If you can, does it give you ideas for how to improve the model?
• Test parameter estimation. Some property-based tests: (1) generate a random \(\boldsymbol{\theta}\), use it to generate a dataset, fit on the dataset. Is \(\lVert \boldsymbol{\theta} - \hat{\boldsymbol{\theta}} \rVert\) small? It should be, especially in the limit of much data. (2) Test that \(\mathcal{L}(\hat{\boldsymbol{\theta}}) \leq \mathcal{L}(\boldsymbol{\theta})\).
• Plot parameter-space. You can plot the contours of the loss as a function of two parameters in \(\boldsymbol{\theta}\) at a time, while fixing the remaining parameters at \(\hat{\boldsymbol{\theta}}\). Sometimes studying the Hessian is worthwhile.
• Plot residuals. When are residuals large? Are residuals homoscedastic? Should you model heteroscedasticity?

Optimization

• Understand the optimization algorithm. Have at least a rudimentary understanding of the algorithms you use. Read the documentation and perhaps a summary paper. Play with the optimizer on toy problems. Get a feeling for how it works, and more importantly how and when it fails. Study the full verbose output of the algorithm.
• Put variables on the same scale. First order methods (gradient methods) depend critically on the condition number. Pure second order methods (Newton’s method) on the other hand are affine invariant and in principle unaffected by the condition number (numerics and quasi-Newton methods make this untrue in reality). Standardize or take logs/exp/sigmoid/etc as needed.
• De-correlate variables. Sometimes you can de-correlate variables through re-parametrization of the curve function.
• Use a prior. Even if you lack prior knowledge in the modeling sense, weak regularization tends to help the optimizer.
• Pick a good starting point. Graphical trial and error often works. You can also solve a simpler problem, e.g. solve a least squares problem first and use that solution as the starting point for your \(p=1\) norm problem. The curve function \(f\) can sometimes be relaxed by removing non-linearities. You could use (y[-1] - y[0])/(x[-1] - x[0]) as initial guess for slope, mean(y) as initial guess for bias, etc. Stress test the solver on bad initial solutions too.
• Implement gradients and/or Hessian. Take the time to work out gradients. Use a CAS like SymPy or WolframAlpha along with pen and paper. Verify gradients numerically with scipy.optimize.check_grad or similar functions. If the gradient is too unwieldy, consider autograd with one of the many available packages: jax, autograd, PyTorch, TensorFlow or Stan.
• Use numerically stable implementations. Watch out for expressions like \(\ln(1 + e^{\theta})\) that overflow, catastrophic cancellation and other phenomena. See The Right Way to Calculate Stuff for more examples.

Curve fitting seems easy, but this perception can be deceptive. It’s not uncommon that people fail to get the desired results, often because they call standard routines without knowing how they work. I hope the tips and tricks above can help. Experience helps too: work on diverse practical problems and try to solve them. Make notes on what works in practice.