I want to implement robust line fitting over a set of $n$ points $(x_i,y_i)$ by means of the Least Absolute Deviation method, which minimizes the sum

$$\sum_{i=1}^n |y_i-a-bx_i|.$$

As described for instance in Numerical Recipes, you can find the optimal slope, by canceling the function

$$s(b):=\sum_{i=1}^n x_i\text{ sign}(y_i-bx_i-\text{med}_k(y_k-bx_k))$$

where $b$ is the unknown. To find the change of sign, a dichotomic search is - appropriately - recommended. Anyway, the starting search interval as well as the termination accuracy are based on purely statistical arguments (using a $\chi^2$ test).

I was wondering if there is an analytical way to derive safe bounds for the slopes (i.e. values such that $s(b)$ is guaranteed positive or negative), as well as the termination criterion (minimum of $|s(b)|$). (In fact, I don't even know if the function is monotonic.)

Note that my question is not related to the simplex approach for this problem.

  • $\begingroup$ Why would you want to guarantee that $s(b)$ is positive or negative? And what do you mean by "canceling the function" just before the equation for $s(b)$? $\endgroup$ – jbowman Jan 11 '18 at 21:44
  • $\begingroup$ @jbowman: I mean finding the zero of the function, hence the change of sign. $\endgroup$ – Yves Daoust Jan 12 '18 at 8:19
  • $\begingroup$ @Yves Daoust, are you looking for a good initial value for your algorithm, or something else? $\endgroup$ – Lucas Roberts Jan 14 '18 at 3:21
  • $\begingroup$ @LucasRoberts: I need an initial bracketing. $\endgroup$ – Yves Daoust Jan 14 '18 at 10:01
  • $\begingroup$ @Yves Daoust, right if the initial bracketing is to choose an initial value then many times you can use the least squares estimate(s) as initial values. $\endgroup$ – Lucas Roberts Jan 15 '18 at 15:17

It sounds like you are assuming that the intercept term is known. You can search over all pairs of observations (all $(i,j)$ pairs) where each observation is a pair $(y_i, x_i)$ and take the smallest and largest slope term. Denote each slope term $\beta_{i,j} = (y_i-y_j)/(x_i-x_j)$, then denote the minimum of these as $m$ and the maximum as $M$. The initial bracket would be: $[m,M]$. The initial scan has complexity $\mathcal{O}(n^2)$ which isn't too great. There might be a more clever way to do the scan but this one will work provided $n$ is not too large.

The intuition for this bracketing is that the slope term is smoothing the observed data and therefore must lie between the two observed extremes.

| cite | improve this answer | |
  • $\begingroup$ I don't assume the intercept term known. On the opposite, it must be computed for any value of $b$, as a median. Since I asked the question, I have understood that it suffices to sort the values and compute the slopes between successive points, which takes only $O(n \log n)$. $\endgroup$ – Yves Daoust May 10 '18 at 21:48

Not the answer you're looking for? Browse other questions tagged or ask your own question.