Best fit plane to 3D data: different results for two different methods

Question

(Although there are a number of questions regarding how to best fit a plane to some 3D data on SO, I couldn't find an answer for this issue.)

Given N (x, y, z) points, I need the best fit plane

a*x + b*y + c*z + d = 0

defined through the a, b, c, d coefficients that minimize the mean of the orthogonal distances from the points to the plane. The point-plane orthogonal distance (for a given (x0, y0, z0) point) is defined as:

d = |a*x0 + b*y0 + c*z0 + d|/sqrt(a^2 + b^2 + c^2)

I set up two methods (code below):

• Singular-Value Decomposition (source)
• Basin-Hopping minimization of the mean orthogonal distances

As I understand it, the SVD method should produce the exact best fit plane by minimizing the orthogonal distances analytically. What I find instead is that the BH method gives better results than the supposedly exact SVD method, even for a low number of BH runs.

By "better" I mean that the final mean orthogonal distance value is smaller with the BH method, than with the SVD method.

What am I missing here?

---

import numpy as np
import scipy.optimize as optimize


def perp_error(params, xyz):
    """
    Mean of the absolute values for the perpendicular distance of the
    'xyz' points, to the plane defined by the coefficients 'a,b,c,d' in
    'params'.
    """
    a, b, c, d = params
    x, y, z = xyz
    length = np.sqrt(a**2 + b**2 + c**2)
    return (np.abs(a * x + b * y + c * z + d) / length).mean()


def minPerpDist(x, y, z, N_min):
    """
    Basin-Hopping method, minimize mean absolute values of the
    orthogonal distances.
    """
    def unit_length(params):
        """
        Constrain the vector perpendicular to the plane to be of unit length.
        """
        a, b, c, d = params
        return a**2 + b**2 + c**2 - 1

    # Random initial guess for the a,b,c,d plane coefficients.
    initial_guess = np.random.uniform(-10., 10., 4)

    # Constrain the vector perpendicular to the plane to be of unit length.
    cons = ({'type': 'eq', 'fun': unit_length})
    min_kwargs = {"constraints": cons, "args": [x, y, z]}

    # Use Basin-Hopping to obtain the best fit coefficients.
    sol = optimize.basinhopping(
        perp_error, initial_guess, minimizer_kwargs=min_kwargs, niter=N_min)
    abcd = list(sol.x)

    return abcd


def SVD(X):
    """
    Singular value decomposition method.
    Source: https://gist.github.com/lambdalisue/7201028
    """
    # Find the average of points (centroid) along the columns
    C = np.average(X, axis=0)

    # Create CX vector (centroid to point) matrix
    CX = X - C
    # Singular value decomposition
    U, S, V = np.linalg.svd(CX)
    # The last row of V matrix indicate the eigenvectors of
    # smallest eigenvalues (singular values).
    N = V[-1]

    # Extract a, b, c, d coefficients.
    x0, y0, z0 = C
    a, b, c = N
    d = -(a * x0 + b * y0 + c * z0)

    return a, b, c, d


# Generate a random plane.
seed = np.random.randint(100000)
print("Seed: {}".format(seed))
np.random.seed(seed)
a, b, c, d = np.random.uniform(-10., 10., 4)
print("Orig abc(d=1): {:.3f} {:.3f} {:.3f}\n".format(a / d, b / d, c / d))

# Generate random (x, y, z) points.
N = 200
x, y = np.random.uniform(-5., 5., (2, N))
z = -(a * x + b * y + d) / c
# Add scatter in z.
z = z + np.random.uniform(-.2, .2, N)

# Solve using SVD method.
a, b, c, d = SVD(np.array([x, y, z]).T)
print("SVD  abc(d=1): {:.3f} {:.3f} {:.3f}".format(a / d, b / d, c / d))
# Orthogonal mean distance
print("Perp err: {:.5f}\n".format(perp_error((a, b, c, d), (x, y, z))))

# Solve using Basin-Hopping.
abcd = minPerpDist(x, y, z, 500)
a, b, c, d = abcd
print("BH   abc(d=1): {:.3f} {:.3f} {:.3f}".format(a / d, b / d, c / d))
print("Perp err: {:.5f}".format(perp_error(abcd, (x, y, z))))

Answer 1

I believe I found the reason for the discrepancy.

When I minimize the perpendicular distance of points to a plane using Basin-Hopping, I am using the absolute valued point-plane distance:

d_abs = |a*x0 + b*y0 + c*z0 + d| / sqrt(a^2 + b^2 + c^2)

The SVD method on the other hand, apparently minimizes the squared point-plane distance:

d_sqr = (a*x0 + b*y0 + c*z0 + d)^2 / (a^2 + b^2 + c^2)

If, in the code shared in the question, I use the squared distance in the perp_error() function instead of the absolute valued distance, both methods give the exact same answer.