# Non-linear least absolute deviation regression with multiple global minima

I am fitting a single exponential decay formula with three parameters (a,b,c): y ~ $a \exp(-xb) + c$ using the LAD cost function: $\min \sum |(y - f(x))|$. $x$ is in units of time (as is $b$), and $a$, $c$ are unitless. Measurements are taken far enough out in time that the last $y$ is effectively only noise from measurement process.

Sometimes we have data points that are corrupted by sporadic errors that are most definitely not normally distributed when they occur at all, so this method seems a very reasonable approach to use. I include an example using simulated data below: Ignoring parameter $c$, I plot the two-parameter cost function graphically, and it is easy to see that there is a plane of points in the parameter space that minimizes the cost function (see figure). I have confirmed this for the three-parameter model by brute force searching the parameter space. My questions is: Is there any way, given the data, the model and assuming I can find all of the reasonable minima, to choose one of those optimal parameter sets as the representative parameter set? Right now I am choosing the centroid of the hyperplane (if that's an appropriate term) of the minimizing parameters, and that's getting me extremely close to R -> quantreg -> nlrq's solutions with tau=0.5 in about 90% of the cases.

• You could just use a gradient-descent algorithm, starting from a random guess, which should often converge to points on your optimal "plane", because there are only three parameters to optimize. Alternatively you could use a non-gradient optimizer (like in package pso). – jubo Aug 14 '14 at 19:01
• jubo -- thanks for the comment. I can find the minima no problem, just that there's a bunch of them I am wondering if there is a way to choose the best, since they all equally minimize the cost fcn. – wtt Aug 14 '14 at 19:23
• wtt if they all have equal values of the cost function, in what other sense do you mean 'best'? What's 'best' for your purpose? – Glen_b -Reinstate Monica Aug 14 '14 at 21:02
• wtt, if you want to find all possible local minima, you could try a hybrid approach of using both genetic algorithm and a gradient descent algorithm. The way this hybrid approach would work is Genetic algortithm finds near global minima, use this as a starting value for a local optimizer for such as gradient descent. Run this multiple times, you might be able to find most of the minima – forecaster Aug 14 '14 at 21:43
• Glen_b, that's a good question. The only other metric I can think of is that the parameter 'b' only makes sense in a certain range, but in this case the range ~ [20,40] is quite reasonable – wtt Aug 15 '14 at 11:35

You get the same issue as you get even in 1-D. Consider trying to find the value of $\stackrel{\sim}{\mu}$ minimizing $\sum_{i=1}^{n}|y_i-\stackrel{\sim}{\mu}|$ where $n$ is even - it's minimized by any value between the two central order statistics of $y$.

The centroid is a perfectly reasonable choice in your case. It may well be what I'd have used. [However with nonlinear regression you might have to consider whether it's possible to get a region that's not convex (and so, in turn the possible, if unlikley, consequence that the centroid isn't actually in the region). It may be that you can't get a non-convex region, though.]

To make some optimal ('best') choice, you have to say what you would like to optimize. You've already optimized your cost function, and have a region of points equal on that criterion. What other criterion would you use to optimize among them that would be 'best' for your purposes?

(There's a quantity that the centroid optimizes, for example; it may well be 'best' in a sense that you might care about.)

A simple solution would be to transform your cost function into:

$cost(\beta, x, y)=\sum_{i=1}^{N} \left|y - \hat y(\beta, x)\right| + \lambda \sum_{k=1}^{P} \left|\beta_k\right|^\alpha$

The second term is a so-called regularization term, the goal of which is to get small values of estimated parameters $\hat\beta$. Common choices of $\alpha$ are 1 and 2 (called the L1 and L2 regularizations, used respectively in Lasso and Ridge regression). A value of 1 will tend to set parameters to 0, a value of 2 tends to just shrink the larger parameters.

The philosophy behind this can be viewed as Occam's razor: small values of $\beta$ yield simpler models, which is desirable (especially in your case). It also has Bayesian justifications (eg. the cost function of ridge regression arises from putting a Gaussian prior on $\beta$ in usual least-squares regression). You then have to find a "good" value of $\lambda$ (the parameter that trades off model complexity and goodness-of-fit), eg. by cross-validation.