# Total least squares curve fit problem

I am trying to fit a quadratic curve across a scatter plot of two variables. Since both variables are noisy I cannot use an ordinary least square regression (OLS) and I would like to have a symmetrical answer regardless which one I choose to be the independent variable.

Total (orthogonal) least squares regression seem to be the answer.

Unfortunately it seems that the two implementations that I have tried to use so far give me very bad fits.

In the figure below the blue line is the OLS fit, which obviously could be improved. The red is the TLS fit using the matlab/octave code below which seems to be the standard approach using single-value decompositions (SVD). I copied it from the corresponding Wikipedia article here.

I can see how it is trying to fit the upper-right denser cloud of points to the (very large) expense of the long tail. I have tested the code with simulated data and it works so is does not seem to be a bug.

I have also tried this other matlab package that implements a different iterative approach but unless I give it very tight bounds to the model parameters also gives me "crazy" fits.

I want to be able to perform TLS between many variables with many bootstrapping repeats to obtain confidence intervals and I cannot afford to come out with bounds for all those repeats eyeballing each graph; that would also tamper with the confidence interval estimation I guess.

I reckon the problem is that the data is a bit too noisy but I can do little about it. You can find the data here.

Any ideas how to go around this issue? Perhaps some know a more robust TLS method with an easy to use implementation available? There seem to be many academic articles but little runnable code.

function B = tls(xdata,ydata)
% xdata = nx3 matrix [x.^2 x ones(size(x))]
% ydata = nx1 matrix
m       = length(ydata);       %number of x,y data pairs
X       = [xdata];
Y       = ydata;
n       = size(X,2);          % n is the width of X (X is m by n)
Z       = [X Y];              % Z is X augmented with Y.
[U S V] = svd(Z,0);           % find the SVD of Z.
VXY     = V(1:n,1+n:end);     % Take the block of V consisting of the first n rows and the n+1 to last column
VYY     = V(1+n:end,1+n:end); % Take the bottom-right block of V.
B       = -VXY/VYY;
end

• Could you post the data ? Does the curve change if you centre to [0 0] and scale X and Y the same ? Sounds like you really want Scale invariant TLS – denis Aug 6 '14 at 15:21
• @denis I have added a link to the data in my last question edit (second last paragraph). x and y are both in mm but corresponding to objects of a rather different size so the 10-fold difference. Thanks for pointing that out! perhaps I need to scale the data into roughly the same size before fitting? – Valentin Ruano Aug 6 '14 at 15:39
• @denis yeah! it looks much better after scaling!!! for each variable I subtract the mean and divide by its sd and then the code above works, well at least for this case.... I will post and updated picture when I have time. Could you please post an answer so I can give you credit? – Valentin Ruano Aug 6 '14 at 21:33

Curves $y = a + b x + c x^2$ are parabolas that point straight up, so cannot match tilted parabolas (think of a big satellite antenna) no matter how you find $a\ b\ c$.

Edit 9 August: the 45° rotation in the original answer below is wrong, @whuber is right.
Consider noisy parabolas, all with $\bar{X} = \bar{Y} = 0$ and $s_X = s_Y$, tilted at various angles: 0° right, 45°, 90° up ... I see no direct way of finding the tilt. Brute force is klutzy, but may be good enough:

XY -= mean
for angle in e.g. [0 10 20 .. 170]:
rotate XY by angle
scale, XY /= std
TLS: Y = XB + Res
fit = rms residual_i = |data_i - nearest point on the parabola|
take the best fit.


If more accuracy is needed, use a 1d optimizer such as golden search.
A quite different approach would be to minimize a nonlinear function of a b c and tilt.

In short, the question of how to use TLS to directly fit a noisy tilted parabola seems to me open, even in 2d.

A method that might work well enough for shallow parabolas:

1. centre the data at [0 0] -- subtract the means
2. scale X and Y the same -- divide by their sd s
3. now the 45° line $y = x$ is pretty good line fit; see methods-for-fitting-a-simple-measurement-error-model (or it might be $y = - x$). Rotate the data 45° clockwise, so that a parabola now points up
4. least-squares fit a parabola, by TLS (or Ordinary ? try-it-and-see).
5. reverse steps 1 - 3: rotate 45° counterclockwise, unscale, shift back.
• I see why this might approximate a solution for datasets that are approximately linear. However, its application to the general problem that was posed seems limited. More importantly, what one hopes to get out of any statistical analysis includes an assessment of the uncertainty in the results. How would you obtain that using the procedure you have outlined? – whuber Aug 7 '14 at 14:11
• @whuber What I have in mind is to compare a strait line with a curve (in this case a poly-2) and was thinking that creating a confidence intervals on 'c' by bootstrapping the data would be sufficient. But honestly not sure of that, any ideas are welcome. – Valentin Ruano Aug 7 '14 at 14:59
• @whuber, I'm afraid I just plot the residuals $X_i, res_i$ from the TLS split Y = XB + res, and print e.g. Totallsq: residuals / Y 3 % = 0.3 / 10 rms. If that helps to improve the fit or understand the data a bit better, I'm satisfied (satisficed ). If there are better ways of analyzing and visualizing the residuals for TLS, please comment; if not, is TLS unsound ? – denis Aug 8 '14 at 9:33
• Valentin Ruano, could you possibly put up a plot of $X_i, res_i$ (just vlines -- don't know if matlab has a function for the lines to nearest-on-curve) ? Thanks – denis Aug 8 '14 at 9:41
• Denis, I'm not even sure what the question is about. To me (and most mathematicians I believe), a "quadratic curve" is the zero set of a function of the form $z(x,y)=ax^2+bxy+cy^2+dx+ey+f.$ This includes more than rectilinear parabolas: it consists of all conic sections. The next uncertainty I have concerns the error model: should we assume the variances in $x$ and $y$ are equal or not? Can we assume the errors in the coordinates are independent or at least uncorrelated? Finally, without further information from the OP it is unclear how to quantify goodness of fit in the first place. – whuber Aug 9 '14 at 19:07

A total least squares solution for fitting a parabola in a 2D point set has not been investigated yet in the literature (or as far as I know). I would propose you to use the traditional geodetic way for solving such non-linear problems. This involves a linearization of the condition equations

i.e. $Yi$ + $VYi$ = $a$($Xi$ + $VXi$) + $b$($Xi$ + $VXi$)$^2$ + $c$) and iterations.

You can find this procedure under the name Gauss-Helmert Model.