Take the 2-minute tour ×
Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It's 100% free, no registration required.

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?

share|improve this question
1  
(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
1  
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
1  
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
show 1 more comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Browse other questions tagged or ask your own question.