Fit plane to 3D data using least squares

Question

I have some samples of data of the form $x,y$ and $z=f(x,y)$. I wish to fit a plane (i.e. $z = Ax + By + C$) to the data with the smallest mean square errors. I have found an "answer" in section 3 of this document, and several other locations too but the answers always end with variations of "now solve these equations and you can find $A$, $B$ and $C$"

I just about have the ability to solve these equations, but the process gets so messy that the likelihood of me making a mistake is quite high (and I need a guaranteed-correct answer). Surely someone has written out the full solution longhand somewhere - i.e. in the form $A=\dots$, $B=\dots$ and $C=\dots$ anyone know where this has been done?

EDIT: I see two answers already which leave out the last step as being trivial... and indeed they don't require any advanced maths... but you do need quite a big whiteboard to work it all out and the scope for making a mistake along the way is quite large. Indeed I noticed that assorted tutorials going through worked examples always have rather neat whole numbers. I find it hard to believe that there is nowhere online where someone has worked out the general case. I.e. find $A$, $B$ and $C$ in the following set of equations...

$$Ad+Be+Cf=g$$

$$Ah+Bi+Cj=j$$

$$Al+Bm+Cn=p$$

...where $d$ to $p$ are all constants.

Answer 1

If you don't feel confident with the resolution of a $3\times3$ system, work as follows:

• take the average of all equations,

$$\bar z=A\bar x+B\bar y+C$$

• subtract from all equations, giving

$$z_i-\bar z=A(x_i-\bar x)+B(y_i-\bar y)$$ or $$\hat z_i=A\hat x_i+B\hat y_i.$$

• solve the least squares system

$$\sum \hat z_i\hat x_i=A\sum \hat x_i^2+B\sum \hat x_i\hat y_i\\ \sum \hat z_i\hat y_i=A\sum \hat x_i\hat y_i+B\sum \hat y_i^2$$

which is $2\times 2$. This gives you $A$ and $B$.

By Cramer, $$A=\frac{\sum \hat z_i\hat x_i\sum \hat y_i^2-\sum \hat z_i\hat y_i\sum \hat x_i\hat y_i}{\sum \hat x_i^2\sum \hat y_i^2-\left(\sum \hat x_i\hat y_i\right)^2}\\ B=\frac{\sum \hat x_i^2\sum \hat z_i\hat x_i-\sum\hat z_i\hat x_i\sum \hat x_i\hat y_i}{\sum \hat x_i^2\sum \hat y_i^2-\left(\sum \hat x_i\hat y_i\right)^2}$$

• C is found from

$$C=\bar z-A\bar x-B\bar y.$$

---

For convenient evaluation of the sums, notice that

$$\sum \hat u_i\hat v_i=\sum u_iv_i-N\bar u\bar v,$$ i.e. $$\overline{\hat u\hat v}=\overline{uv}-\bar u\bar v.$$ Then

$$\begin{align}A&=\frac{(\overline{zx}-\bar z\bar x)(\overline{y^2}-\bar y^2)-(\overline{zy}-\bar z\bar y)(\overline{xy}-\bar x\bar y)}{(\overline{x^2}-\bar x^2)(\overline{y^2}-\bar y^2)-(\overline{xy}-\bar x\bar y)^2}\\ B&=\frac{(\overline{x^2}-\bar x^2)(\overline{zy}-\bar z\bar y)-(\overline{zx}-\bar z\bar x)(\overline{xy}-\bar x\bar y)}{(\overline{x^2}-\bar x^2)(\overline{y^2}-\bar y^2)-(\overline{xy}-\bar x\bar y)^2}\\ C&=\bar z-A\bar x-B\bar y.\end{align}$$

Answer 2

To do least square fit, you simply follow these three steps.

1.) Define the least square error function, for example, you could do

$ e(A, B, C) = \sum(Ax + By + C - z)^2 $

2.) Assuming you have no special constraint you want on $ A $, $ B $ and $ C $, so it is simply an unconstrained optimization problem. You differentiate the error function and set zero.

$ \frac{\partial e}{\partial A} = \sum 2(Ax + By + C - z)x = 0 $

$ \frac{\partial e}{\partial B} = \sum 2(Ax + By + C - z)y = 0 $

$ \frac{\partial e}{\partial C} = \sum 2(Ax + By + C - z) = 0$

3.) At this point, I guess your reference just ask you to solve for $ A, B, C $ and you are unsure about how to do that.

The key observation is that these are just linear equations!

Ley say, for example, that you have these 4 data points.

(12, 31, 27) (22, 32, 37) (13, 33, 17) (0, 0, 0)

You put in to $ \frac{\partial e}{\partial A} $ and see?

That gives

$ 2(12A + 31B + C - 27) + 2(22A + 32B + C - 37) + 2(13A + 33B + C - 17) + 2(0A + 0B + C - 0) $

$ = 94A + 192B + 7C - 162 = 0 $

The other equations can be simplified similarly. Then the rest is just solving 3 linear equations with 3 unknowns, any matrix method would work.

Answer 3

Say you have $n$ data sets:

$( x_0 , y_0 , z_0 ), ( x_1 , y_1 , z_1 ), ( x_2 , y_2 , z_2 ), ( x_3 , y_3 , z_3 ),...( x_n , y_n , z_n )$

Create a matrix $A\in R^{(n+1)x(n+1)}$ s.t.

$$ \begin{matrix} 1 & x_0 & y_0 \\ 1 & x_1 & y_1 \\ 1 & x_2 & y_2 \\ 1 & x_3 & y_3 \\ . & . & . \\ 1 & x_n & y_n \\ \end{matrix} $$

Afterwards create a column (result) vector $b\in R^{(n+1)x1}$ s.t.

$$ \begin{matrix} z_0 \\ z_1 \\ z_2 \\ z_3 \\ . \\ z_n \\ \end{matrix} $$

Then we get the equation $A*u = b$ where $u$ is your unknown vector $$ \begin{matrix} C \\ A \\ B \\ \end{matrix} $$

Final step is to obtain the normal equation $A^T(A*u) = A^T(b)$ which is a 3x3 system of linear equations:

$(A^TA)*u = A^T(b)$

I think the rest is trivial.