# Least Squares “analytic expression” for fitting a 2D quadratic function to measurements

I have n scattered elevation measurements: $\{x_i,y_i,z_i\}_{i=1..n}$ that I want to fit a quadratic function to: $z = ax^2 + by^2 + cxy + dx + ey + f$.

The problem can be written as a vector equation: $$\vec{z}= \bf{X}\vec{\beta}$$ where the 'ith row of X is: $$Xi = [x^2, x, xy, y^2, y, 1]$$ and $$\vec{\beta} =[a, b, c, d, e, f]^T$$

and solved by QR decomposition of X.

Instead I am considering another approach (the approach one typically uses for y = ax +b):

$E = \Sigma (ax^2 + by^2 + cxy + dx + ey + f - z)^2$

$\partial E / \partial a=\Sigma 2(ax^2 + by^2 + cxy + dx + ey + f - z)x^2=0$

$a\overline{x^4} + b\overline{x^2y^2} + c\overline{x^3y} + d\overline{x^3} + e\overline{x^2y} + f\overline{x^2} - \overline{x^2z} = 0$

where the overline indicates average.

Likewise demanding that the derivative of E wrt to b, c, d, e and f are zero results in a linear equation for $\beta$ where the matrix is a 6x6 matrix of average "moments". This can then be solved by Cramer's rule.

There are some practical problems with this approach:

1. There are many permuations of the "moments", finding them all is a hazzle.
2. Cramer's rule requires that I find the determinant of a 6x6 matrix. The analytical expression for this is extremely long and difficult to get right, and I want to avoid using decompositions.

Therefore I was considering simplifying the approach by doing it in two stages:

1. Do a Least Square Error fit of a plane to the data: $z = dx + ey + f$
2. Subtract the plane z values from all measurements.
3. Fit a quadratic function to the data: $z = ax^2 + by^2 + cxy$

Will the simplified approach give the same result as my "complete" approach?

Are there any issues with numerical stability that I must be aware of?

-
(1) By the way, the $6\times6$ matrix you get from your intermediate approach is just (some multiple of) the matrix $\mathbf X^T\mathbf X$ in the normal equations, $\mathbf X^T\mathbf X\mathbf \beta=\mathbf X^T\mathbf z$. (2) No, your simplified approach will not work. Suppose the data fits $z=x^2$ but all the $x_i$'s are positive. A plane fit will incorrectly give you positive coefficient for the linear term in $x$. – Rahul Aug 16 '13 at 11:31
Do all of this symbolically, in Mathematica. This will give you a closed-form solution, but it will be a gigantic mess of course. Then use the CForm command to generate code. – bubba Aug 16 '13 at 11:38
@Rahul: thanks! How about just finding f as the average of z first? This will make the problem slightly easier (5x5). – Andy Aug 16 '13 at 12:15
@bubba: thanks! I was not aware that Mathematica had this function. I tried wolfram alpha but it complained about long too long input text. I am therefore considering buying a Mathematica license. Do you have any experience with wolfram alpha professional? Maybe that would suffice. – Andy Aug 16 '13 at 12:18
What I have is a very old Mathematica license (version 6, I think). I'm not familiar with Wolfram's newer products. Their products are expensive, of course, so think this through before you spend the money. I have done similar things in the past, but your particular problem may be different in some way that I'm not seeing. – bubba Aug 17 '13 at 0:03