三维最小二乘平面

Question

在给定一组3D数据点的情况下，计算(x，y，z)空间中的最小二乘平面的算法是什么？换句话说，如果我有一堆像(1，2，3)，(4，5，6)，(7，8，9)等的点，如何计算最佳拟合平面f(x，y) = ax + by + c？从一组3D点中获取a、b和c的算法是什么？

Answer 1

如果您有n个数据点(xi，yi，zi)，则计算3x3对称矩阵A，其条目为：

sum_i x[i]*x[i],    sum_i x[i]*y[i],    sum_i x[i]
sum_i x[i]*y[i],    sum_i y[i]*y[i],    sum_i y[i]
sum_i x[i],         sum_i y[i],         n

还要计算3个元素的向量b：

{sum_i x[i]*z[i],   sum_i y[i]*z[i],    sum_i z[i]}

然后对给定的A和b求解Ax =b。解向量的三个分量是最小二乘拟合平面{a，b，c}的系数。

请注意，这是“普通最小二乘”拟合，只有当期望z是x和y的线性函数时才适用。如果你在寻找三维空间中的“最佳拟合平面”，你可能想要学习“几何”最小二乘。

另请注意，如果您的点在一条线上，则此操作将失败，就像您的示例点一样。

Answer 2

一个平面的方程式是: ax + by +c= z。所以用你所有的数据建立如下矩阵：

    x_0   y_0   1  
A = x_1   y_1   1  
          ... 
    x_n   y_n   1  

和

    a  
x = b  
    c

和

    z_0   
B = z_1   
    ...   
    z_n

换句话说: Ax = B，现在求解x，它是你的系数。但是由于(我假设)你有超过3个点，系统是过度确定的，所以你需要使用左边的伪逆。所以答案是：

a 
b = (A^T A)^-1 A^T B
c

下面是一些简单的Python代码和一个示例：

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np

N_POINTS = 10
TARGET_X_SLOPE = 2
TARGET_y_SLOPE = 3
TARGET_OFFSET  = 5
EXTENTS = 5
NOISE = 5

# create random data
xs = [np.random.uniform(2*EXTENTS)-EXTENTS for i in range(N_POINTS)]
ys = [np.random.uniform(2*EXTENTS)-EXTENTS for i in range(N_POINTS)]
zs = []
for i in range(N_POINTS):
    zs.append(xs[i]*TARGET_X_SLOPE + \
              ys[i]*TARGET_y_SLOPE + \
              TARGET_OFFSET + np.random.normal(scale=NOISE))

# plot raw data
plt.figure()
ax = plt.subplot(111, projection='3d')
ax.scatter(xs, ys, zs, color='b')

# do fit
tmp_A = []
tmp_b = []
for i in range(len(xs)):
    tmp_A.append([xs[i], ys[i], 1])
    tmp_b.append(zs[i])
b = np.matrix(tmp_b).T
A = np.matrix(tmp_A)
fit = (A.T * A).I * A.T * b
errors = b - A * fit
residual = np.linalg.norm(errors)

print "solution:"
print "%f x + %f y + %f = z" % (fit[0], fit[1], fit[2])
print "errors:"
print errors
print "residual:"
print residual

# plot plane
xlim = ax.get_xlim()
ylim = ax.get_ylim()
X,Y = np.meshgrid(np.arange(xlim[0], xlim[1]),
                  np.arange(ylim[0], ylim[1]))
Z = np.zeros(X.shape)
for r in range(X.shape[0]):
    for c in range(X.shape[1]):
        Z[r,c] = fit[0] * X[r,c] + fit[1] * Y[r,c] + fit[2]
ax.plot_wireframe(X,Y,Z, color='k')

ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
plt.show()

​

​

Answer 3

除非有人告诉我如何在这里输入方程式，否则让我写下您必须进行的最终计算：

首先，给定点r_i \n \R，i=1..N，计算所有点的质心：

r_G = \frac{\sum_{i=1}^N r_i}{N}

然后，计算法向量n，它与基向量r_G一起通过计算3x3矩阵A来定义平面：

A = \sum_{i=1}^N (r_i - r_G)(r_i - r_G)^T

有了这个矩阵，法向量n现在由A的特征向量给出，对应于A的最小特征值。

要了解特征向量/特征值对，请使用您选择的任何线性代数库。

此解基于Hermitian矩阵A的Rayleight-Ritz定理。