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:
- There are many permuations of the "moments", finding them all is a hazzle.
- 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:
- Do a Least Square Error fit of a plane to the data: $ z = dx + ey + f$
- Subtract the plane z values from all measurements.
- 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?