"Least square root" fitting? A fitting method with multiple minima

Question

I am looking for the name of a fitting method that will work even if points from multiple dataseries are meshed together.

As far as I understand there are two major methods, least squares and least absolute error; the method I am looking for would be the next step in the direction from LS to LAE:

• Least Squares fitting: Minimize the square of the distances of all data points to a function. The further a point is, the stronger it's deviation is weighted. May be strongly influenced by outliers. Will always converge to one globally unique solution (at least for simple curves), easy to find with gradient descent methods.

---

• Least absolute deviations: Minimize the sum of the distances of all datapoints. All deviations are weighted equal. More robust to outliers, may converge to potentially infinite solutions that all have the same total deviation

(generated by hand; images may not be accurate)

---

• "Least Square Root" ??? (Google does not really have many hits for this name). I am looking for a function that may converge towards multiple local minima, and weights deviations stronger the closer a point is. Which solution it would converge would depend on the starting parameters. I think a least-square-root of all deviations would do this, but I cannot find any reference to this sort of fit.

Is there an official name for a method like this, and what would it be called? Or is this sort of fit fundamentally not possible with curve fitting methods?

I am not looking for a method that would find all local minima, one local minimum depending on starting parameters would suffice.

---

(note that in the example pictures the functions were all lines, but ideally I am looking for a method that is able to iteratively fit points towards arbitrary n-dimensional functions)

Answer 1

It seems to me you are actually looking for a regression models with re-descending ("far away points are weighted less than close ones") loss function. Such loss functions --for example the Tukey biweight-- lead to highly non-convex optimization problems, meaning that there are, potentially, a finite but factorial-order increasing number of potential solution. Clearly, which ever one you will pick up will depend on your "starting weights" (you initial guess about which observations belong to the model you are looking for).

R has some facilities to fit such models. An example is the lmrob function in the robustbase package. For illustrative reasons I use a univariate example below but in principle the logic extends to higher dimensions.

library(robustbase)
#constructing the data --one group at a time
x1<-cbind(1,rnorm(100,0,5))
x2<-cbind(1,rnorm(100,0,5))
b1<-c(-3,2)
b2<-c(4,-2)
y1<-x1%*%b1+rnorm(100)
y2<-x2%*%b2+rnorm(100)
#mixing the data
x<-rbind(x1,x2)
y<-c(y1,y2)
sam1<-sample(1:100,3)  
bet1<-lm(y[sam1]~x[sam1,]-1)
ini1<-list(coefficients=as.numeric(coef(bet1)),scale=sd(bet1$resid))
ctrl<-lmrob.control()
mod1<-lmrob(y~x-1,init=ini1)
    
plot(x[,2],y,pch=16)
points(x[sam1,2],y[sam1],col="red",pch=16)
abline(mod1$coef,col="red")

$d+1$ points from the same relationship. The red line is the final fitted regression function" />

Now, as I mentioned, fitting these types of model is actually a highly non-convex problem. The final solutions will depend on the quality of your starting points. In the example above, my "guess" that the original points (marked in red in the plot above) belonged to the same relationship was correct. In general, specially in higher dimensions, this is a highly questionable assumption. None-withstanding this, if your initial guess is correct, the biweight loss function enabled an appreciable gain over the naive solution (the OLS line passing by the 3 initial points).

EDIT:

I forgot to mention that you have a parameter that determines the width of the area for which the bi-weight gives non-zero weights. This is the tuning.psi parameter of the lmrob.control object.

There is a tradeoff between efficiency and the width of the biweight: the wider the biweight, the more points you effectively use, the more efficiency you get. A limiting case is a biweight that never re-descends (maximal width) which would give you the OLS solution. Adjusting the width is done by:

ctrl$tuning.psi<-1.35
mod1<-lmrob(y~x-1,init=ini1,control=ctrl)

If you set the ctrl$tuning.psi too small, you will get convergence problems. This can be solved by increasing the max.it value of the control object:

ctrl$max.it<-500

There is a whole theory on optimal values of the tunning constant, but it only applies if the sub-population you are targeting has more than half of the observations. I gather this is not the case you are concerned with. If this is the case, I think it is best is to play with to get a handle on it.

Answer 2

I have worked on this problem over the last two years. My co-authors and I called it a "data association problem", as part of the problem consists in finding the correct association between data points that belong to the same curve. The other part would be performing the regression itself on both data sets.

As user user603 correctly points out, in general this problem corresponds to a difficult non-convex optimization problem. In your case it may be quite simple, as your data correspond to straight lines. If you simply want to do this - detect lines - you could try the Hough transform.

I worked on the scenario where the data are nonlinear curves, without restrictions, as in this example figure:

Check out my recent publication "Overlapping Mixtures of Gaussian Processes for the data association problem" if you'd like, it contains much background information on this kind of problem. I have some other articles on this, but they have not been accepted for publication yet.

By the way I'd love to hear where you got this data.