Orthogonal Distance Regression Planes

Date: 07/30/2003 at 12:30:04
From: R.P.
Subject: Fitting a plane into a set of points

I have a set of data points that I have collected from an experiment. 
I want to fit a 3D plane (best-fit) into these points (the points are 
in the form (x1,y1,z1), (x2,y2,z2),...) in order to evaluate my 
results. I know doing this calculation is not going to be easy. 
However, I want to see the approach so i can get an idea. Even a link 
to a website that can help would be enough.

Date: 07/30/2003 at 13:12:45
From: Doctor George
Subject: Re: Fitting a plane into a set of points

Hi R.P.,

Thanks for writing to Doctor Math.

There are different kinds of best fit planes. One best fit plane 
minimizes the maximum distance from the points to the plane. Another 
minimizes the sum of squared distances to the plane.

Actually, it gets more complex than that. For the least squares plane, 
there is one kind in which the x and y values are fixed, and the 
measured error is in z alone. This is called a regression plane. For 
this plane the minimized distance is only in the z direction. There is 
also a orthogonal distance regression plane that minimizes the 
perpendicular distances to the plane. This is used when there is 
measurement error in all three coordinates.

Does this help you pin down what kind of best fit you are looking 
for? We can get into some details once we get this sorted out.

- Doctor George, The Math Forum
  http://mathforum.org/dr.math/ 

Date: 07/30/2003 at 13:22:38
From: R.P.
Subject: Fitting a plane into a set of points

Thanks for the quick reply. The type of best fit I am interested in 
is orthogonal distance regression planes.

Date: 07/30/2003 at 14:03:01
From: Doctor George
Subject: Re: Fitting a plane into a set of points

Hi R.P.,

Finding the orthogonal distance regression plane is an eigenvector 
problem. It is quite involved, as you suspected. The best solution 
utilizes the Singular Value Decomposition (SVD). I hope you have 
access to a good linear algebra software package.

Starting with the distance from a point to a plane, we wish to find 
a, b, c and d such that we minimize

  f(a,b,c,d) = Sum [|axi + byi + czi + d|^2 / (a^2 + b^2 + c^2)]

If we set the partial derivative with respect to d equal to zero, we 
can solve for d to get

        d = -(a*x0 + b*y0 + c*z0)

where (x0, y0, z0) is the centroid of the data. This means that the 
least squares plane contains the centroid. If we substitute it back 
into the equation for the plane we get

        a(x - x0) + b(y - y0) + c(z - z0) = 0

We can rewrite f(a,b,c,d) like this

  f(a,b,c) = Sum [|a(xi-x0) + b(yi-y0) + c(zi-z0)|^2 / (a^2+b^2+c^2)]

Now we are going to switch over to a matrix representation. I will 
be using both the upper and lower cases of some letters. I hope that 
does not cause confusion.

Let's define v and M such that

             T
            v  = [a  b  c]

                 --                             --
                 | x1 - x0    y1 - y0    z1 - z0 |
                 | x2 - x0    y2 - y0    z2 - z0 |
            M =  |    .          .          .    |
                 |    .          .          .    |
                 |    .          .          .    |
                 | xn - x0    yn - y0    zn - z0 |
                 --                             --

If you multiply the matrices out, you will see that f(a,b,c) becomes

               T T          T
      f(v) = (v M )(Mv) / (v v)

              T  T        T
           = v (M M)v / (v v)

                  T
Let's define A = M M, which when divided by the number of data points 
becomes the covariance matrix of the data.

f(v) is called a Rayleigh Quotient. It is minimized by the eigenvector 
of A that corresponds to its smallest eigenvalue.

We could compute the eigenvectors of A, but this is not needed. The 
SVD of M is
                    T
             M = USV 

where S is a diagonal matrix containing the singular values of M, the 
columns of V are its singular vectors, and U is an orthogonal matrix.

                    T
Now            A = M M

                       T T    T
                 = (USV ) (USV )

                       T T      T
                 = (V S U ) (USV )

                      2 T 
                 = V S V

This decomposition of A diagonalizes the matrix and provides an 
eigenvector decomposition. It means that the eigenvalues of A are the 
squares of the singular values of M, and the eigenvectors of A are the 
singular vectors of M.

To conclude, the orthogonal least squares 3D plane contains the 
centroid of the data, and its normal vector is the singular vector of 
M corresponding to its smallest singular value.

Does this give you enough to go on? Write again if you need more 
help.

- Doctor George, The Math Forum
  http://mathforum.org/dr.math/ 

Date: 06/29/2005 at 18:41:25
From: Russ
Subject: Best Fit Plane for a set of points

I have read your excellent answer.  In it you state that, "If we 
set the partial derivative with respect to d equal to zero, we can 
solve for d to get d = -(ax0 + by0 + cz0) where (x0,y0,z0) is the
centroid of the data.  This means that the least squares plane contains
the centroid."

When I take the partial derivative I get d = -(axi + byi + czi) but 
I can not see why this leads to the conclusion that this is the 
centroid. 
               
It seems logical that this would be true, however I would like to 
know why does the best fit plane HAVE to contain the centroid?  I am 
doing some work with best fit planes and have even solved some 
problems using your article.  However my coworkers are very skeptical 
of what I am doing and I'd like to convince them this is in fact true.  

Date: 06/30/2005 at 08:04:41
From: Doctor George
Subject: Re: Best Fit Plane for a set of points

Hi Russ,

Thanks for writing to Doctor Math. I'm glad that you found my article
to be helpful.

   f(a,b,c,d) = Sum [|axi + byi + czi + d|^2 / (a^2 + b^2 + c^2)]

I'll use p for partial derivative since the text character set does
not contain the usual symbol for partial derivative.

   pf/pd = 2 Sum [(axi + byi + czi + d) / (a^2 + b^2 + c^2)] = 0

leads to

           Sum (axi + byi + czi + d) = 0

           a Sum xi + b Sum yi + c Sum zi + Nd = 0

where N is the number of points. Now we continue,

           Nd = -(a Sum xi + b Sum yi + c Sum zi)

           d = -(a Sum xi + b Sum yi + c Sum zi) / N

           d = -[a (Sum xi)/N + b (Sum yi)/N + c (Sum zi)/N]

           d = -(a*x0 + b*y0 + c*z0)

Does that make sense? Write again if you need more help.

- Doctor George, The Math Forum
  http://mathforum.org/dr.math/